2.1. Experimental set-up: FOCUS

The experimental set-up, named FOCUS and depicted in figure 1, consists of an upright acrylic cylindrical enclosure with a vertical main axis and an aspect ratio $\Gamma =R/H=4$, with $R$ the tank radius and $H$ its height. Inside this enclosure, a gas (air or air/${\rm SF}_6$ mixture) is heated from below and cooled along the side.

Figure 1. $(a)$ Photo of the experimental set-up with the laser sheet inclined for better visualisation (it is horizontal during acquisition). The dotted yellow rectangular region shows the PIV camera field of view. $(b)$ Sketch of a vertical cross-section of the experimental set-up with magnified views of the bottom plate configuration below the working tank and below the side boundary, respectively.

The choice of a gas as a simulant was motivated by the Prandtl number characteristic of the metal layer in the nuclear application, which is less than 1 ($Pr = \nu /\kappa \simeq 0.1$, with $\nu$ the kinematic viscosity and $\kappa$ the thermal diffusivity). When $Pr \lt 1$, the viscous boundary layer is nested within the thermal boundary layer. This structure plays a major role in heat transfers (Rein et al. Reference Rein, Carénini, Fichot, Favier and Le Bars2023a). Using gases with $Pr\simeq 0.7 \lt 1$ ensures maintaining a similar configuration, thus preserving similar heat transfer mechanisms, without the cost and many difficulties of using liquid metals (including ,for instance, their opacity or the necessity to work at higher than ambient temperature). With ${\rm SF}_6$ being denser than air, it allows more turbulent regimes to be reached by significantly increasing the Rayleigh number: indeed, the thermal conductivity and the dynamic viscosity of air and ${\rm SF}_6$ are mostly similar, so the Rayleigh number (see the exact definition (3.5)) scales as the gas density squared. Note that such gases have already been successfully used to study high-Rayleigh-number thermal convection in the past; see, e.g. Ashkenazi & Steinberg (Reference Ashkenazi and Steinberg1999); Ahlers et al. (Reference Ahlers, Funfschilling and Bodenschatz2009a); and more recently Cierpka et al. (Reference Cierpka, Kästner, Resagk and Schumacher2019).

In our set-up (see figure 1), the gas layer is heated from below through eight silicone heating elements ($450\,\text {mm}\times 250\,\text {mm}$ each) with a combined maximum power of 7.2 kW, resting on a cylindrical aluminium plate with a radius of $R= {44}\,\text {cm}$, for better heat flux homogenisation. These heating elements are positioned on a 5 cm thick calcium silicate insulating plate to limit heat losses and to concentrate heat into the aluminium plate. Cooling on the side is provided by a cylindrical polymethyl methacrylate (PMMA) torus with a rectangular section of 11 cm in height and 40 mm in width, through which water flows at an average rate of 15 L min^–1. The water temperature is regulated by a thermostatic bath with a power of 1.2 kW. Three water inlets and outlets located at the bottom and top of the torus, arranged in a staggered pattern and spaced at ${60}^\circ$ intervals, promote water mixing through natural and forced convection. The thickness of the inner wall of the torus is only 4 mm (minimum thickness achievable through manufacturing processes), so as to ensure a temperature as constant as possible at the wall. Note that on the edges of the bottom plate, a thin 3 mm cork plate was added to prevent direct contact between the heating elements and the bottom of the cooling torus, while holes were drilled in the insulating plate to dissipate heat. Note also that insulated walls can introduce side wall effects, especially when gas is used (Roche et al. Reference Roche, Castaing, Chabaud, Hébral and Sommeria2001). In our set-up however, we expect that the temperature-regulated water inside the tank, controlled by the powerful thermal bath, largely compensates for any residual heat flux.

A transparent polycarbonate lid with a thickness of 35 mm is used to maintain the gas in the tank and to constrain the outgoing heat flux, while allowing flow visualisation from above. This outgoing heat flux through the lid is a parameter of our study and accounts for approximately 50 % of the input heat flux across all experiments (see Appendix A.1). The choice of using a polycarbonate cover was guided by the temperature range that can be reached with the experimental set-up, in order to prevent any changes in lid opacity and ensure optimal flow visualisation.

Finally, to minimise the ${\rm SF}_6$ mixing with the external environment, the entire set-up is enclosed in a large square PMMA tank measuring ($1220\,\text {mm}\times 1220\,\text {mm}\times 470\,\text {mm}$).

The operational temperature range of the device has been limited to less than $1000\,^\circ{\rm C}$, reaching a maximal bulk temperature of $64\,^\circ{\rm C}$. Indeed, since we use gas as a simulant for the incompressible metal layer, we want to avoid excessive temperature differences that would result in significant and undesirable compressibility effects. For an ideal gas, the density variation can be estimated as the ratio of the largest temperature change between the initial and statistically stationary states to the mean temperature in the statistically stationary state. Throughout all the experiments, we measured an average relative density variation of 5 %, peaking at 11.5 % for the highest Rayleigh number obtained.

2.2. Instrumentation

2.2.1. Heat flux sensor

In order to measure heat flux fluctuations through the side wall, we use a thermal flux sensor (‘gSkin’ manufactured by ‘greenTEG’) measuring $18\,\text {mm}\times 18\,\text {mm}\times 0.5\,\text {mm}$ and positioned on the inner wall of the torus. It has a response time of the order of 0.7 s and its sensitivity depends on temperature via the relation

(2.1) \begin{equation} s=S_0 + (\theta _0-\theta _S)S_c, \end{equation}

with $S_0\,=\,{39.92}\,{\unicode{x03BC} }\text {V}\,\text {W}^{-1}\,\text {m}^2$, $\theta _S=22.5\,{}^\circ \text {C}$, $S_c\,=\, 0.0499\,{\unicode{x03BC} }\text {V}\,\text {W}^{-1}\,\text {m}^2\,{}^\circ \text {C}^{-1}.$ Here $\theta _0$ corresponds to the temperature value of the substrate on which the flux sensor is placed, hence, the cooling temperature here. During the measurement campaign, the sensor was mainly placed at mid-height of the torus, using thermal paste to avoid any parasitic thermal resistance and to ensure good adhesion to the wall over the entire surface despite the local curvature.

