1 Introduction
Many turbulent thermal convection phenomena in nature and technology are present in low-Prandtl-number fluids for which the kinematic viscosity of the fluid,
$\unicode[STIX]{x1D708}$
, is much smaller than the thermal diffusivity of the temperature field,
$\unicode[STIX]{x1D705}$
, or in other words, for which the ratio of both, the Prandtl number,
$$\begin{eqnarray}Pr=\frac{\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D705}},\end{eqnarray}$$
is much smaller than one. Applications range from turbulent convection in the Sun (Hanasoge et al.
Reference Hanasoge, Gizon and Sreenivasan2015) or the liquid-metal core of the Earth (King & Aurnou Reference King and Aurnou2013) to nuclear engineering (Grötzbach Reference Grötzbach2013) or liquid-metal batteries (Kelley & Sadoway Reference Kelley and Sadoway2014). In all these examples
$Pr\lesssim 2\times 10^{-2}$
, and in the solar convection case it could be smaller than
$Pr\sim 10^{-4}$
. Convective turbulence in these examples appears mostly in combination with other physical processes, such as rotation, radiation or the generation of magnetic fields. Rayleigh–Bénard convection (RBC), a fluid flow in a layer which is cooled from above and heated from below, is the simplest example of a turbulent convective flow.
The record of experimental and numerical data for low-Prandtl-number convection is much smaller than for convection in air or water. Laboratory experiments are challenging in several ways. For Prandtl numbers
$Pr<0.2$
, opaque liquid metals such as liquid mercury or gallium at
$Pr\simeq 0.025$
(Rossby Reference Rossby1969; Cioni et al.
Reference Cioni, Ciliberto and Sommeria1996, Reference Cioni, Ciliberto and Sommeria1997a
; Takeshita et al.
Reference Takeshita, Segawa, Glazier and Sano1996; Glazier et al.
Reference Glazier, Segawa, Naert and Sano1999; Aurnou & Olson Reference Aurnou and Olson2001; King & Aurnou Reference King and Aurnou2013) or liquid sodium at
$Pr\simeq 0.005$
(Cioni et al.
Reference Cioni, Horanyi, Krebs and Müller1997b
; Horanyi et al.
Reference Horanyi, Krebs and Müller1999; Frick et al.
Reference Frick, Khalilov, Kolesnichenko, Mamykin, Pakholkov, Pavlinov and Rogozhkin2015) have to be used. Additionally, experiments in sodium require a working temperature of more than
$100\,^{\circ }$
C. An optical access to the flow by laser-imaging techniques is impossible. The heating and cooling at the bottom and top walls are usually established with copper plates. In liquid sodium experiments, the thermal conductivity of the plates is thus not significantly higher than that of the working fluid. This introduces significant variations of the temperature at the plates and can alter the magnitude of the heat transfer since the boundary conditions will be a mixture of prescribed temperature and prescribed normal temperature derivative, known as Robin boundary conditions. For an experimental analysis we refer to Horanyi et al. (Reference Horanyi, Krebs and Müller1999). In low-Prandtl-number convection, the thermal boundary layer at the solid heating and cooling plates is much thicker than the velocity boundary layer. This has important implications for the turbulent heat and momentum transfer which is one focus of the research in turbulent convection (e.g. in Grossmann & Lohse Reference Grossmann and Lohse2011; Chillà & Schumacher Reference Chillà and Schumacher2012; Scheel & Schumacher Reference Scheel and Schumacher2014). The global heat transfer is quantified by the dimensionless Nusselt number,
$Nu$
, while the global momentum transfer is measured by the Reynolds number,
$Re$
. Detailed definitions of both parameters are given in § 2.1. Both the Nusselt and Reynolds numbers are functions of the Prandtl number
$Pr$
and a second dimensionless parameter which quantifies the thermal driving of the turbulence, namely the Rayleigh number
$Ra$
. Different scaling laws
$Nu(Ra)$
have been reported in different convection experiments for mercury and gallium, i.e. at fixed Prandtl number. In liquid sodium, only one experimental data record by Horanyi et al. (Reference Horanyi, Krebs and Müller1999) exists which is similar to the parameters in this paper. This discussion implies that direct numerical simulations (DNS) are important for gaining access to the full three-dimensional turbulent fields in low-
$Pr$
convection and to assure the accuracy of the boundary conditions.
In the present work, we study the local and global statistical properties of turbulent convection flows at very low Prandtl numbers as present in liquid mercury and liquid sodium. We analyse a series of very high-resolution DNS obtained for moderate Rayleigh numbers in a cylindrical cell with an aspect ratio of one. We obtain global scaling laws for the turbulent heat and momentum transfer and compare our findings with available data as well as the scaling theory of Grossmann & Lohse (Reference Grossmann and Lohse2001) and a recent update by Stevens et al. (Reference Stevens, van der Poel, Grossmann and Lohse2013). In addition, we analyse the turbulence statistics for a series of three DNS runs which have been conducted at different Rayleigh and Prandtl numbers but with the same Grashof number,
$Gr=Ra/Pr$
, for all three runs (for definition see § 2.2). This perspective (Schumacher et al.
Reference Schumacher, Götzfried and Scheel2015) leaves the dimensionless Navier–Stokes equations unchanged and thus provides a more consistent comparison of convective flow at different Prandtl numbers from the point of view of fluid turbulence.
A second focus of the present work is a detailed local statistical analysis of the dynamics in the boundary layers of the velocity and temperature fields. While the heat transport is reduced in low-
$Pr$
convection, production of vorticity and shear are strongly enhanced, both in the bulk and in the boundary layers. This is the reason why simulations of turbulent convection at very low Prandtl numbers become very demanding, as is reported for homogeneous turbulent flows (Mishra & Verma Reference Mishra and Verma2010), in planar turbulent flows between parallel free-slip or no-slip walls (Kerr & Herring Reference Kerr and Herring2000; Breuer et al.
Reference Breuer, Wessling, Schmalzl and Hansen2004; Verma et al.
Reference Verma, Mishra, Pandey and Paul2012; Petschel et al.
Reference Petschel, Stellmach, Wilczek, Lülff and Hansen2013) or in closed cylindrical cells (Camussi & Verzicco Reference Camussi and Verzicco1998; van der Poel et al.
Reference van der Poel, Stevens and Lohse2013; Schumacher et al.
Reference Schumacher, Götzfried and Scheel2015). We focus our interest on the boundary layer dynamics and determine that the enhanced intermittency in the bulk is connected to a stronger intermittent behaviour in the boundary layer, in particular for the velocity field. But, while the vorticity and shear are enhanced in low-Prandtl-number convection, we also find that the scaling of many diagnostic quantities is weaker. We postulate that this is due to the reduced heat transport and weaker dependence of the Nusselt number on the Rayleigh number.
The outline of the paper is as follows. In the next section, we present the model equations, the numerical method and the data sets. Section 3 discusses the scaling laws for global heat and momentum transfer. Section 4 turns to the statistics at constant Grashof number and is followed by a detailed boundary-layer analysis in § 5. Finally we summarize our findings and give a brief outlook for future work.
2 Numerical model
2.1 Equations of motion
We solve the three-dimensional equations of motion in the Boussinesq approximation. The equations are made dimensionless by using height of the cell
$H$
, the free-fall velocity
$U_{f}=\sqrt{g\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}TH}$
and the imposed temperature difference
$\unicode[STIX]{x0394}T$
. The equations contain the three control parameters: the Rayleigh number
$Ra$
, the Prandtl number
$Pr$
and the aspect ratio
$\unicode[STIX]{x1D6E4}=D/H$
with the cell diameter
$D$
. The equations are given by
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{u}=-\unicode[STIX]{x1D735}p+\sqrt{\frac{Pr}{Ra}}\unicode[STIX]{x1D6FB}^{2}\boldsymbol{u}+T\boldsymbol{e}_{z}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}t}+(\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735})T=\frac{1}{\sqrt{RaPr}}\unicode[STIX]{x1D6FB}^{2}T, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}Ra=\frac{g\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}TH^{3}}{\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}},\quad Pr=\frac{\unicode[STIX]{x1D708}}{\unicode[STIX]{x1D705}}.\end{eqnarray}$$
The free-fall Reynolds number
$Re_{f}$
is given by
$$\begin{eqnarray}Re_{f}=\frac{U_{f}H}{\unicode[STIX]{x1D708}}=\sqrt{\frac{Ra}{Pr}}.\end{eqnarray}$$
The variable
$g$
stands for the acceleration due to gravity,
$\unicode[STIX]{x1D6FC}$
is the thermal expansion coefficient,
$\unicode[STIX]{x1D708}$
is the kinematic viscosity and
$\unicode[STIX]{x1D705}$
is the thermal diffusivity. We use an aspect ratio of
$\unicode[STIX]{x1D6E4}=1$
here. No-slip boundary conditions for the fluid (
$\boldsymbol{u}=0$
) are applied at the walls. The side walls are thermally insulated (
$\unicode[STIX]{x2202}T/\unicode[STIX]{x2202}\boldsymbol{n}=0$
) and the top and bottom plates are held at constant dimensionless temperatures
$T=0$
and 1, respectively. In response to the input parameters
$Ra$
,
$Pr$
and
$\unicode[STIX]{x1D6E4}$
, turbulent heat and momentum fluxes are established. The turbulent heat transport is determined by the Nusselt number which is defined as
$$\begin{eqnarray}Nu=1+\sqrt{RaPr}\langle u_{z}T\rangle _{V,t}.\end{eqnarray}$$
The turbulent momentum transport is expressed by the Reynolds number which is defined as
$$\begin{eqnarray}Re=u_{rms,V}\sqrt{\frac{Ra}{Pr}}=u_{rms,V}Re_{f}\quad \text{with}~u_{rms,V}=\sqrt{\langle u_{i}^{2}\rangle _{V,t}}.\end{eqnarray}$$
In both definitions
$\langle \cdot \rangle _{V,t}$
stands for a combined volume and time average.
The equations are numerically solved by the Nek5000 spectral element method package which has been adapted to our problem. The code employs second-order time stepping using a backward difference formula. The whole set of equations is transformed into a weak formulation and discretized with a particular choice of spectral basis functions (Fischer Reference Fischer1997; Deville et al. Reference Deville, Fischer and Mund2002). The resulting linear, symmetric Stokes problem is solved implicitly. This system is split, decoupling the viscous and pressure steps into independent symmetric positive definite subproblems which are solved either by Jacobi (viscous) or multilevel Schwartz (pressure) preconditioned conjugate gradient iteration. Fast parallel solvers based on direct projection or a more scalable algebraic multigrid are used for the coarse-grid solve that is part of the pressure preconditioner. For further numerical details and comprehensive tests of the sufficient spectral resolution, we refer to detailed investigations in Scheel et al. (Reference Scheel, Emran and Schumacher2013).
Prandtl number–Rayleigh number parameter space. The solid diagonal line corresponds to data at constant Grashof number. The different working fluids are indicated in the legend. The data for air are from Scheel & Schumacher (Reference Scheel and Schumacher2014).
Parameters of the different spectral element simulations. We show the Rayleigh number
$Ra$
, the Prandtl number
$Pr$
, the total number of spectral elements
$N_{e}$
, the polynomial order
$N$
and the number of grid points inside the velocity boundary layer,
$N_{bl}$
. The latter is calculated by (5.5). We also show the Nusselt number
$Nu$
(see (2.6)), the thermal boundary layer thickness
$\unicode[STIX]{x1D6FF}_{T}$
(see (5.7)), the root-mean-square velocity
$u_{rms,V}$
and the Reynolds number
$Re$
(see (2.7)). The error bars have been obtained by evaluating the results over the first and second halves of the corresponding data set separately and taking the difference of both results.
2.2 Data sets
The
$Pr{-}Ra$
parameter space covered by the data in table 1 is shown schematically in figure 1. We have a fairly complete data set for thermal convection in air at
$Pr=0.7$
(see Scheel & Schumacher Reference Scheel and Schumacher2014), and have added new data at the very low Prandtl numbers for liquid mercury,
$Pr=0.021$
, and liquid sodium,
$Pr=0.005$
. We have chosen the Rayleigh numbers in data sets 1, 4 and 8 of table 1 to keep the Grashof number
$Gr$
, defined as
$$\begin{eqnarray}Gr=\frac{g\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}TH^{3}}{\unicode[STIX]{x1D708}^{2}}=\frac{Ra}{Pr}=Re_{f}^{2},\end{eqnarray}$$
constant at
$4.76\times 10^{8}$
. Consequently,
$Gr=Re^{2}/u_{rms,V}^{2}$
with (2.7).
Snapshots of horizontal cross-sections at mid-plane of temperature field (a,c) and velocity magnitude (b,d) for
$Ra=10^{8}$
and
$Pr=0.021$
(a,b) and
$Pr=0.7$
(c,d). The data for
$Ra=10^{8}$
and
$Pr=0.7$
have been taken from Scheel & Schumacher (Reference Scheel and Schumacher2014).
Snapshots of horizontal cross-sections at
$\unicode[STIX]{x1D6FF}_{T}/2$
of temperature field (a,c) and velocity magnitude (b,d) for
$Ra=10^{8}$
and
$Pr=0.021$
(a,b) and
$Pr=0.7$
(c,d). The data for
$Ra=10^{8}$
and
$Pr=0.7$
have been taken from Scheel & Schumacher (Reference Scheel and Schumacher2014).
A comparison of the instantaneous data for the same Rayleigh number of
$Ra=10^{8}$
and
$Pr=0.021$
(a,b) versus
$Pr=0.7$
(c,d) is shown in figures 2 and 3. Colour density plots for both the temperature (a,c) and velocity magnitude (b,d) are shown for a representative snapshot. A horizontal cut through the cell (
$z=1/2$
) and inside the thermal boundary layer (
$z=\unicode[STIX]{x1D6FF}_{T}/2$
) are shown in figures 2 and 3, respectively. We see in both cases that the temperature is much more diffuse in the lower
$Pr$
case, but the velocity is higher in magnitude and exhibits finer structure and larger fluctuations. This combination of a diffuse temperature field coupled with a stronger fluctuating velocity field give rise to rather different dynamics.
3 Global turbulent heat and momentum transfer
In figure 4 we plot the Nusselt number as a function of Rayleigh number for all three Prandtl numbers (circles for
$Pr=0.7$
, stars for
$Pr=0.021$
and diamonds for
$Pr=0.005$
). We also added the data for
$Pr=6$
which we also have for
$Ra=1\times 10^{7}$
(see Scheel et al.
Reference Scheel, Emran and Schumacher2013). As expected there is a decrease in heat transport for a given Rayleigh number as the Prandtl number decreases. We find the scaling exponent decreases from
$0.29\pm 0.01$
for
$Pr=0.7$
to
$0.26\pm 0.01$
for
$Pr=0.021$
and to
$0.265\pm 0.01$
for
$Pr=0.005$
(see table 2). It should be stressed that the present data base for the low-Prandtl-number simulations is still rather small.
Nusselt number versus Rayleigh number for
$Pr=0.005$
(diamonds),
$Pr=0.021$
(stars),
$Pr=0.7$
(open circles) and
$Pr=6$
(filled square). The lines are fits to the data and are given in table 2. The inset replots all data compensated by
$Ra^{0.29}$
.
Reynolds number versus Rayleigh number for
$Pr=0.005$
(diamonds),
$Pr=0.021$
(stars) and
$Pr=0.7$
(open circles). The lines are fits to the data and are given in table 2. The inset replots all data compensated by
$Ra^{0.49}$
.
Nusselt number versus Prandtl number for three different Rayleigh numbers which are indicated in the legend. The solid lines are the fit to the Grossmann–Lohse equations as given in Stevens et al. (Reference Stevens, van der Poel, Grossmann and Lohse2013).
The scaling results are compared with previous published results in table 3, where the coefficients are given by
$$\begin{eqnarray}Nu=\unicode[STIX]{x1D6FC}Ra^{\unicode[STIX]{x1D6FD}}.\end{eqnarray}$$
For
$Pr=0.021$
, our exponent
$\unicode[STIX]{x1D6FD}$
and prefactor
$\unicode[STIX]{x1D6FC}$
agree within our uncertainty with experiments by Rossby (Reference Rossby1969), Cioni et al. (Reference Cioni, Ciliberto and Sommeria1996), Takeshita et al. (Reference Takeshita, Segawa, Glazier and Sano1996), King & Aurnou (Reference King and Aurnou2013) and simulations by Camussi & Verzicco (Reference Camussi and Verzicco1998) and Kerr & Herring (Reference Kerr and Herring2000). However, our exponent is on the high end of the range and our prefactor is on the low end of the range in table 3.
To our knowledge, we can compare our results for
$Pr=0.005$
to only one other convection experiment in liquid sodium which has been conducted at the Forschungszentrum Karlsruhe by Horanyi et al. (Reference Horanyi, Krebs and Müller1999) and is shown in table 3. Very recent experiments by Frick et al. (Reference Frick, Khalilov, Kolesnichenko, Mamykin, Pakholkov, Pavlinov and Rogozhkin2015) have been conducted in cylindrical cells with a very low aspect ratio of
$\unicode[STIX]{x1D6E4}=0.2$
. Such a geometry will generate a complex large-scale circulation (LSC) pattern that consists of several rolls on top of each other rather than a single roll as for
$\unicode[STIX]{x1D6E4}\sim 1$
. Therefore, these experiments are not used for comparison. Our Nusselt number value at
$Ra=1\times 10^{6}$
is
$Nu=4.45$
(see table 1) which is higher than the
$Nu\simeq 4.0$
which is reported by Horanyi et al. (Reference Horanyi, Krebs and Müller1999) although our results agree well for both the exponent and prefactor, within our numerical error. The conclusion which we can draw from the limited amount of data sets is that the scaling exponents
$\unicode[STIX]{x1D6FD}$
we find from our DNS decrease as
$Pr$
decreases, and agree with experiments and simulations when the uncertainty is taken into account. However, our measured exponents are 5 % larger than the one-fourth scaling reported in some experiments, (such as King & Aurnou Reference King and Aurnou2013), and 30 % larger than the one-fifth scaling predicted theoretically by Grossmann & Lohse (Reference Grossmann and Lohse2000). However, the difficulty of maintaining isothermal boundary conditions for liquid metals, as noted in Horanyi et al. (Reference Horanyi, Krebs and Müller1999), should be taken into consideration for the experiments. It was found numerically by Verzicco & Sreenivasan (Reference Verzicco and Sreenivasan2008) that when constant heat flux was used as the plate boundary conditions instead of isothermal, the Nusselt number for the largest Rayleigh numbers was reduced. This was attributed to the inhibition of thermal plume growth with constant heat flux. This could account for the smaller experimentally measured scaling exponents.
Comparison of scaling coefficients for Nusselt number with Rayleigh number (see (3.1)). We only selected cases which overlapped with our data range. The asterisks indicate numerical simulations; the others are experiments.
We then investigate the dependence of Reynolds number
$Re$
on Rayleigh number in figure 5 using (2.7) to compute our Reynolds numbers. We find that the Reynolds number increases rather significantly as Prandtl number decreases, which is a reflection of the much stronger velocity field. The scaling exponent is similar, but has decreased with Prandtl number from
$0.49\pm 0.01$
for
$Pr=0.7$
to
$0.44\pm 0.01$
for
$Pr=0.021$
, but then returns to
$0.49\pm 0.01$
for
$Pr=0.005$
. Experimentally Takeshita et al. (Reference Takeshita, Segawa, Glazier and Sano1996) found
$Re=6.24Ra^{0.46\pm 0.02}$
for
$Pr=0.025$
which agrees with our exponent and prefactor. In the simulations done by Verzicco & Camussi (Reference Verzicco and Camussi1999) a scaling law
$Re\simeq Ra^{0.53}$
is found which does not agree with our results.
In figure 6 we show the Nusselt number as a function of Prandtl number for three Rayleigh numbers. Also included as solid lines are fits to the scaling theory by Grossmann & Lohse (Reference Grossmann and Lohse2001) with the updated prefactors from Stevens et al. (Reference Stevens, van der Poel, Grossmann and Lohse2013). The fits are very good, especially for
$Pr=0.7$
and 6. However we do see a consistent disagreement for the lower Prandtl numbers, with the Grossmann–Lohse theory slightly underestimating the value for the Nusselt number. This was also seen in the simulations by van der Poel et al. (Reference van der Poel, Stevens and Lohse2013).
4 Turbulence statistics for constant Grashof number
In Schumacher et al. (Reference Schumacher, Götzfried and Scheel2015), we emphasized that a change in the parameter variation from constant
$Ra$
to constant
$Gr$
sheds a different light on the Prandtl number dependence in turbulent convection. Studies at constant Grashof number require a simultaneous variation of both the Rayleigh and Prandtl numbers as seen in figure 1. In this case, the momentum equation (2.2) remains unchanged. The explicit Prandtl number dependence appears only in the advection–diffusion equation (2.3). There is however an indirect Prandtl number dependence in the momentum equation as well. This is manifest in the buoyancy term. A more diffusive temperature field injects kinetic energy at a larger scale into the convection flow, as demonstrated in Schumacher et al. (Reference Schumacher, Götzfried and Scheel2015).
In figures 7 and 8 we compare root-mean-square profiles for different Prandtl numbers but fixed Grashof number for runs 1, 4 and 8. The root-mean-square values are defined as follows:
$$\begin{eqnarray}u_{rms}(z)=\sqrt{\langle u_{i}^{2}\rangle _{A,t}},\quad w_{rms}(z)=\sqrt{\langle u_{z}^{2}\rangle _{A,t}},\quad \unicode[STIX]{x1D703}_{rms}(z)=\sqrt{\langle (T-\langle T\rangle _{A,t})^{2}\rangle _{A,t}}.\end{eqnarray}$$
Figure 7(a) displays the expected increase of the velocity fluctuations, and thus of the turbulent kinetic energy, with decreasing Prandtl number. The vertical velocity fluctuations,
$w_{rms}$
yield a significant contribution to
$u_{rms}$
at all heights. As visible in the rescaled plot of figure 7(c), the profiles of
$w_{rms}$
always monotonically increase towards the cell centre plane at
$z=0.5$
. They form an increasingly wider plateau around the centre plane when the Prandtl number increases. Furthermore, in figure 7(b), for the lowest Prandtl number case, we detect that the maximum of the profiles of
$u_{rms}$
is found at
$z=0.5$
. This is different from the case for
$Pr=0.7$
where the fluctuations obey a maximum at
$z\approx \unicode[STIX]{x1D6FF}_{T}$
which can be attributed to the finer thermal plumes which detach from both boundary layers at higher
$Pr$
.
Root-mean-square velocity profiles for various Prandtl numbers. (a) The total root-mean-square velocity,
$u_{rms}$
, and the vertical root-mean-square velocity,
$w_{rms}$
, are given by (4.1). The thermal boundary layer thicknesses
$\unicode[STIX]{x1D6FF}_{T}$
obtained from (5.7) are also shown. All profiles are additionally averaged over the lower and upper halves of the cell. (b) Profiles of
$u_{rms}$
collapsed with respect to their amplitude in the mid-plane. (c) Profiles of
$w_{rms}$
collapsed with respect to their amplitude in the mid-plane. All data sets are for
$Gr=4.76\times 10^{8}$
.
Temperature field statistics for various Prandtl numbers. (a) The time-averaged temperature deviation
$\unicode[STIX]{x1D703}_{rms}$
is given by (4.1). The thermal boundary layer thicknesses
$\unicode[STIX]{x1D6FF}_{T}$
obtained from (5.7) are also shown. All profiles are additionally averaged over the lower and upper halves of the cell. (b) Scaled mean temperature profiles
$\unicode[STIX]{x1D6E9}(y)$
versus
$y=z/\unicode[STIX]{x1D6FF}_{T}$
(see (4.2)). All data sets are for
$Gr=4.76\times 10^{8}$
.
Energy dissipation rate profiles for various Prandtl numbers. The kinetic energy dissipation rate
$\unicode[STIX]{x1D716}$
and thermal dissipation rate
$\unicode[STIX]{x1D716}_{T}$
are given by area–time averages of (4.4) and (4.3), respectively. All profiles are additionally averaged over the lower and upper halves of the cell. All data sets are for
$Gr=4.76\times 10^{8}$
.
The profiles of the root-mean-square of the temperature fluctuations are shown in figure 8(a). All three profiles obey a local maximum that is close to the thermal boundary layer thickness. The maximum amplitude of the root-mean-square temperature fluctuations is of approximately the same magnitude for the two larger Prandtl numbers and drops significantly for the smallest accessed Prandtl number of
$Pr=0.005$
. For the latter case, the temperature field is so diffusive such that a detachment of several separate thermal plumes is not observed anymore. Rather one pronounced hot upwelling on one side of the cell is found in combination with a pronounced cold downwelling on the diametral side. In other words, the LSC is driven by a single pair of plumes which causes the drop in fluctuations.
Figure 8(b) shows the rescaled mean temperature profiles for the same data in the vicinity of the bottom plate. Following Shishkina & Thess (Reference Shishkina and Thess2009) we define
$$\begin{eqnarray}\unicode[STIX]{x1D6E9}(y)=2-2\langle T(z)\rangle _{A,t}\quad \text{with}~y=2Nu\,z=z/\unicode[STIX]{x1D6FF}_{T}.\end{eqnarray}$$
The dotted line indicates the slope of one of the profiles close to the wall. We see that the profiles for the low Prandtl numbers start to differ from the one for convection in air further away from the wall. The data for the lowest Prandtl number display an increase of
$\unicode[STIX]{x1D6E9}$
above one which indicates a positive slope of the unscaled mean temperature profile in the bulk. This behaviour is most probably attributed to the small Rayleigh number that was used in the sodium case.
The vertical profiles of the plane–time averages of kinetic energy dissipation and thermal dissipation rate fields are plotted in figure 9. They are denoted by
$\unicode[STIX]{x1D700}(z)=\langle \unicode[STIX]{x1D700}\rangle _{A,t}$
and
$\unicode[STIX]{x1D700}_{T}(z)=\langle \unicode[STIX]{x1D700}_{T}\rangle _{A,t}$
, respectively. The thermal dissipation rate field is given by
$$\begin{eqnarray}\unicode[STIX]{x1D700}_{T}(\boldsymbol{x},t)=\frac{1}{\sqrt{Gr}\,Pr}\left(\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}x_{j}}\right)^{2},\end{eqnarray}$$
and the kinetic energy dissipation rate field is defined as
$$\begin{eqnarray}\unicode[STIX]{x1D700}(\boldsymbol{x},t)=\frac{1}{2\sqrt{Gr}}\left(\frac{\unicode[STIX]{x2202}u_{i}}{\unicode[STIX]{x2202}x_{j}}+\frac{\unicode[STIX]{x2202}u_{j}}{\unicode[STIX]{x2202}x_{i}}\right)^{2}.\end{eqnarray}$$
The profiles imply that both dissipation rates increase in magnitude as the Prandtl number decreases. While the thermal dissipation rate grows as a result of the enhanced thermal diffusivity, the kinetic energy dissipation is increased as a result of the enhanced small-scale intermittency. The latter result has to be a consequence of enhanced local velocity derivatives since the prefactor in (4.4) is unchanged. The ratios
$Q_{T}=\unicode[STIX]{x1D700}_{T}(z=0)/\unicode[STIX]{x1D700}_{T}(z=0.5)$
and
$Q_{v}=\unicode[STIX]{x1D700}(z=0)/\unicode[STIX]{x1D700}(z=0.5)$
give
$Q_{T}\simeq 840,100,60$
and
$Q_{v}\simeq 19,16,19$
, respectively, for Prandtl numbers
$Pr=0.7,0.021,0.005$
. While the ratio of kinetic energy dissipation between boundary layer and bulk remains more or less unchanged, the dominance of thermal dissipation in the boundary layers for large Prandtl numbers is significantly diminished.
Mean kinetic energy dissipation rate and rescaled local enstrophy profiles for various Prandtl numbers. All profiles are additionally averaged over the lower and upper halves of the cell. All data sets are for
$Gr=4.76\times 10^{8}$
.
The local enstrophy can be defined by means of the vorticity field
$\unicode[STIX]{x1D714}_{i}=\unicode[STIX]{x1D716}_{ijk}\unicode[STIX]{x2202}_{j}u_{k}$
with
$i,j,k=1,2,3$
and is given by
$$\begin{eqnarray}\unicode[STIX]{x1D6FA}(\boldsymbol{x},t)=\unicode[STIX]{x1D714}_{i}(\boldsymbol{x},t)\unicode[STIX]{x1D714}_{i}(\boldsymbol{x},t).\end{eqnarray}$$
In homogeneous isotropic turbulence, the ensemble averages of the mean dissipation and the local enstrophy are connected by an exact relation which translates in our units into
$$\begin{eqnarray}\langle \unicode[STIX]{x1D700}\rangle =\frac{1}{\sqrt{Gr}}\langle \unicode[STIX]{x1D6FA}\rangle .\end{eqnarray}$$
Figure 10 shows that this relation is satisfied to a good approximation even in our closed cylindrical cell when we compare the vertical mean profiles
$\unicode[STIX]{x1D716}(z)$
and
$\unicode[STIX]{x1D6FA}(z)/\sqrt{Gr}$
. The dashed lines for the energy dissipation rate and the corresponding solid lines for the rescaled local enstrophy are found to collapse almost perfectly.
5 Local statistical analysis
Isosurfaces of the local enstrophy are plotted in figure 11 for
$Pr=0.021$
for two different Rayleigh numbers, along with plots of the temperature field at two horizontal cross-sections in the thermal boundary layers. The plot, which excludes the field at the side walls, displays the fine structure of the characteristic vortex tubes, which are well known from homogeneous isotropic turbulence (Kerr Reference Kerr1985). Enhanced small-scale intermittency was found to be in line with ever finer vortex tubes as discussed by Ishihara et al. (Reference Ishihara, Gotoh and Kaneda2009). This is exactly what we detect in the centre of the convection cell when comparing the two isosurface plots. In the following we want to investigate in detail if this enhanced intermittency in the bulk is in line with enhanced intermittency in the boundary layers, in particular in the velocity boundary layer.
Isosurfaces of constant enstrophy for
$Pr=0.021$
. (a)
$Ra=10^{7}$
and
$1\times 10^{4}<\unicode[STIX]{x1D6FA}<1\times 10^{5}$
. (b)
$Ra=10^{8}$
and
$3.5\times 10^{4}<\unicode[STIX]{x1D6FA}<5\times 10^{4}$
. Also shown are the temperature field at cuts through the bottom and top boundary layers. See supplementary movies available at http://dx.doi.org/10.1017/jfm.2016.457.
5.1 Skin-friction scaling
We start with the definition of the skin-friction coefficient (see e.g. Schlichting & Gersten Reference Schlichting and Gersten2000) which is defined as
$$\begin{eqnarray}c_{\unicode[STIX]{x1D70F}}=\frac{\langle \hat{\unicode[STIX]{x1D70F}}_{w}\rangle }{\hat{\unicode[STIX]{x1D70C}}_{0}\hat{U} ^{2}/2},\end{eqnarray}$$
where
$\langle \cdot \rangle$
is an ensemble average. The coefficient relates the mean wall shear stress to the dynamic pressure. The hat indicates a quantity with dimension. Here,
$\hat{\unicode[STIX]{x1D70C}}_{0}$
is the constant mass density of the fluid and
$\hat{U}$
a characteristic velocity which will be specified further below. The following definition of the wall shear stress field at the plate is used:
$$\begin{eqnarray}\hat{\unicode[STIX]{x1D70F}}_{w}(\hat{x},{\hat{y}},\hat{z}=0,\hat{t})=\hat{\unicode[STIX]{x1D70C}}_{0}\hat{\unicode[STIX]{x1D708}}\sqrt{\left|\frac{\unicode[STIX]{x2202}\hat{u} _{x}}{\unicode[STIX]{x2202}\hat{z}}\right|_{\hat{z}=0}^{2}+\left|\frac{\unicode[STIX]{x2202}\hat{u} _{y}}{\unicode[STIX]{x2202}\hat{z}}\right|_{\hat{z}=0}^{2}}.\end{eqnarray}$$
Following Scheel & Schumacher (Reference Scheel and Schumacher2014) this can be related to a local thickness of the velocity boundary layer by
$$\begin{eqnarray}\hat{\unicode[STIX]{x1D706}}_{v}^{o}(\hat{x},{\hat{y}},\hat{z}=0,\hat{t})=\frac{\hat{\unicode[STIX]{x1D70C}}_{0}\hat{\unicode[STIX]{x1D708}}\hat{u} _{rms,V}}{\hat{\unicode[STIX]{x1D70F}}_{w}(\hat{x},{\hat{y}},\hat{z}=0,\hat{t})},\end{eqnarray}$$
excluding zero-stress events. If we assume that the
$\hat{U}$
in (5.1) is also given by
$\hat{u} _{rms,V}$
(see (2.7)) and express all variables in units of
$H,U_{f}$
and
$\unicode[STIX]{x0394}T$
(see § 2), we obtain the following definition for
$c_{\unicode[STIX]{x1D70F}}$
:
$$\begin{eqnarray}c_{\unicode[STIX]{x1D70F}}=\frac{2}{Re\langle \unicode[STIX]{x1D706}_{v}^{o}\rangle },\end{eqnarray}$$
where
$Re$
is defined in (2.7) and the mean outer velocity boundary layer thickness is found by
$$\begin{eqnarray}\langle \unicode[STIX]{x1D706}_{v}^{o}\rangle =\int _{0}^{\infty }\,\unicode[STIX]{x1D706}_{v}^{o}\,p(\unicode[STIX]{x1D706}_{v}^{o})\,\text{d}\unicode[STIX]{x1D706}_{v}^{o},\end{eqnarray}$$
with
$p(\unicode[STIX]{x1D706}_{v}^{o})$
being the probability density function (PDF) of the outer velocity boundary layer thickness.
Rather than an ensemble average, a plane–time average is performed. Averages over the entire plate area
$r\leqslant 0.5$
will be denoted by
$\langle \cdot \rangle _{A,t}$
or over the interior plate area
$r_{c}<0.3$
by
$\langle \cdot \rangle _{b,t}$
. Thus the two friction coefficients are given as
$$\begin{eqnarray}c_{\unicode[STIX]{x1D70F}}(A)=\frac{2}{Re\langle \unicode[STIX]{x1D706}_{v}^{o}\rangle _{A,t}}\quad \text{and}\quad c_{\unicode[STIX]{x1D70F}}(b)=\frac{2}{Re\langle \unicode[STIX]{x1D706}_{v}^{o}\rangle _{b,t}}.\end{eqnarray}$$
For the interior average, we selected a range of cutoff radii from
$0.2<r_{c}<0.4$
and found that the results for scaling exponents and prefactors do not change significantly. However, there is a systematic trend towards smaller boundary layer thicknesses (and hence larger values of
$c_{\unicode[STIX]{x1D70F}}(b)$
) as
$r_{c}$
decreases. We chose
$r_{c}=0.3$
as our cutoff radii for all plots. The uncertainty in the values in table 4 for interior averages is based on how the results vary as
$r_{c}$
varies. In other words, the scaling exponent was found for
$0.2<r_{c}<0.4$
and the difference between the largest and smallest scaling exponent is given as the uncertainty in scaling exponent. Likewise for the prefactor.
We have plotted
$c_{\unicode[STIX]{x1D70F}}$
as a function of Rayleigh number and Reynolds number for
$Pr=0.005,0.021,0.7$
, and also compare the whole area and the interior averages in figure 12. The resulting power-law fits are given in table 4. In all cases the skin friction magnitude and scaling exponent are reduced when the side-wall regions are included. For the scaling with respect to the Reynolds number (2.7) is used.
Comparing friction coefficient
$c_{\unicode[STIX]{x1D70F}}(A)$
as defined in (5.6) versus Rayleigh number in (a) and versus Reynolds number in (b), for
$Pr=0.005$
(filled circles with black outline),
$Pr=0.021$
(filled circles) and
$Pr=0.7$
(open circles). Also plotted is
$c_{\unicode[STIX]{x1D70F}}(b)$
as defined in (5.6) versus Rayleigh number (or Reynolds number) for
$Pr=0.005$
(filled squares with black outline),
$Pr=0.021$
(filled squares) and
$Pr=0.7$
(open squares). The lines are fits to the data and are given in table 4.
Fits to the data in figure 12 for friction coefficient
$c_{\unicode[STIX]{x1D70F}}$
versus Rayleigh number
$Ra$
averaged over the entire plate area (
$A$
) or over only the interior (
$b$
),
$r<0.3$
. Also given are fits to the data in figures 13 and 14 for local boundary layer thickness versus
$Ra$
averaged over time and the entire plate area (
$\langle \unicode[STIX]{x1D706}_{T}^{o}\rangle _{A}$
, and
$\langle \unicode[STIX]{x1D706}_{v}^{o}\rangle _{A}$
) or over only the interior plate area (
$\langle \unicode[STIX]{x1D706}_{T}^{o}\rangle _{b}$
and
$\langle \unicode[STIX]{x1D706}_{v}^{o}\rangle _{b}$
). Also fit are the boundary layer thicknesses
$\unicode[STIX]{x1D6FF}_{T}$
obtained from Nusselt number and
$\unicode[STIX]{x1D6FF}_{v}$
obtained from the Reynolds number. The uncertainty for all area-averaged quantities and
$\unicode[STIX]{x1D6FF}_{T}$
,
$\unicode[STIX]{x1D6FF}_{v}$
is determined by the scatter in the data and is the result of performing a least squares fit. The uncertainty in the interior averages is discussed just after (5.6).
Which scaling exponents can be expected for a velocity boundary layer in turbulent convection in the present Rayleigh number range? One limiting case for the skin-friction coefficient could be that of a laminar, stationary, two-dimensional, zero-pressure-gradient flat plate boundary layer (Blasius Reference Blasius1908). There
$c_{f}\sim Re_{L}^{-1/2}$
with
$Re_{L}$
being the Reynolds number which is composed of the constant outer inflow velocity
$\hat{U} _{\infty }$
and the downstream length
$\hat{L}$
of the flat plate. The other limiting case for our present flow could be that of a fully developed turbulent flat plate boundary layer which yields a scaling law of the skin-friction coefficient which can be approximated by
$c_{f}\sim Re_{L}^{-1/5}$
for
$5\times 10^{5}\lesssim Re_{L}\lesssim 10^{7}$
(Coles & Hirst Reference Coles and Hirst1969).
Clearly, these limiting cases neglect buoyancy effects which are present in the RBC flow. Furthermore, the flat plate Reynolds number
$Re_{L}$
cannot be directly associated with the present
$Re$
. We see that for the
$Pr=0.7$
case, the magnitude of the scaling exponent for
$c_{f}$
taken for
$A$
is slightly smaller than the laminar limit. When the side-wall regions at the plate are excluded and the average is taken over
$b$
, the scaling exponent is found to be very close to the laminar limit. For the lower Prandtl numbers, the scaling with Reynolds number when the side-wall regions are included agrees, within the error bars, with the exponent of
$-$
0.2 for a turbulent boundary layer although the Reynolds numbers are, by more than an order of magnitude, smaller than
$Re_{L}$
. When the side walls are excluded, the scaling with Reynolds number is found to take values around
$-$
0.25. Not only does the bulk of the convection flow display more vigorous turbulence as
$Pr$
decreases, but the same seems to hold for the velocity boundary layer as the global scaling of the skin-friction coefficient suggests.
5.2 Boundary layer thickness scaling
We can define boundary layer thicknesses in terms of the global heat and momentum transfer (see Grossmann & Lohse Reference Grossmann and Lohse2000) by
$$\begin{eqnarray}\hat{\unicode[STIX]{x1D6FF}}_{T}=\frac{{\hat{H}}}{2Nu},\end{eqnarray}$$
and
$$\begin{eqnarray}\hat{\unicode[STIX]{x1D6FF}}_{v}=\frac{a{\hat{H}}}{\sqrt{Re}}.\end{eqnarray}$$
Note that
$a$
is a free parameter (Prandtl Reference Prandtl1905; Blasius Reference Blasius1908) which is set to
$a=1/4$
in the present case following Grossmann & Lohse (Reference Grossmann and Lohse2001). We define the (dimensionless) mean thermal boundary layer thicknesses by Scheel & Schumacher (Reference Scheel and Schumacher2014) to be
$$\begin{eqnarray}\unicode[STIX]{x1D706}_{T}^{o}=\frac{1}{2}\left|\frac{\unicode[STIX]{x2202}T}{\unicode[STIX]{x2202}z}\right|_{z=0}^{-1},\quad \langle \unicode[STIX]{x1D706}_{T}^{o}\rangle =\int _{0}^{\infty }\,\unicode[STIX]{x1D706}_{T}^{o}\,p(\unicode[STIX]{x1D706}_{T}^{o})\,\text{d}\unicode[STIX]{x1D706}_{T}^{o}.\end{eqnarray}$$
Figure 13 summarizes our findings for the series of low-
$Pr$
runs and compares them with the data obtained by Scheel & Schumacher (Reference Scheel and Schumacher2014). The scaling exponents for convection in air agree to a good approximation with those for liquid-metal convection, but with slightly decreasing exponent magnitude as
$Pr$
decreases, which is consistent with the scaling for
$Nu(Ra)$
. The value for
$\langle \unicode[STIX]{x1D706}_{T}\rangle$
is always larger for the low-
$Pr$
case. The boundary layer thickness is also larger when the side-wall regions are included. It is observed that this difference is consistently larger for the low-
$Pr$
case, even with increasing
$Ra$
. The plumes that detach from the plates are coarser due to enhanced thermal diffusion and cause larger variations in the local boundary layer thickness.
(a) Scaling of mean thermal boundary layer thickness with Rayleigh number for different analysis methods and different Prandtl number. The filled data points with black outlines are for
$Pr=0.005$
, the filled data points are for
$Pr=0.021$
and the open symbols are for
$Pr=0.7$
. The squares for
$\langle \unicode[STIX]{x1D706}_{T}^{o}\rangle _{A}$
are taken over the whole plate area, the diamonds are for
$\langle \unicode[STIX]{x1D706}_{T}^{o}\rangle _{b}$
taken over the interior area with
$r<0.3$
. Circles correspond to the theoretical values for
$\unicode[STIX]{x1D6FF}_{T}=1/(2Nu)$
. All data are also averaged with respect to time. The fits are given in table 4. (b) The same data sets replotted where all the mean thermal boundary layer thickness data (represented by ‘
$\unicode[STIX]{x1D706}_{T}$
’) are compensated by
$Ra^{-0.29}$
.
Figure 14 summarizes our findings for
$\langle \unicode[STIX]{x1D706}_{v}^{o}\rangle$
, which is defined in (5.5). We find that
$\langle \unicode[STIX]{x1D706}_{v}^{o}\rangle$
is significantly smaller for the low-
$Pr$
case. We also find that the scaling exponent for
$\langle \unicode[STIX]{x1D706}_{v}^{o}\rangle$
for the low-
$Pr$
cases becomes larger in magnitude as
$Pr$
decreases. Since
$\langle \unicode[STIX]{x1D706}_{v}^{o}\rangle$
is determined from the inverse of the velocity gradients, this suggests that the magnitudes of these gradients become larger with
$Ra$
at a faster rate for the lower
$Pr$
case. However, the scaling exponent for
$\unicode[STIX]{x1D6FF}_{v}$
as measured from the Reynolds number for
$Pr=0.021$
case is slightly smaller in magnitude than that for convection in air (
$Pr=0.7$
) or liquid sodium (
$Pr=0.005$
), consistent with our results from figure 5.
Scaling of mean velocity boundary layer thickness with Rayleigh number for different analysis methods and different
$Pr$
. The filled data points with black outlines are for
$Pr=0.005$
, the filled data points are for
$Pr=0.021$
and the open symbols are for
$Pr=0.7$
. The squares are for
$\langle \unicode[STIX]{x1D706}_{v}^{o}\rangle _{A}$
taken over the whole cell and the diamonds are for
$\langle \unicode[STIX]{x1D706}_{v}^{o}\rangle _{b}$
taken over the interior plate area. The circles are the theoretical values for
$\unicode[STIX]{x1D6FF}_{v}=a/\sqrt{Re}$
with the Reynolds numbers taken from table 1 and
$a=1/4$
. All data are averaged with respect to time. The fits are again listed in table 4.
5.3 Local shear Reynolds number
We define the local shear Reynolds number as
$$\begin{eqnarray}Re_{sh}(x,y)=\sqrt{\frac{Ra}{Pr}}\unicode[STIX]{x1D706}_{v}^{o}(x,y)u_{rms}(z=\unicode[STIX]{x1D706}_{v}^{o}(x,y)).\end{eqnarray}$$
This is a good measure of the stability of the flow in the boundary layer. It is similar to (41.1) in Landau & Lifshitz (Reference Landau and Lifshitz1987) but we pick
$u_{rms}$
at
$\unicode[STIX]{x1D706}_{v}^{o}$
because the boundary layer thickness varies locally (see Scheel & Schumacher Reference Scheel and Schumacher2014), hence it is more straightforward than selecting
$U$
, the fluid velocity outside the space- and time-averaged boundary layer.
According to Tollmien (Reference Tolmien1929) and Landau & Lifshitz (Reference Landau and Lifshitz1987), when the value of
$Re_{sh}\simeq 420$
is exceeded, then the boundary layer is turbulent. The boundary layer scale that was used for this estimate is the displacement thickness scale, and this applies to two-dimensional turbulent shear flow. We could compute the area- and time-averaged
$Re_{sh}$
to assess when the boundary layer becomes turbulent, which was done in Wagner et al. (Reference Wagner, Shishkina and Wagner2012) and Scheel & Schumacher (Reference Scheel and Schumacher2014). But, since the locally calculated
$Re_{sh}$
(5.10) varies strongly and reaches values greater than 420 for local regions of the plate and for large enough Rayleigh number, we can also assess this spatial intermittency by finding the fractional area which is defined as
$$\begin{eqnarray}\mathscr{T}=\frac{A(Re_{sh}>Re_{sh}^{c})}{A},\end{eqnarray}$$
where the numerator of (5.11) describes the area of the plate where
$Re_{sh}(x,y,t)>Re_{sh}^{c}$
and
$Re_{sh}^{c}$
is a critical shear threshold Reynolds number, i.e. 420 in the case of Landau & Lifshitz (Reference Landau and Lifshitz1987). We can compute
$\mathscr{T}$
very precisely with the code and our area average uses the local volume element for our collocation data points. Also note that the total area
$A$
is given by
$\unicode[STIX]{x03C0}/4$
.
Threshold plots demonstrating the area of the bottom plate whose shear Reynolds number as calculated in (5.10) is greater than 250 (red). The grey area denotes shear Reynolds numbers less than 250. (a,c,e) is for
$Pr=0.7$
and (a)
$Ra=1\times 10^{8}$
(c)
$Ra=1\times 10^{9}$
, and (e)
$Ra=1\times 10^{10}$
and (b,d,f) is for
$Pr=0.021$
and (b)
$Ra=1\times 10^{6}$
, (d)
$Ra=1\times 10^{7}$
. (f)
$Ra=1\times 10^{8}$
. All panels display snapshots.
Since there is no theoretical value for the transition to turbulence in the boundary layer for Rayleigh–Bénard convection, we take three values for
$Re_{sh}^{c}:250,350,450$
. This brackets a range of possible turbulent onset values (including the value of 420 from Landau & Lifshitz (Reference Landau and Lifshitz1987)) and provides calculation of the uncertainty. We plot snapshots of
$\mathscr{T}$
for various Rayleigh and Prandtl numbers in figure 15 to highlight the intermittency. We chose
$Re_{sh}>250$
as our cutoff in the figures as this provides a larger region for
$\mathscr{T}$
. We see in figure 15 that the structures where
$Re_{sh}>250$
are much finer for the higher
$Pr$
case and increase more in number as
$Ra$
increases.
In figure 16 we look at the time evolution of the fractional area and note that it varies by approximately 25 %, as is expected for the intermittency in the boundary layer. We also see variation between the top and bottom plate, with some anticorrelation which we attribute to the large-scale circulation. If, for example, a plume is ejected from the bottom plate, causing a large value of
$\mathscr{T}$
on the bottom plate, it takes some time for this plume to be swept up to the top plate by the large-scale circulation. Eventually this plume hits the top plate, which could cause another plume to be ejected by the top plate, causing
$\mathscr{T}$
to now be larger at the top plate. However, the anticorrelation is not perfect and the situation is more complex since plume ejection occurs intermittently.
Finally we plot the time-averaged
$\mathscr{T}$
in figure 17 and see that the turbulent fraction increases as a function of Rayleigh number for all Prandtl numbers, with the turbulent fraction always being larger for the lower Prandtl number.
In figure 17 we also can determine the critical Rayleigh number
$Ra_{Ac}$
where
$\langle \mathscr{T}\rangle _{t}=1$
for our Prandtl numbers and
$Re_{sh}^{c}$
values, which would indicate a transition to a fully turbulent boundary layer. We find, using
$Re_{sh}^{c}=350$
as our midrange value, and
$Re_{sh}^{c}=250$
, 450 as upper and lower bounds, that
$Ra_{Ac}=(3\pm 2)\times 10^{13}$
for
$Pr=0.7$
,
$Ra_{Ac}=(2\pm 10)\times 10^{13}$
for
$Pr=0.021$
. We did not fit the data or estimate an
$R_{Ac}$
for
$Pr=0.005$
since the data cover an interval which is too small and the Rayleigh numbers are too low for a reasonable extrapolation.
The uncertainty is high in these predicted critical Rayleigh numbers, since we are using data at Rayleigh numbers that are three orders of magnitude smaller than the predicted transition values. However, the critical Rayleigh number for
$Pr=0.7$
is close to the value of
$Ra_{1}^{\ast }\simeq 2\times 10^{13}$
found by Ahlers et al. (Reference Ahlers, He, Funfschilling and Bodenschatz2012), He et al. (Reference He, Funfschilling, Bodenschatz and Ahlers2012) for similar parameters. In these papers the authors noted that
$Ra_{1}^{\ast }$
indicated a transition from the classical laminar scaling of Nusselt number with Rayleigh number, but the fully turbulent scaling of Nusselt number with Rayleigh number was not seen until
$Ra_{2}^{\ast }\geqslant 5\times 10^{14}$
.
Even though
$\langle \mathscr{T}\rangle _{t}$
increases as
$Pr$
decreases, we still find
$Ra_{Ac}$
occurs at about the same Rayleigh number for
$Pr=0.021$
, when the uncertainty is taken into account. This is because the growth of
$\langle \mathscr{T}\rangle _{t}$
with Rayleigh number is so much smaller for the lower
$Pr$
. The theoretically predicted trend by Grossmann & Lohse (Reference Grossmann and Lohse2000) is for
$Ra_{Ac}$
to decrease as
$Pr$
decreases, and the values as estimated from figure 2 of Grossmann & Lohse (Reference Grossmann and Lohse2000) are around
$10^{13},10^{14},10^{16}$
for
$Pr=0.005,0.021,0.7$
, respectively. This disagreement may be due to the higher uncertainty in these data points, or it may be indicative of something more fundamental, i.e. that the transition from a laminar to a turbulent boundary layer may be more universal than previously predicted. However, more data at higher Rayleigh numbers would provide a more precise value for
$Ra_{Ac}$
.
Time dependence of fractional area
$\mathscr{T}$
(defined in (5.11) greater than 250 (black line
$=$
top plate, magenta (light grey) line
$=$
bottom plate). The fractional area is given for each statistically independent snapshot. (a)
$Pr=0.7$
and
$Ra=1\times 10^{9}$
. (b)
$Pr=0.021$
and
$Ra=1\times 10^{7}$
.
Time-averaged turbulent fraction
$\langle \mathscr{T}\rangle _{t}$
of top and bottom plate area whose shear Reynolds number defined in (5.10) is greater than 350 (black circles,
$Pr=0.7$
, magenta stars,
$Pr=0.021$
, yellow filled diamonds with black outline,
$Pr=0.005$
). The error bars correspond to the cases where
$Re_{sh}^{c}=250$
(upper) and
$Re_{sh}^{c}=450$
(lower). Note this fraction has been averaged over many statistically independent time slices. The thick black line corresponds to
$\mathscr{T}=1$
. The fit to the last four data points for
$Pr=0.7$
is
$\langle \mathscr{T}\rangle _{t}=(7\pm 1\times 10^{-7})Ra^{0.5\pm 0.1}$
(dashed line). The fit to all data points for
$Pr=0.021$
, is
$\langle \mathscr{T}\rangle _{t}=(1.6\pm 1\times 10^{-4})Ra^{0.28\pm 0.05}$
(solid line).
Kurtosis (5.12) versus Rayleigh number for
$Pr=0.7$
(circles),
$Pr=0.021$
(stars) and
$Pr=0.005$
(diamonds).
We can also check for intermittency by evaluating the kurtosis, as defined in (5.12). For a Gaussian distribution the kurtosis should be equal to three. Anything larger than this indicates wider tails to the distribution, and hence intermittency. The kurtosis is plotted in figure 18 for all Rayleigh and Prandtl numbers. We see that the kurtosis is much larger than three for all data. For
$Pr=0.7$
, the kurtosis stays fairly level at 500 until
$Ra=1\times 10^{8}$
and then it increases with
$Ra$
. However, for the lower
$Pr$
data the kurtosis increases with
$Ra$
for all data and the kurtosis values are also very close to one another:
$$\begin{eqnarray}\text{Kurtosis}={\displaystyle \frac{\langle (Re_{sh}-\langle Re_{sh}\rangle )^{4}\rangle }{\langle (Re_{sh}-\langle Re_{sh}\rangle )^{2}\rangle ^{2}}}.\end{eqnarray}$$
6 Summary and discussion
A comparison of the global and local statistics has been presented for low
$Pr$
numbers (
$Pr=0.005,0.021$
) and those of air (
$Pr=0.7$
) with very finely resolved DNS. The low-
$Pr$
cases are in a parameter regime that has not been studied numerically at such fine resolution. It is found that the temperature field is more diffuse while the velocity field is more vigorous as the Prandtl number is lowered. It is also found that as the Prandtl number decreases, there is a slight overall decrease in the value of the scaling exponents of diagnostic quantities, such as Nusselt number and thermal boundary layer thicknesses, with Rayleigh number. We attribute this to the more diffuse temperature field and resulting reduced heat transport as
$Pr$
decreases to these very low values.
The global heat transfer scaling exponents are found to have values of
$0.29\pm 0.01$
,
$0.26\pm 0.01$
and
$0.265\pm 0.01$
for
$Pr=0.7,0.021$
and
$0.005$
respectively. When compared with experimental results of Horanyi et al. (Reference Horanyi, Krebs and Müller1999), the present exponents are found to agree quite well for
$Pr=0.005$
. The exponent for
$Pr=0.021$
agrees with the experiments by Rossby (Reference Rossby1969), Cioni et al. (Reference Cioni, Ciliberto and Sommeria1996), Takeshita et al. (Reference Takeshita, Segawa, Glazier and Sano1996), King & Aurnou (Reference King and Aurnou2013). However our exponent is on the high end of the range, and disagrees with the theory which predicts a scaling exponent of one-fifth for this parameter regime in Grossmann & Lohse (Reference Grossmann and Lohse2000). The global momentum scaling exponents are
$0.49\pm 0.01,0.44\pm 0.01,0.49\pm 0.01$
for
$Pr=0.7,0.021,0.005$
. For the case of
$Pr=0.021$
, they agree well with the experimental results by Takeshita et al. (Reference Takeshita, Segawa, Glazier and Sano1996). Finally it is found that when the Nusselt number is plotted versus Prandtl number, the agreement with the Grossmann–Lohse theory is quite good, except for a slight disagreement in the lower Prandtl numbers (
$0.021,0.005$
). It is stressed once more that the number of data points of our DNS is small and that the predictions of the scaling exponents for global heat and mass transfer should be considered with caution.
The root-mean-square profiles of velocity show a significant enhancement of velocity with decreasing
$Pr$
and a shift of the maxima towards the centre of the container. The root-mean-square profiles of temperature show similar values for the location of the maxima for
$Pr=0.7,0.021,0.005$
, but a reduction in overall root-mean-square temperature for
$Pr=0.005$
. Finally we see that the magnitude of the vertical mean profiles of
$\unicode[STIX]{x1D716}$
and
$\unicode[STIX]{x1D716}_{T}$
systematically increases as
$Pr$
decreases, as is expected. The kinetic dissipation rate increases due to the enhanced velocity fluctuations and the thermal dissipation rate increases due to the higher diffusivity of the temperature field. The agreement between profiles of
$\unicode[STIX]{x1D716}(z)$
and
$\unicode[STIX]{x1D6FA}(z)/\sqrt{Gr}$
is also very good, similar to that found for homogeneous isotropic turbulence.
When we compute the thermal boundary layer thicknesses we find that they are also consistently larger as
$Pr$
decreases, as expected. Also, the magnitude of the scaling exponents decrease slightly as
$Pr$
decreases, consistent with the
$Nu$
scaling. In contrast, the velocity boundary layer thicknesses are consistently smaller as
$Pr$
decreases and the scaling exponent magnitudes increase as
$Pr$
decreases. The effect of the side walls is found to remain significant for the lower
$Pr$
cases (
$0.021$
and
$0.005$
) as Rayleigh number increases but it becomes less significant for
$Pr=0.7$
.
Turning now to the local statistics, it is found that the skin-friction coefficient
$c_{\unicode[STIX]{x1D70F}}$
is significantly reduced as
$Pr$
decreases for all Rayleigh numbers, averaging either over the central region or the whole plate. It is also found that the scaling exponent with Reynolds number when the side-wall regions are included increases from
$-0.54\pm 0.01$
for
$Pr=0.7$
(consistent with a laminar boundary layer) to
$-0.23\pm 0.04$
for the lowest
$Pr$
which bears a closer resemblance to an intermittently turbulent boundary layer. Inspired by this, we compute a local shear Reynolds number and resulting turbulent fraction
$\langle \mathscr{T}\rangle _{t}$
. We find that this turbulent fraction increases with Rayleigh number and is larger for the lower Prandtl numbers. By extrapolating to
$\langle \mathscr{T}\rangle _{t}=1$
we can predict the Rayleigh number
$Ra_{Ac}$
at which the boundary layer would become fully turbulent. We find
$Ra_{Ac}=(3\pm 2)\times 10^{13}$
for
$Pr=0.7$
and
$Ra_{Ac}=(2\pm 10)\times 10^{13}$
for
$Pr=0.021$
. It is stressed once more that our prediction is based on an extrapolation which might be wrong once the gap to the high Rayleigh numbers is filled with further DNS data points.
The present analysis has confirmed the results of Schumacher et al. (Reference Schumacher, Götzfried and Scheel2015), namely that the fluid turbulence in low-
$Pr$
convection is much more strongly driven and reaches a more vigorous turbulence level than comparable convection flows in air or water. The reason for this observation is that the bottleneck in the RBC system – the thermal boundary layer – is thicker due to enhanced diffusion. Coarser thermal plumes can inject kinetic energy at a larger scale as shown in Schumacher et al. (Reference Schumacher, Götzfried and Scheel2015) and extend the Kolmogorov-like cascade range in the bulk. This is in line with an enhanced energy flux in the inertial range and a smaller Kolmogorov length. Coarser plumes intensify also the large-scale circulation in the cell and thus amplify the transient and intermittent behaviour in the velocity boundary layer which was shown here, e.g. by the skin friction or the local shear Reynolds number. Our current studies indicate that low-Prandtl-number convection at the same Grashof number may be more susceptible to a transition to turbulence in the tiny boundary layer when compared with convection in air.
Given the very low-
$Pr$
data records, a next step would be to refine the analysis in the transient boundary layers and to develop alternative measures that quantify the increasing number of turbulent patches in the velocity boundary layer. These studies are currently under way and will be reported elsewhere.
Acknowledgements
Helpful discussions with J. Aurnou, R. du Puits, P. Frick and B. Eckhardt are acknowledged. The work is supported by Research Unit 1182 and the Research Training Group 1567 of the Deutsche Forschungsgemeinschaft. We acknowledge supercomputing time at the Blue Gene/Q JUQUEEN at the Jülich Supercomputing Centre which was provided by grant HIL09 of the John von Neumann Institute for Computing. Furthermore, we acknowledge an award of computer time provided by the INCITE program. This research used resources of the ALCF at ANL, which is supported by the DOE under contract DE-AC02-06CH11357.
Supplementary movies
Supplementary movies are available at http://dx.doi.org/10.1017/jfm.2016.457.