2.3. Experimental protocol

We now detail the experimental protocol implemented during the measurement campaign. We begin by initiating the acquisition system for temperature and flux. In the case of ${\rm SF}_6$ usage, plastic sheets are placed above the set-up and fixed to the outer basin to contain the gas. The outer basin is first filled completely with ${\rm SF}_6$ from the bottom, until gas is detected on the surface of the sheets using a gas detector. Then, the experimental set-up is filled with ${\rm SF}_6$ while aiming to maximise density. There, we temporarily use a densimeter (WIKIA Tech, Avenisense NorthDome) placed inside the system, with its sensor positioned 2 cm above the aluminium plate, near the edge opposite to the flux sensors.

Water circulation and water cooling via the thermostated bath are initiated. The bath temperature set point is chosen to maintain a difference of approximately $5\,^\circ \text {C}$ with the ambient temperature, aiming to keep the same external heat losses between each experiment (see details in Appendix A.3). Finally, the bottom heating system is activated by selecting the electrical power to be supplied to the heating elements.

We now enter the transient regime. When using ${\rm SF}_6$, it is expected that the density inside the set-up decreases as it is not gas tight, and ${\rm SF}_6$ diffuses into the ambient air. If the density drops below the desired threshold, more ${\rm SF}_6$ is injected into the system until the desired density is reached. The gas in the system is then a mixture of air and ${\rm SF}_6$. The mass fraction of ${\rm SF}_6$, denoted $y$, is determined via a density measurement further detailed in Appendix A.2.

Figure 2. Sketch of the modelled fluid layer in a vertical plane through the cylinder. The top boundary condition is no slip for the gas experiment simulations, representing the lid used in the set-up, while a free-slip condition is applied to model the metal layer in the nuclear application.

During this transient regime, the characteristic time and thermal losses are evaluated (see Appendix A.3). The steady-state regime is achieved when the heating time exceeds 3 characteristic times, i.e. in about approximately $9\,\text{h}$. At this stage, it must be ensured that the temporal variation of temperature is less than $1\,^\circ\text{C}\,\text{h}^{-1}$.

Once the steady state is reached, flux data are collected for 3–6 diffusive times (approximately $1\,\text{h}$). During this period, the necessary elements for PIV measurements are prepared, including the water/UCON oil mixture to generate PIV particles. After the flux data acquisition is completed, the seeding phase begins. A low-power horizontal laser sheet (0.2 W) is set up, an air/${\rm SF}_6$ injection flow rate through the atomizer is fixed, and the set-up is seeded for a predefined duration to control the quantity of injected PIV particles, ensuring experimental reproducibility. The injection of these particles disturbs the flow and the thermal fluxes involved. It is necessary to wait for the particle jet to homogenise and the thermal signals to return to steady state before starting the flow video acquisition. The laser power, gain, sampling frequency and exposure time are then adjusted to identify and track the movements of the PIV particles. A series of video acquisitions are then conducted, each for a duration of approximately one minute (a duration for which we verify a constant sampling frequency).

3.1. Governing equations

Following Rein et al. (Reference Rein, Carénini, Fichot, Favier and Le Bars2023a), we model the experimental set-up considering an incompressible fluid in the Boussinesq approximation, confined in a cylinder of thickness $H$ and radius $R$, with gravity directed downward (see figure 2). Heating from below is applied with a uniform heat flux $\phi _{{in}}$, while cooling from above occurs with a uniform outgoing flux $\phi _{{out}}$. Note that this upper boundary condition is an approximation, since the experimental outgoing flux is dynamically constrained by the convective heat transfers within the system and may change radially. Nonetheless, the thermal budget described in Appendix A.3 confirms that the surface average measurement and the local, centre plate measurement of the outgoing flux have close values, hence validating this assumption with an outgoing flux $\phi _{{out}}$ ranging between $40\,\%\;\text{and}\;60\,\%$ of the incoming one. The study hence focuses on cases where $\phi _{{in}}$ and $\phi _{{out}}$ are unequal, resulting in residual heat flux escaping through the side boundary, whose temperature is fixed at $\theta _0$. No-slip conditions are applied to the bottom and side boundaries, while two configurations are considered for the top boundary in the DNS: a no-slip condition to represent the experimental lid and a free-slip condition to model the liquid metal layer in the nuclear application.

Lengths are rescaled using the height of the cylinder $H$, while time is rescaled using the vertical diffusive dash–dot time scale $H^2/\kappa$. In our set-up, the two main control parameters are the bottom heat flux and the side temperature, which is the only imposed temperature. Hence, the non-dimensional temperature $T$ is defined relative to the imposed side temperature and is rescaled using the imposed bottom flux $\phi _{{in}}$,

(3.1) \begin{equation} T=\frac {k}{\phi _{{in}}H}(\theta -\theta _0) \ , \end{equation}

where $k$ is the thermal conductivity. The non-dimensional conservation equations of momentum, mass and energy are then

(3.2) \begin{equation} \begin{aligned} \frac {1}{Pr}\left (\frac {\partial \boldsymbol {u}}{\partial t}+ \boldsymbol {u\cdot \nabla }\boldsymbol {u}\right ) =-\boldsymbol {\nabla }P +Ra_{\phi } T\boldsymbol {e}_z+\boldsymbol {\nabla }^2\boldsymbol {u} \ , \end{aligned} \end{equation}

(3.3) \begin{equation} \boldsymbol {\nabla \cdot u}=0 \ , \end{equation}

(3.4) \begin{equation} \frac {\partial T}{\partial t}+ \boldsymbol {u}\cdot \boldsymbol {\nabla } T=\nabla ^2 T \ . \end{equation}

Here $\boldsymbol {u}$ and $P$ are the non-dimensional velocity and pressure (including the hydrostatic contribution) fields. The problem is characterised by four non-dimensional parameters: the aspect ratio $\Gamma$, the flux ratio $R_F$, the Rayleigh number $Ra_{\phi }$ based on the heat flux imposed at the bottom $\phi _{{in}}$, the Prandtl number $Pr$ set to $0.1$ to mimic the properties of the metal layer or $0.7$ to replicate those of the gas used in the experiment. They are defined by

(3.5) \begin{equation} \Gamma = \frac {R}{H},\quad R_F=\frac {\phi _{{out}}}{\phi _{{in}}},\quad Ra_{\phi }=\frac {\beta g \phi _{{in}} H^4}{k \nu \kappa },\quad Pr=\frac {\nu }{\kappa }=\left \{\begin{array}{ll} 0.1 & {(\rm application)},\\[4pt] 0.7 & {(\rm experiment)}, \end{array} \right . \ \end{equation}

where $\beta$ is the thermal expansion coefficient, assumed to be constant. The non-dimensional boundary conditions can be written as

(3.6) \begin{equation} \begin{aligned} \boldsymbol {u}(z=0)=\boldsymbol {0}\,\,\,\mbox {and}\,\,\,\frac {\partial T}{\partial z}\bigg|_{z=0}=-1 \ ,\\ \boldsymbol {u}(r=\Gamma )=\boldsymbol {0}\,\,\,\mbox {and}\,\,\, T(r=\Gamma )=0 \ , \\ \left \{\begin{array}{ll} {\rm Application:} & \frac {\partial u}{\partial z}\big|_{z=1}=\frac {\partial v}{\partial z}\big|_{z=1}=w(z=1)=0\\[3pt] {\rm Experiment :} & \boldsymbol {u}(z=1)=\boldsymbol {0} \end{array} \right . \ \mbox {and}\,\,\,\frac {\partial T}{\partial z}\bigg|_{z=1}=-R_F, \end{aligned} \end{equation}

where $\boldsymbol {u}= (u,v,w )$ are the velocity components in cylindrical coordinates along the unit vectors $\boldsymbol {e}_r,\boldsymbol {e}_{\varphi },\boldsymbol {e}_z$.

3.2. Numerical approach

The governing equations (3.2)--(3.4) with corresponding boundary conditions (3.6) are computationally solved using the spectral element code Nek5000 (Fischer Reference Fischer1997; Deville, Fischer & Mund Reference Deville, Fischer and Mund2002). The cylindrical geometry is discretised with up to $\mathcal {E}=9984$ hexahedral elements, refined near all boundaries to accurately resolve viscous and thermal boundary layers, with typically 15 grid points across the viscous boundary layers, reaching a minimum of 10 grid points for the highest Rayleigh number.

Numerical simulations are initiated with a quiescent fluid and a uniform temperature field $T=0$ throughout the domain. Small temperature perturbations of amplitude $10^{-3}$ are introduced, leading to thermal convection growth during a transient phase lasting approximately $5$ vertical diffusive times, with longer durations observed for higher aspect ratios ($\Gamma$). Once the system reaches a statistically stationary state, various spatio-temporal averages are computed. Temporal and volume average operators $\left \lt \cdot\right \gt$ are defined over the entire fluid domain volume $V$ and time $\tau$ as

(3.7) \begin{equation} \left \lt .\right \gt =\frac {1}{\tau V}\int _{t_0}^{t_0+\tau }\int _{V}\cdot\mathrm {d}V\mathrm {d}t ,\ \end{equation}

where $\tau$ typically ranges between $2$ and $0.1$ diffusive times for the lowest and highest Rayleigh numbers, respectively. Additionally, adding specific variables as a subscript means that an average along those specific directions is made. We always consider temporal averages during the statistically steady state so that the time variable is never explicitly written. For instance, $\left \lt .\right \gt _{\varphi }$ indicates an average in time and along the azimuthal direction only.

Although most results are derived from DNS, some extreme cases necessitated filtered simulations, as described in Fischer & Mullen (Reference Fischer and Mullen2001). A viscous dissipation criterion was employed to distinguish between DNS and filtered simulations. For further details, readers are directed to Rein et al. (Reference Rein, Carénini, Fichot, Favier and Le Bars2023a).

4.1. Scalings

In figure 3(a) the non-dimensional mean temperature is plotted against the Rayleigh number. To do so, we averaged in time each thermocouple signal and then averaged over the six thermocouples within the system. Using DNS, we compute the mean temperature following the same set of discrete locations as in the experiments, considering both $R_F=0.5$ and $R_F=0.6$ (which encompasses most of the experimental values for the flux ratio; see table 1), and spanning a range of $Ra_{\phi }$ from $10^5$ to $10^8$.

Figure 3. $(a)$ Non-dimensional temperature averaged in time and over the six thermocouples in the bulk of the experiment, as a function of $Ra_{\phi }$. The dash–dot line shows the $Ra_{\phi }^{-1/5}$ asymptotic scaling from Rein et al. (Reference Rein, Carénini, Fichot, Favier and Le Bars2023a) (the inset shows the compensated plot). $(b)$ Non-dimensional horizontal r.m.s. velocity in the $(x,y)$ plane at height $z= {6}\,\text {cm}$ as a function of $Ra_{\phi }$. The dash–dot line and continuous line show a $1/3$ and $2/5$ best fit power law based on experimental data. For both plots, square ($\square$), diamond ($\Diamond$) and circle ($\circ$) symbols represent respectively data from DNS simulations at $R_F=0.5$ and $R_F=0.6$, and the experiments. Empty circles indicate when ${\rm SF}_6$ gas is used. For all DNS simulations, $\Gamma =4$ and $Pr=0.7$.

Overall, a satisfactory agreement between DNS and experimental measurements is observed. One can notice here the slight influence of the flux ratio. Indeed, the two lowest experimental measurement points in terms of $Ra_{\phi }$ (close to $Ra_{\phi } \sim 10^7$) correspond to flux ratios $R_F\simeq 0.5$ (see Appendix A.3) and align well with the numerical trend observed for $R_F=0.5$. Similarly, the air experimental measurements near $Ra_{\phi } \sim 10^8$, associated with flux ratios close to $0.6$, show excellent agreement with the numerical trend at $R_F=0.6$. The experimental measurements conducted with ${\rm SF}_6$ have a ratio close to $0.6$, and extend the numerical trend by nearly two orders of magnitude in terms of $Ra_{\phi }$.

Also shown is the slope of the asymptotic scaling in $Ra_{\phi }^{-1/5}$ predicted by Rein et al. (Reference Rein, Carénini, Fichot, Favier and Le Bars2023a) and which corresponds to the classical scaling of vertical convection predicted by Batchelor (Reference Batchelor1954); the inset graph shows the results compensated by this scaling with a good agreement above $Ra_{\phi } =10^8$. The scaling, as detailed in Rein et al. (Reference Rein, Carénini, Fichot, Favier and Le Bars2023a), comes from the balance at a statistically stationary state between the top/bottom flux mismatch, and the lateral outgoing diffusive flux through the thermal boundary layer at the vertical boundary. As the Rayleigh number increases, the system becomes more turbulent leading to more efficient bulk mixing limited by the diffusive boundary layer, and hence, a better agreement with the scaling law.

In figure 3(b) we plot the norm of the non-dimensional horizontal velocity value $U_h$ as a function of the Rayleigh number for numerical simulations at $R_F =0.5$ and $0.6$, and the experimental data. This velocity is computed as follows:

(4.1) \begin{equation} U_h=\left [\left \lt u^2\right \gt _{r \varphi }+\left \lt v^2\right \gt _{r \varphi }\right ]^{1/2}. \end{equation}

The averaging procedure along the azimuthal and radial directions involves the same spatial domain for the experiments and DNS (i.e. the PIV field of view; see § 2.2.3). The time averaging is processed during a characteristic time of 1 min for the experiments (0.1 diffusive time for air, 0.02 for ${\rm SF}_6$) and over 0.1 diffusive time for DNS. Each point represents a measurement of the characteristic horizontal fluid velocity localised in the PIV measurement zone. First, it is worth noting the excellent agreement between the experimental measurements and DNS over the range of Rayleigh numbers common to both approaches. In particular, the experimental points between $Ra_{\phi }=10^7$ and $Ra_{\phi }=10^8$, conducted in air, compare favourably with the numerical trend obtained at $Ra_{\phi }=10^8$. Additionally, the experimental measurements involving the use of ${\rm SF}_6$ extend the trend over almost two orders of magnitude in $Ra_{\phi }$. A power law relationship compatible with $Ra_{\phi }^{1/3}$ is observed, with a best fit estimate of $Ra_{\phi }^{0.36}$ using the least squares method based on the experimental points. Note that in the numerical simulations, the root-mean-square (r.m.s.) velocity exhibits minimal dependence on the flux ratio. This observation aligns with our previous findings using $Pr=0.1$ and free-slip boundary conditions applied to the upper surface (not shown). In the following section we use dimensional analysis to gain insight into the physics underlying this scaling law.

Figure 4. $(a)$ Compensated horizontal advective flux following (4.4) as a function of $Ra_{\phi }$. $(b)$ Compensated radial r.m.s. velocity following (4.5) as a function of $Ra_{\phi }$. For both plots, different symbols correspond to different aspect ratios $\Gamma$ and colour variation from red to green involves $R_F$ variation from 0.1 to 0.9. Data are from DNS with $Pr=0.1$ and a free-slip upper boundary, for which our data set is more complete in terms of $\Gamma$ and $R_F$ systematics, but the same behaviour is expected with the experimental configuration. Empty/full symbols indicate respectively filtered/DNS simulations. The shadings represent three times the standard deviation of the horizontal advective flux time series at the statistically stationary state.

4.2. Dimensional analysis

Let us focus on the heat transport by considering the following non-dimensional energy equation derived from (3.4) using cylindrical coordinates:

(4.2) \begin{equation} \frac {\partial T}{\partial t} + \frac {1}{r}\frac {\partial }{\partial r}\left [r \left (uT - \frac {\partial T}{\partial r}\right )\right ] + \frac {1}{r}\frac {\partial }{\partial \varphi }\left [vT - \frac {1}{r}\frac {\partial T}{\partial \varphi }\right ]+ \frac {\partial }{\partial z}\left [wT - \frac {\partial T}{\partial z}\right ]=0. \end{equation}

The horizontal transport of the heat flux at a radius $r$ is expressed by averaging the energy equation (4.2) in time, vertical and azimuthal directions, which leads to

(4.3) \begin{equation} \left \lt uT \right \gt _{\varphi z} - \left \lt \frac {\partial T}{\partial r} \right \gt _{\varphi z} =\frac {(1-R_F)}{2} r. \end{equation}

This equation illustrates that the difference in heat flux between the top and bottom (i.e. the right-hand side of (4.3)) is transported horizontally through a combination of diffusion ($\left \lt \partial _r T \right \gt _{\varphi z}$) and advection ($\left \lt uT \right \gt _{\varphi z}$) processes. If we now assume the system to be fully turbulent, the radial diffusive gradient matters only within the thermal boundary layer. Within the bulk, we expect the balance $\left \lt uT \right \gt _{\varphi z} \sim ({(1-R_F)}/{2}) r$. By averaging along the radial direction, the horizontal advective flux scales as follows:

(4.4) \begin{equation} \langle uT \rangle = \frac {1}{\Gamma }\int ^{\Gamma }_0 \left \lt uT \right \gt _{\varphi z} \mathrm {d}r \sim (1-R_F)\Gamma . \end{equation}

This scaling is confirmed in figure 4(a), where we use DNS data with $\Gamma$ ranging from $4$ to 16, $R_F$ from 0.1 to 0.9 and $Pr=0.1$. The higher the Rayleigh number, the better the agreement with the scaling given by (4.4). In particular, above $Ra_{\phi } \sim 10^5$, $\langle uT \rangle$ becomes independent of the Rayleigh number and all the data collapse on the same straight line.

We now seek to explain the scaling of the horizontal velocity. Following the determination of the free fall velocity scaling in classical Rayleigh–Bénard convection, we assume that the dimensional horizontal velocity scale is independent of viscosity and diffusivity, but depends on its driving buoyancy, i.e. the horizontal turbulent buoyancy flux $\beta g F_h \sim \beta g \phi /\rho c_p$, related to $\langle uT \rangle$ in non-dimensional form.

Scaling analysis gives $\hat {U}^{{rms}}_h\sim (\beta g F_h H)^{1/3}$, where $\hat {\cdot}$ denotes dimensional variable. In non-dimensional form, this becomes ${U}^{{rms}}_h\sim (Ra_{\phi } \Pr \langle uT \rangle )^{1/3}$. Hence, with (4.4),

(4.5) \begin{equation} {U}^{{rms}}_h\sim (Ra_{\phi } \Pr (1-R_F)\Gamma )^{1/3}. \end{equation}

We recover the scaling law in $Ra_{\phi }^{1/3}$ identified in figure 3(b). We further validate this scaling law in figure 4(b), which shows the compensated r.m.s. horizontal velocity over a large range of $(\Gamma , R_F, Ra_{\phi })$ at $Pr=0.1$. From our DNS data set, we consider the inner part of the cylinder ($r\leq \Gamma /2$) in order to focus on the turbulent bulk, and we compute

(4.6) \begin{equation} U_{h}^{{rms}}=\left \lt \left (u - \langle u\rangle _{\varphi z r_f}\right )^2 + \left (v - \langle v\rangle _{\varphi z r_f}\right )^2 \right \gt _{\varphi z r_f}^{1/2}, \end{equation}

where the $r_f$ subscript indicates an average along the radial direction over $r\in [0,\Gamma /2]$. The agreement with (4.5) is very good with all the data converging at large $Ra_{\phi }$. This is particularly clear for the lower flux ratios. As the value of $R_F$ increases, the relevance of a characteristic horizontal velocity derived from horizontal heat transport becomes less meaningful. In the asymptotic case $R_F=1$, this typical velocity scale effectively disappears since all the injected heat flux is evacuated at the top.

Note that in the context of vertical convection (flow between two differentially heated vertical walls), similar dimensional arguments have been put forward in the seminal work of George & Capp (Reference George and Capp1979). This scaling has been notably revisited and corroborated by subsequent studies such as Versteegh & Nieuwstadt (Reference Versteegh and Nieuwstadt1999); Hölling & Herwig (Reference Hölling and Herwig2005); Ng, Chung & Ooi (Reference Ng, Chung and Ooi2013). Note also that even if we observe the same $1/3$ exponent in $Ra_{\phi }$ as for horizontal convection with imposed flux conditions (Mullarney, Griffiths & Hughes Reference Mullarney, Griffiths and Hughes2004), the velocity corresponds in this case to the horizontal velocity within the thermal boundary layers (Rossby Reference Rossby1965; Griffiths, Hughes & Gayen Reference Griffiths, Hughes and Gayen2013), whereas in our case, it corresponds to the bulk velocity driven by the imposed buoyancy flux. Finally, we pointed out in our previous paper that our system shares similarity with the Grossmann & Lohse (Reference Grossmann and Lohse2000) regime I–l (Rein et al. Reference Rein, Carénini, Fichot, Favier and Le Bars2023a). Hence, one might expect a typical r.m.s. velocity scaling like $Ra_{\phi }^{2/5}$ (i.e. the imposed flux version of the classical $Ra^{1/2}$ scaling in the imposed temperature convection). This does not appear to be consistent with our data, as show in figure 3(b). Presumably, the experimental range of the flux ratio observed ($R_F \in [0.4 , 0.6]$) does not reach the asymptotically high value necessary to manifest the fully developed I–l regime described by Grossmann & Lohse (Reference Grossmann and Lohse2000), which would require larger $Ra_{\phi }$ and/or $R_F$.

4.3. Persistent drifting pattern in the turbulent regime

4.3.1. Characteristics

In Rein et al. (Reference Rein, Carénini, Fichot, Favier and Le Bars2023a) a drifting pattern was observed at high Rayleigh numbers for low flux ratios (see their figure 2(c) for example), but no detail was given about its properties. During the experimental campaign discussed here, the system also showed the spontaneous emergence of such a branch pattern. A video of experiment $3$ from table 1 conducted in air at $Ra_{\phi } \sim 10^7$ can be downloaded as supplementary material available at available at https://doi.org/10.1017/jfm.2025.2. In this video we have removed the average image spanning the entire duration of the footage and applied a moving average over five frames across all images. This adjustment allows for a clearer visualisation of the movement of the PIV particles, represented by white streaks. The video is in real time, except for the clear appearance of a branch where a slowdown was applied. It shows a branch of the pattern, visible as extended white streaks moving from right to left.

Figure 5. $(a)$ Hovmöller (space–time) diagram of the radial velocity averaged along the radial direction. The radial velocity is measured at $z=6$ cm (i.e. $z/H=0.54$). $(b)$ Average of the radial velocity time spectra from each PIV measurement location. The vertical dash–dot line shows the frequency ($f_{{max}} =6.7\times 10^{-2}$ Hz) of the peak amplitude. For both, data correspond to experiment $9$ from table 1 where $Ra_{\phi } = 1.26\times 10^8$, using air.

The drifting branches are also visible through the velocity fields obtained from PIV processing, as, for instance, on the Hovmöller diagram shown in figure 5. Thick bands associated with the highest radial velocities indicate the passage of the branches in the camera’s field of view. The global angular drift frequency of the whole pattern $\omega _p$ can be deduced from such diagrams by tracking a branch of the pattern. The observed drift can occur in both clockwise and anticlockwise directions in independent realisations of the same configuration. This indicates that there is no preferred rotation direction as expected from the problem axisymmetry. A local frequency $\omega$ can also be computed from the temporal spectrum of the radial velocity at each spatial location. The spectra are then averaged over the PIV measurement area for experiments and over the whole domain for DNS. An example based on the PIV measurment from experiment 9 is shown figure 5(b).

Figure 6. $(a)$ Number of branches as a function of $Ra_{\phi }$. The insets show two snapshots of the radial velocity maps from DNS at $z=0.66$ for $Ra_{\phi }=10^5$ and $10^8$, $R_F=0.6$, $Pr=0.7$ and $\Gamma =4$. $(b)$ Angular drift frequency as a function of $Ra_{\phi }$. The dashed–dotted line represents a best fit based on the experimental data for $Ra_\phi \gt 5 \times 10^7$. In both graphs, square symbols ($\square$), diamond symbols ($\Diamond$) and circle symbols ($\circ$) respectively represent data for DNS at $R_F=0.5$ and $R_F=0.6$, and experimental data. Empty circles indicate the use of ${\rm SF}_6$. All DNS were performed with $Pr=0.7$ and $\Gamma =4$.

The two determined frequencies are related by $\omega =m\,\omega _p$, where $m$ is the number of branches or azimuthal wavenumber. Figure 6 shows $m$ and $\omega$ as functions of the Rayleigh number for all our experimental and numerical results. We note a good overall qualitative agreement, even if DNS exhibit a stronger dependence of $m$ on the flux ratio. Regarding $\omega$, a power law emerges, indicating a dependence $\omega \sim Ra_{\phi }^{1/3}$. This relationship was determined through a best fit using the least squares method, applied to the experimental data while excluding values below $Ra_{\phi } = 5\times 10^{7}$ to focus on the asymptotic behaviour at high Rayleigh numbers. This $1/3$ exponent is reminiscent of the one identified for the horizontal velocity in figure 3. However, we argue that both are not directly related. Indeed, from the DNS at $R_F=0.5$ and 0.6, we can compute the mean angular velocity averaged in $\varphi$ and $z$ at $r=3\Gamma /2$ (our conclusion remains qualitatively unchanged for other radii). For the low $Ra_{\phi }$ values of DNS, this advection angular velocity is significantly smaller than the measured angular drift velocity. Hence, we argue that the drifting pattern is not caused by the advection of some structure by a mean azimuthal velocity, but rather seems related to a wave-like phenomenon where the fluid, on average, does not undergo significant azimuthal movement. We also note from our systematic DNS study that this pattern appears only above a given threshold in $Ra_{\phi }$. Initially, prograde and retrograde similar structures superimpose. By increasing $Ra_{\phi }$, one of the two directions is randomly selected and then persists up to the largest $Ra_{\phi }$ investigated. A detailed study of the onset of this oscillatory instability is necessary. This will be the subject of future work, including comparisons with other large-scale structures observed in turbulent convection, such as the jump rope vortex observed in liquid metal by Akashi et al. (Reference Akashi, Yanagisawa, Sakuraba, Schindler, Horn, Vogt and Eckert2022); Cheng et al. (Reference Cheng, Mohammad, Wang, Keogh, Forer and Kelley2022); Teimurazov et al. (Reference Chapman and Cowling2023).

4.3.2. Analysis in the reference frame of the drifting pattern

To gain deeper insights into the temperature and velocity field structure associated with the drifting pattern, we analyse our DNS in the frame of reference rotating with the pattern. To do so, we first determine the angular drift frequency ($\omega$) and azimuthal wavenumber ($m$) from the time and spatial spectra of the azimuthal velocity. Then, we spectrally interpolate, at each time step, all variables on a grid ($X(t), Y(t), Z(t)$) rotating with the phase speed $\omega /m$ of the pattern. We finally introduce a phase-averaging operator denoted by $\langle \cdot \rangle _p$, which is effectively a time average over a time $\tau _p$ within the reference frame of the drifting pattern, defined as

(4.7) \begin{equation} \langle T \rangle _p = \frac {1}{\tau _p}\int ^{t_0+\tau _p}_{t_0} T\left (t,X,Y,Z\right ) \mathrm {d}t \ . \end{equation}

The typical value of $\tau _p$ ranges between 0.2 and 1 diffusive times for the largest and the lowest Rayleigh numbers, respectively. To ensure the convergence of the data obtained from this averaging, we verify that the probability density functions (PDF) of the radial temperature gradient and velocity fields remain consistent even when considering only half of the data.

Figure 7. Phase-averaged temperature, radial velocity and azimuthal velocity fields, processed during 1 diffusive time. Results are presented for three horizontal slices at $z=0.1$, $0.5$ and $0.9$. Parameters are $\Gamma =4$, $R_F=0.1$, $Ra_{\phi }=10^5$ and $Pr=0.1$. In the present DNS, the pattern drifts anticlockwise.

Figure 7 shows an example of the obtained temperature, radial velocity and azimuthal velocity fields at three different depths. Thanks to the phase average, we now clearly identify the six branches pattern on all fields, while it would otherwise disappear when using a standard time average in the laboratory frame. Furthermore, in the upper part of the domain, we observe positive radial velocities moving from the centre towards the edge, while in the lower part, the radial velocities are negative, moving from the edge towards the centre. This aligns with the cooling effect applied at the edge, causing the fluid to cool, increasing its density, and subsequently falling down along the wall.

Concerning the azimuthal velocity, a clear shear is evident between the upper ($z \gt 0.5$) and lower ($z \lt 0.5$) regions. This shear was unnoticed without the phase average, because it was masked by temporal fluctuations of similar magnitude to the average flow (not shown). Note that we observe a shear in the opposite direction compared with figure 7 when the pattern drifts clockwise. In the pattern frame of reference, phase shifts between the fields are observed. The temperature and radial velocity seem to be in phase opposition, while the azimuthal velocity appears to be in phase quadrature with respect to both the temperature and radial velocity (not shown).

We see on the temperature field a well-defined six branches hot pattern at $z=0.1$, which becomes more diffuse as $z$ increases. Yet, for any given depth and radius, large azimuthal fluctuations are observed, including close to the outer wall. This thermal pattern thus plays a crucial role in the transport of the heat flux to the side, which will be the focus of the next section.

4.4. Wall flux fluctuations

In line with Rayleigh–Bénard studies (e.g. Shang et al. (2003, 2004); Gasteuil et al. (Reference Gasteuil, Shew, Gibert, Chillá, Castaing and Pinton2007); Lakkaraju et al. (Reference Lakkaraju, Stevens, Verzicco, Grossmann, Prosperetti, Sun and Lohse2012); Labarre, Fauve & Chibbaro (Reference Labarre, Fauve and Chibbaro2023)), we want to quantify heat flux fluctuations on the lateral wall and to understand their connection with both turbulent fluctuations and the large-scale flow structures. To do so, we define the outgoing flux at the wall $r=\Gamma$ by

(4.8) \begin{equation} \Phi (t,\varphi ,z)=-\frac {\partial T}{\partial r} \bigg |_{r=\Gamma }. \end{equation}

Once the system reaches a statistically stationary state, the global flux balance (by integrating (3.4) over the volume) implies that the flux mismatch between the top and the bottom $(1-R_F)\pi \Gamma ^2$ has to be balanced on time average by the outgoing flux at the side denoted by $\Phi _{{side}}2\pi \Gamma$. Hence,

(4.9) \begin{equation} \Phi _{{side}}=\frac {(1-R_F)\Gamma }{2}. \end{equation}

For the following subsections, the wall heat flux $\Phi$ is normalised by this mean side flux $\Phi _{{side}}$.

4.4.1. The DNS results

In figure 8 we analyse the flux fluctuations in DNS, considering all azimuthal and vertical positions along the side wall over a period of $0.2$ thermal diffusive time. The PDFs of the heat flux are depicted for $Ra_{\phi }$ ranging from near the instability threshold at $Ra_{\phi }=2\times 10^4$ to the fully turbulent regime at $Ra_{\phi }=10^8$. Note that the PDFs remain mostly unchanged when varying the aspect ratio from 4 to 16 at $Ra_{\phi }=10^6$, $Pr=0.1$ and $R_F=0.1$ (not shown).

Figure 8. Wall heat flux PDF from DNS with different values of $Ra_{\phi }$. The fluxes have been scaled by the time-averaged flux $\Phi _{{side}}$. Input parameters are $Pr=0.1$, $R_F=0.1$ and $\Gamma =4$.

With increasing $Ra_{\phi }$, more intense fluctuations across a wider range of values are observed. An exponential tail seems to emerge notably at $Ra_{\phi }=10^8$. It is worth noting that Lakkaraju et al. (Reference Lakkaraju, Stevens, Verzicco, Grossmann, Prosperetti, Sun and Lohse2012) highlighted a similar exponential tail pattern in the heat flux along the sidewall of a cylindrical Rayleigh–Bénard system, with the slope increasing with $Ra_{\phi }$. Additionally, Shang et al. (2003, 2004), in a Cartesian Rayleigh–Bénard configuration, identified a similar exponential behaviour in the distribution of both vertical and horizontal bulk heat fluxes. This exponential tail is not observed close to the onset of instability, but requires reaching sufficiently turbulent conditions.

Figure 9. (a) Three-dimensional view of the temperature field at the bottom plate (note the upside-down axes system with $z$ pointing downward) with $\Gamma = 4$, $Ra_{\phi }=10^7$, $R_F=0.1$ and $Pr=0.1$. (b) Hovmöller (space–time) diagram from this simulation of the wall heat flux averaged in $z$ and scaled by the mean wall flux $\Phi _{{side}}$.

The large-scale drifting pattern is also responsible for intense fluctuations in heat flux near the wall. To illustrate this point, we performed a DNS at $Ra_{\phi }=10^7$, $\Gamma =4$, with $Pr$ = 0.1 and the flux ratio at $R_F=0.1$. Figure 9(a) shows a snapshot of the temperature field, seen from below. The largest temperatures are localised along four drifting branches, as more clearly identified in the space--time diagram in figure 9(b). Within these branches, intense fluctuations involve heat fluxes up to two times the average flux. Once scaled for our air experiment, the characteristic size and duration of one of these fluctuations are approximately 4 cm and 1.5 s, respectively.

Figure 10. $(a)$ The PDF of the phase-averaged wall heat flux as a function of $Ra_{\phi }$. The other input parameters are $Pr=0.1$, $R_F=0.1$ and $\Gamma =4$. All depths and azimuths along the side wall have been considered. $(b)$ Same but for the wall heat flux fluctuations related to the phase average. For both plots, the fluxes have been scaled by the time-averaged flux $\Phi _{{side}}$.

We now want to disentangle the contributions to the side heat flux fluctuations from the large-scale drifting pattern and those from the small-scale turbulent convective fluctuations. Figure 10(a) shows the wall heat flux PDF after applying phase averaging: it thus focuses on the effect of the large-scale drifting pattern only. Note that similar analyses have been conducted in turbulent Rayleigh--Bénard systems in cylindrical geometry by Lakkaraju et al. (Reference Lakkaraju, Stevens, Verzicco, Grossmann, Prosperetti, Sun and Lohse2012) to assess how the local heat flux varies depending on the measurement location relative to the orientation plane of the large-scale circulation. We observe that the flux distribution carried by the pattern extends over a range from 0.5 to 2 times the mean flux $\Phi _{{side}}$. This range barely depends on the Rayleigh number. Indeed, the drifting pattern is essentially responsible for the transport of the mean flux ($\Phi _{{side}}$) out of the system, which is independent of $Ra_{\phi }$. Increasing the $Ra_{\phi }$ increases the rapid turbulent fluctuations, which will be considered next. However, the shape of the PDFs changes with $Ra_{\phi }$. Indeed, the observed pattern evolves as $Ra_{\phi }$ increases, notably by varying the number of branches and increasing the skewness of the flux profile associated with a branch (not shown).

Figure 10(b) shows the PDF of the wall heat flux fluctuations related to the phase-averaged operator: we thus focus here on the turbulent fluctuations taking place around the phase average field. As the Rayleigh number increases, the range of heat flux fluctuations widens. We also see the emergence of the exponential tail at $Ra_{\phi }=10^8$, previously identified in figure 8. This means that the largest and rarest heat flux events are associated with turbulent fluctuations, and not with the pattern itself. Their probability will keep increasing with $Ra_{\phi }$. It is also worth noting that the left portion of the fluctuations, with the amplitude ranging from –2 to 0, appears to be unaffected by changes in the Rayleigh number, showing a consistent overlap. Since $\phi ^\prime = \phi - \langle \phi \rangle _p$ and the heat flux is consistently positive for $R_F=0.1$, negative $\phi ^\prime$ values suggest low heat fluxes $\phi$ and high heat fluxes caused by the pattern’s drift $\langle \phi \rangle _p$. These occurrences indicate scenarios where the observed heat flux primarily arises from the pattern’s drift. This finding aligns with the observation that the pattern’s drift consistently induces similar fluctuations regardless of the Rayleigh number (refer to figure 10a).

In summary, the heat flux fluctuations within the system exhibit a broad distribution that can be attributed to two main contributions: the pattern contribution, which transports the mean flux $\Phi _{{side}}$ with a fixed amplitude independent of $Ra_{\phi }$, and the fluctuations, which exhibit wider distributions and exponential tails as $Ra_{\phi }$ increases. In light of these numerical observations, how does experimental analysis contribute to our understanding of the heat flux fluctuations?

4.4.2. Experimental results

In figure 11 we plot the time evolution of the wall heat flux measured experimentally as well as the radial velocity at a location situated 10 cm from the flux sensor in the radial direction towards the centre. The purpose of this measurement is to explore the relationship between radial velocity and heat flux measured at the wall. Upon examination, we observe that when we shift the radial velocity by a characteristic time $\tau _{{adv}}= {2.6}\,\text {s}$, the two rescaled signals reasonably match. Moreover, the characteristic time $\tau _{{adv}}$ aligns with an advective time scale, namely $d/U$, where $d$ represents the distance between the flux sensor and the location of the velocity measurement and $U$ corresponds to the time-average radial velocity at that specific location. This alignment suggests a correlation between the radial velocity and the wall heat flux. This is compatible with the numerical results shown in figure 9(b), and indicates that the drifting branches associated with higher radial velocity values act as pathways for transporting heat fluctuations towards the wall.

Figure 11. Temporal evolution of the wall heat flux measured by our probe spanning a width of 18 mm and a depth $z\in [{4,6}\,\text {cm}, {6,4}\,\text {cm}]$ (orange), and radial velocity measured by PIV at $r\,=\,{34}\,\text {cm}$, $z\,=\,{6}\,\text {cm}$, in the direction aligned with the flux sensor (green). The time is shifted by 2.6 s for the velocity. For each signal, the mean temporal component has been removed to keep only the fluctuating part, which is then normalised by the maximum value of the signal. Data from experiment 3 in table 1.

Figure 12. (a) Wall heat flux PDF for experiment 9$n^o 9$ in table 1 (orange $\circ$) and for a DNS with $Ra_{\phi }=10^8$, $\Gamma = 4$, $Pr=0.7$, $R_F=0.6$ (green $\Diamond$). The continuous blue thick line represents the low-pass filtered DNS results. The thin red line represents a fit by a gamma distribution. The fluxes are all scaled by the time-averaged flux $\Phi _{{side}}$. (b) Same for experiment 9$n^o 9$ and for DNS heat flux fluctuations due to the drifting pattern only (see § 4.4.1).

In figure 12(a) we compare the PDF of the wall heat flux derived from experimental measurements and from numerical simulations. The heat flux is measured experimentally at mid-height and averaged over the entire surface of the sensor. It is positioned halfway up the sidewall, extending end to end from $z={4.6}\,\text {cm}$ to $z={6.4}\,\text {cm}$, with a width of ${1.8}\,\text {cm}$. We consider experiment 9$n^o 9$ from table 1, where the non-dimensional parameters have been estimated by $Ra_{\phi }=1.26\times 10^8$ and $R_F=0.54$. Numerically, we fix $Ra_{\phi }=10^8$, $\Gamma =4$, $\Pr =0.7$ and consider the flux ratio $R_F=0.6$. Given the numerical difficulty to get statistics over long periods, we assume system ergodicity. As a result, the numerical flux calculations consider all azimuthal events at each output snapshot, limited in height to match the depth of the flux sensor ($z\in [{4.6}\,\text {cm}$ to $ {6.4}\,\text {cm})]$).

When comparing the PDFs, we initially note that the numerical PDF predicts much larger fluctuations than those measured experimentally. However, it should be noted that the sensor has a response time ($T_r$) of 0.7 s, meaning that any event occurring over a period shorter than this will not be directly measured. To account for this experimental aspect, we apply a first-order low-pass temporal filter, with a cutoff frequency ($f_c$) corresponding to the flux sensor’s response time ($f_c=1/T_r$). The application of this filter is done in Fourier space by multiplying the Fourier transform of the experimental signal by the transfer function of a first-order low-pass filter. We then note a good qualitative concordance between DNS and experimental results. A priori, a two-dimensional spatial filter should be considered, since the heat sensor averages over its entire surface. However, the good agreement suggests this is unnecessary, as short events correlate with small-scale phenomena. Note that we did the same numerical computation but for $R_F=0.5$, and the same qualitative agreement was found. As a result, the most extreme events then appear to have a relatively brief duration, at least shorter than the response time associated with the sensor.

The shape of the numerical PDF of figure 12(a) is compatible with a gamma distribution, denoted $\mathcal {P}_{\gamma }[\Phi ]$ and defined by

(4.10) \begin{equation} \mathcal {P}_{\gamma }[\Phi ] = \frac {1}{\gamma (N) \Psi ^N}\,\Phi ^{N-1} e^{-\Phi /\Psi }, \end{equation}

where $\gamma (N)$ represents the gamma function with $N$ and $\Psi$ the two fitting parameters. Here $N$ is an integer without dimensionality, while $\Psi$ has a dimension of a heat flux. Since in our case the mean flux is imposed, $\Psi$ is constrained to $\Phi _{{side}}/N$. The distribution now has only one free parameter, $N$, which is set to 6 in figure 12(a), but depends on the input parameters $Ra_{\phi }$, $R_F$, $Pr$ and $z$ due to vertical inhomogeneity. In figure 8 for example, the PDF accounting for all events along the vertical direction $z$ at $R_F=0.1$ and $Ra_{\phi }=10^8$ is best fitted by a gamma distribution with $N=4$.

In figure 12(b) we analyse the PDF of wall heat flux, comparing experiment 9 with the same DNS conducted at $Pr=0.7$, $\Gamma =4$, $R_F=0.6$ and $Ra_{\phi }=10^8$, but considering only the heat flux distribution due to the drifting pattern. Both PDFs are similar. This makes sense since the drifting pattern is persistent and drifts slowly: it is hence easily catched by the sensor despite its limited response time. The heat sensor however misses the rapid, turbulent fluctuations described in § 4.4.1. The experimental PDF also appears to have fewer values with $\Phi /\Phi _{{side}} \lt 1$ compared with the DNS results. While a perfect match between the two PDFs is not expected, the remaining discrepancy likely stems from the fact that the flux sensor averages over its surface, a factor not accounted for in the numerical data.