2. Oberbeck–Boussinesq approximation

The OB approximation is a simplified model, incorporating buoyancy effects. It assumes an incompressible flow and introduces two key assumptions. Firstly, it posits that all fluid properties remain constant, except for density in the buoyancy force term of the momentum equation, which is assumed to depend linearly on temperature with the isobaric thermal expansion being the relevant coefficient. Secondly, it disregards the contributions of pressure work and viscous dissipation in the heat equation (see also Roche et al. Reference Roche, Gauthier, Kaiser and Salort2010, appendix A2). The OB governing equations for the velocity field $\boldsymbol {u} (\boldsymbol {x}, t)$ (with $u_j$ the velocity components in spatial directions $x_j$), the temperature field $T(\boldsymbol {x}, t)$ and the hydrodynamic pressure $p (\boldsymbol {x}, t)$ include the following continuity equation, and momentum and heat equations:

(2.1)$$\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$

(2.2)$$\begin{gather}\partial_t \boldsymbol{u} + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{u} + {\boldsymbol{\nabla} p/\rho} = \nu\boldsymbol{\nabla}^2 \boldsymbol{u} +\alpha(T-T_0)g\boldsymbol{e}_z, \end{gather}$$

(2.3)$$\begin{gather}\partial_t T + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) T = \kappa\boldsymbol{\nabla}^2 T. \end{gather}$$

Here $\rho$ is the density, $\nu \equiv \mu /\rho$ the kinematic viscosity, $\kappa \equiv k/(\rho c_p)$ the thermal diffusivity, with $\mu$ being the dynamic viscosity, $k$ the thermal conductivity, and $c_p$ the specific heat at constant pressure, and $\partial _t$ denotes the partial derivative in time $t$, and $\boldsymbol {e}_z$ the unit vector pointing upwards. The reference temperature $T_0$ is taken at the arithmetic mean of the top and bottom temperatures, $T_0\equiv (T_++T_-)/2$.

Starting from the works of Spiegel & Veronis (Reference Spiegel and Veronis1960) and Veronis (Reference Veronis1962), the validity of the OB approximation (2.1)–(2.3) has been scrutinized in various theoretical studies, where the magnitudes of specific terms in the governing equations were examined. A rigorous mathematical variational approach, employing expansions of fluid properties as power series, was initially introduced by Mihaljan (Reference Mihaljan1962), who primarily focused on the temperature dependency of the density. Subsequently, other researchers extended the method to incorporate pressure and temperature dependencies of other fluid parameters.

When moving beyond the OB approximation, a simple exploration involves considering a first-order linear dependence on both temperature and pressure, of each material property $\varphi$,

(2.4)\begin{equation} \frac{\varphi }{\varphi_0}=1+\varepsilon_{\varphi,T}\frac{T-T_0}{\varDelta}+\varepsilon_{\varphi,P}\frac{P-P_0}{\rho_0 gL}, \end{equation}

including the density ($\varphi =\rho$), absolute viscosity ($\varphi =\mu$), specific heat at constant pressure ($\varphi =c_p$), thermal expansion coefficient ($\varphi =\alpha$), and thermal conductivity ($\varphi =k$), and $\kappa \equiv k/(\rho c_p)$ is a dependent parameter, fully determined by $k$, $\rho$ and $c_p$ (Gray & Giorgini Reference Gray and Giorgini1976). Thus the manifestation of non-Oberbeck–Boussinesqness (NOBness) is determined by different dimensionless factors $\varepsilon _{\varphi, T}$ and $\varepsilon _{\varphi, p}$. When dealing with a specific fluid under given operational conditions of reference temperature $T_0$ and pressure $P_0$, it is not immediately evident which among these parameters $\varepsilon _{\varphi,T}$ and $\varepsilon _{\varphi,p}$ is most relevant in inducing NOBness in the experimental set-up. (The fluid properties at the reference $T_0$ and $P_0$ will be denoted in the following with a subscript ‘$0$’.) Unravelling the diverse NOB effects is possible in direct numerical simulations (DNS) which allow for the use of artificial fluids that exhibit only a single specific source of NOBness. For instance, Ahlers et al. (Reference Ahlers, Calzavarini, Araujo, Funfschilling, Grossmann, Lohse and Sugiyama2008) conducted such an investigation for ethane. Various studies have explored NOB effects in RB convection in cryogenic helium and pressurized SF$_6$ (Roche et al. Reference Roche, Gauthier, Kaiser and Salort2010; Shishkina, Weiss & Bodenschatz Reference Shishkina, Weiss and Bodenschatz2016; Weiss et al. Reference Weiss, He, Ahlers, Bodenschatz and Shishkina2018; Roche Reference Roche2020; Yik, Valori & Weiss Reference Yik, Valori and Weiss2020), as well as in water and glycerol (Manga & Weeraratne Reference Manga and Weeraratne1999; Ahlers et al. Reference Ahlers, Brown, Araujo, Funfschilling, Grossmann and Lohse2006; Sugiyama et al. Reference Sugiyama, Calzavarini, Grossmann and Lohse2009; Horn, Shishkina & Wagner Reference Horn, Shishkina and Wagner2013; Horn & Shishkina Reference Horn and Shishkina2014).

In experiments, the NOB effects are always present (although they can be negligible) and currently, there is no established method for precisely estimating their influence in heat transport measurements. However, it is possible to estimate a priori the values of $\varepsilon _{\varphi,T}$ and $\varepsilon _{\varphi,p}$ and set thresholds for their maximal values in an experiment. Note that in some studies only the temperature variation of the density ($\alpha \varDelta$) is considered as a measure of NOBness (e.g. Niemela & Sreenivasan Reference Niemela and Sreenivasan2003); however, it is only one ($\varepsilon _{\rho,T}$) of in total ten sources of NOBness (all parameters $\varepsilon _{\varphi,T}$ and $\varepsilon _{\varphi,P}$).

To study the validity of the OB approximation, we follow a comprehensive variational approach suggested by Gray & Giorgini (Reference Gray and Giorgini1976). We start with the governing equations for a flow of a Newtonian fluid of variable properties (§ 15, Chapter II of Landau & Lifshitz Reference Landau and Lifshitz1987):

(2.5)$$\begin{gather} D_t \rho + \rho \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$

(2.6)$$\begin{gather}\rho D_t \boldsymbol{u} + \boldsymbol{\nabla} P= \boldsymbol{\nabla}\boldsymbol{\cdot}(\mu {\boldsymbol{\mathsf{S}}}) -\rho g \boldsymbol{e}_z +\boldsymbol{\nabla}(\eta\,\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}), \end{gather}$$

(2.7)$$\begin{gather}\rho c_p D_t T= \boldsymbol{\nabla}\boldsymbol{\cdot}\left(k\boldsymbol{\nabla} T\right) +\alpha T D_tP +\mu \varPhi, \end{gather}$$

where $P=p-\rho _0gz\boldsymbol {e}_z$ is the (full) pressure, $\boldsymbol{\mathsf{S}}$ the deformation rate tensor with components, $S_{ij}\equiv \partial _j u_i+\partial _i u_j -(2\delta _{ij}/3)\partial _k u_k$, $\varPhi \equiv (S_{ij}/2)(\partial _j u_i+\partial _i u_j)$ the viscous dissipation function, $D_t\equiv \partial _t+\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }$ denotes the full (material) derivative, $\delta _{ij}$ the Kronecker symbol, and $\eta$ the second viscosity. As it has been shown by Spiegel & Veronis (Reference Spiegel and Veronis1960), the last term in (2.6) can be neglected if pressure fluctuations are smaller then static variation, which is the case for relatively slow convective flows. Therefore, we assume that the last term in (2.6), which is associated with the second viscosity $\eta$, vanishes.

Substituting the linear representations (2.4) of the fluid properties $\varphi =\rho$, $\mu$, $c_p$, $\alpha$ and $k$ into (2.5)–(2.7), one requires that the residuals, i.e. the terms that distinguish the resulting equations from their OB approximation (2.1)–(2.3), are negligibly small. The OB requirements are fulfilled if not only all $\varepsilon _{\varphi,T}$ and $\varepsilon _{\varphi,P}$ are negligibly small, but also the magnitudes of the pressure work term ($\alpha T\, D_tP$) and of the dissipation term ($\mu \varPhi$) in (2.7) are negligibly small compared with the magnitudes of the other terms in (2.7). Comparing $\rho c_p D_t T$ with $\alpha T\, D_tP$ in (2.7), we conclude that the pressure work is negligible if $\rho _0 c_{p,0} \varDelta \gg \alpha _0T_0\rho _0gL$ (here we estimate the pressure magnitude as $P\sim \rho _0gL$ and the magnitudes of the temperature variation as $\varDelta$ and of the absolute value of the temperature as $T_0$). Comparing $\boldsymbol {\nabla }\boldsymbol {\cdot }(k\,\boldsymbol {\nabla } T)$ with $\mu \varPhi$, we conclude that the dissipation term is negligible if $k_0\varDelta \gg \mu _0\alpha _0 g \varDelta L$ (here the velocity magnitude is estimated as the free-fall velocity $\sqrt {\alpha _0gL\varDelta }$). The last two inequalities can be reformulated, respectively, as $bT_0/\varDelta \ll 1$ and $b{{Pr}}_0\ll 1$, where ${{Pr}}_0\equiv \nu _0/\kappa _0$ is the reference Prandtl number, $\kappa _0\equiv k_0/(\rho _0\,c_{p, 0})$ is the reference thermal diffusivity and $b\equiv \alpha _0gL/c_{p,0}$.

Introducing a certain small threshold on the degree of NOBness $\hat {\sigma }$, $0<\hat {\sigma }\ll 1$, we say that the OB approximation is valid with the accuracy $\hat {\sigma }$, if requirements (2.8)–(2.9) are fulfilled:

(2.8)$$\begin{gather} \varepsilon_{\varphi,T}<\hat{\sigma}, \quad \varepsilon_{\varphi,P}<\hat{\sigma},\quad \text{for } \varphi=\rho, \mu, c_p, k, \end{gather}$$

(2.9)$$\begin{gather}bT_0/\varDelta<\hat{\sigma},\quad b{{Pr}}_0<\hat{\sigma}. \end{gather}$$

Relations (2.8)–(2.9) are derived by substituting (2.4) for all fluid properties $\varphi$ into (2.5)–(2.7) and requiring that the terms, which are not present in the OB approximation (2.1)–(2.3), become negligibly small when the measure of the NOBness $\hat {\sigma }$ becomes infinitesimal, $\hat {\sigma }\rightarrow 0$.

3. The validity of the Oberbeck–Boussinesq approximation for different set-ups

For any given fluid and threshold $\hat {\sigma }$ for the residuals, from (2.8)–(2.9) one can derive the OB-validity region in terms of the maximally possible temperature difference $\varDelta$ in the system and height $L$ of the container. Here, to capture the NOB effects, we confine our examination to the lowest order, (2.4), focusing solely on a linear term in the temperature and pressure expansions of the fluid properties. This is sufficient, as shown in a recent study by Macek et al. (Reference Macek, Zinchenko, Musilova, Urban and Schumacher2023) for the case of cryogenic helium.

In figure 1 the OB-regions are calculated for operational conditions typically employed in RB experiments and some fluids, namely, for water, air, ethane, helium, pressurized gas sulphur hexafluoride (SF$_6$) at room temperatures, and cryogenic helium, using REFPROP (2013). For any given fluid and reference temperature $T_0$ and pressure $P_0$, the OB-validity regions depend on the OB threshold $\hat {\sigma }$. The choices $\hat {\sigma }=5\,\%$, $\hat {\sigma }=10\,\%$ and $\hat {\sigma }=20\,\%$ give the embedded OB-validity domains coloured, respectively, green, blue and red in figure 1 (see also sketch in figure 4(a) in Ecke & Shishkina Reference Ecke and Shishkina2023). There, each triangularly shaped OB-validity domain is bounded from the right and left by maximum admissible variations of the fluid properties with the temperature (vertical line) and by maximum admissible pressure work term in the heat equation (inclined line), respectively. The lower boundary is determined by the onset of convection, which is calculated in figure 1 for an infinite horizontally extended fluid layer, such that ${{Ra}}_c\approx 1708$. For any laterally bounded domain, this boundary moves up, since the critical ${{Ra}}_c$ for the onset of convection scales as $\sim \varGamma ^{-4}$ for $\varGamma \rightarrow 0$ (Shishkina Reference Shishkina2021; Ahlers et al. Reference Ahlers2022).

Figure 1. Regions of the validity of the OB approximation, in terms of the maximum temperature difference, $\varDelta$, and container height, $L$, according to (2.8)–(2.9), for different fluids: (a) water at $T_0=40\,^\circ {\rm C}$ and $P_0=1$ bar, (b) air at $T_0=40\,^\circ {\rm C}$ and $P_0=1$ bar, (c) ethane at $T_0=40\,^\circ {\rm C}$ and $P_0=1$ bar, (d) helium at $T_0=40\,^\circ {\rm C}$ and $P_0=1$ bar, (e) helium at $T_0=-268.15\,^\circ {\rm C}=5$ K and $P_0=1$ bar, (f) SF$_6$ at $T_0=30\,^\circ {\rm C}$ and $P_0=20$ bar. The nested green, blue and red OB-validity regions correspond, respectively, to the thresholds on the degree of NOBness $\hat {\sigma }=5\,\%$, $10\,\%$ and $20\,\%$. The lower boundaries of these regions have the slopes $L\propto \varDelta ^{-1/3}$, the left boundaries $L\propto \varDelta$, and the right ones are vertical. The values of the maximum achievable Rayleigh numbers, ${{Ra}}_{max,\hat {\sigma }}$, are written with the corresponding colours.

The boundedness of the OB-validity region restricts the maximal ${{Ra}}$, which can be reached in almost-OB experiments. This means that for any chosen fluid and threshold on the degree of NOBness (parameter $\hat {\sigma }$), Rayleigh numbers ${{Ra}}$ larger than a certain maximum value ${{Ra}}_{max,\,\hat {\sigma }}$ can in principle not be realized experimentally. In figure 1, these ${{Ra}}_{max,\hat {\sigma }}$ values are achieved in the upper corners of each OB-validity region, and each value of ${{Ra}}_{max,\hat {\sigma }}$ is written in the plots with the colour of the corresponding OB-validity region.

By comparing the values of ${{Ra}}_{max, \hat {\sigma }}$ among the analysed fluids in figure 1, it can be inferred that water emerges as the optimal fluid for investigating RB convection under OB conditions at extreme ${{Ra}}$, owing to its highest ${{Ra}}_{max, \hat {\sigma }}$. Nonetheless, we note that achieving such exceedingly high ${{Ra}}$ is feasible only when the water layer's depth is several hundred metres. Considering this limitation, the utilization of cryogenic helium and pressurized SF$_6$ becomes more advantageous for experimental inquiries into the transition to the ultimate regime. In figure 2, by examples of pressurized SF$_6$ (used in Göttingen experiments) and cryogenic helium (used in other studies of the ultimate regime, including Oregon, Grenoble and Brno) we show the dependencies of ${{Ra}}_{max, \hat {\sigma }}$ on the reference temperature $T_0$ and pressure $P_0$. One sees, in particular, that keeping the operating pressure $P_0$ in a range between 10 bar and 15 bar in experiments with SF$_6$ is favourable in the sense that this allows higher values of ${{Ra}}$ under (almost) OB conditions to be achieved. In all studied cases, we found that pressure variations of the fluid properties are negligible compared with their temperature variations.

Figure 2. Maximum achievable Rayleigh numbers ${{Ra}}_{max, \hat {\sigma }}$, as functions of the reference pressure $P_0$, for different thresholds on the degree of NOBness: $\hat {\sigma }=5\,\%$ (green), $\hat {\sigma }=10\,\%$ (blue) and $\hat {\sigma }=20\,\%$ (red), for (a) gaseous SF$_6$ (solid lines for $T_0=20\,^\circ {\rm C}$ and dashed lines for $T_0=30\,^\circ {\rm C}$) and (b) cryogenic helium (thick lines for $T_0=2.5$ K, thin lines for $T_0=4.0$ K and dashed lines for $T_0=5.5$ K). The continuous lines show the liquid phase only, while the dashed lines correspond to both phases, gas and liquid, showing a V-shape near the critical pressure $P_c\approx 2.27$ bar, where the fluid properties become very sensitive to the temperature and pressure variations.

4. The effect of the temperature dependency of $c_p$ on the measured ${{Nu}}$

Among all fluid properties, the specific heat capacity $c_p$ was found to have the strongest variations with the temperature (i.e. $\varepsilon _{c_p,T}$ is significantly larger than any other $\varepsilon _{\varphi,T}$ or $\varepsilon _{\varphi,P}$), for both considered fluids, which are pressurized SF$_6$ and cryogenic helium.

For each of the two considered fluids, we selected two cases (characterized by $P_0$ and $[T_-,T_+]$) that correspond to some representative Göttingen or Brno measurements: one case is the ‘classical’ one, where the heat transport scaling exponent $\gamma <1/3$ was measured, and the other case we call ‘ultimate’, where $\gamma \gtrsim 0.4$ was measured, see figure 3. In the classical cases, $c_p(T)$ is almost flat, while in the ultimate cases it rapidly decreases with increasing temperature, for both, SF$_6$ and helium. By example of these four cases (classical and ultimate, for both, pressurized SF$_6$ and cryogenic helium) we want to investigate in DNS whether the strong temperature variation of $c_p$ alone (as in the ultimate cases) can significantly increase the heat transport and with this alter the scaling exponent $\gamma$ in the ${{Nu}}\sim {{Ra}}^\gamma$ relation.

Figure 3. (a,c) Specific heat capacity $c_p$ as functions of the temperature $T$, in the classical case (blue symbols) and ultimate case (red symbols), for (a) SF$_6$ for $T\in [13.49\,^\circ {\rm C}, 29.49\,^\circ {\rm C}]$ and pressure $P_0=15.28$ bar (classical) and $P_0=17.68$ bar (ultimate), and for (c) cryogenic helium for $T\in [4.305\,\text {K},\,4.409\,\text {K}]$ and $P_0=0.8137$ bar (classical), and $T\in [5.097\,\text {K},\,5.214\,\text {K}]$ and $P_0=2.0724$ bar (ultimate). (b,d) Non-dimensionalized data from (a,c), respectively.

Assuming that for any small $\hat {\sigma }>0$, all requirements (2.8)–(2.9) are fulfilled apart from the variation of $c_p$ with the temperature, so that $\varphi _{c_p,T}$ is not negligible, we obtain the governing equations, which are similar to the OB equations (2.1)–(2.3), but with

(4.1)\begin{equation} (c_p(T)/c_{p, 0})[ \partial_t T + (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) T ] = \kappa_0\boldsymbol{\nabla}^2 T \end{equation}

instead of (2.3). (In terms of independent fluid parameters, (4.1) can be rewritten as $c_p(T)[ \partial _t T + (\boldsymbol {u}\boldsymbol {\cdot }\boldsymbol {\nabla }) T ] = (k_0/\rho _0)\boldsymbol {\nabla }^2 T$.)

The four profiles of $c_p(T)/c_{p, 0}$ (classical and ultimate, for pressurized SF$_6$ and cryogenic helium, cf. figure 3b,d) are approximated by polynomials

(4.2)\begin{equation} \frac{c_p(T) }{c_{p, 0}}=1+\sum_{j=1}^{3}c_j\left(\frac{T-T_0}{\varDelta}\right)^j, \end{equation}

with $c_j$ being the respective coefficients. This together with (2.1)–(2.2) and (4.1) leads to the Nusselt number

(4.3)\begin{equation} {{Nu}}=\frac{L}{\kappa_0\varDelta}\left\langle u_z(T-T_0)\left( 1+\sum_{j=1}^{3}\frac{c_j} {j+1}\left(\frac{T-T_0}{\varDelta}\right)^j \right)-\kappa_0\frac{\partial T}{\partial z}\right\rangle_{S_z,t}, \end{equation}

where $u_z$ is a component of the velocity field in the vertical direction $z$, and $\langle \cdots \rangle _{S_z,t}$ means the averaging in time and over any horizontal cross-section $S_z$.

We conducted DNS of RB convection in a cubic domain, for ${{Ra}}=10^6$, $10^7$, $10^8$, $10^9$ and $10^{10}$, and Prandtl number ${{Pr}}=1$, according to (2.1)–(2.2) and (4.1) and for the four temperature profiles of $c_p(T)/c_p(T_0)$ as in figure 3(b,d), to check whether the ultimate profiles of $c_p(T)$ can significantly increase the Nusselt numbers compared with the classical profiles. For this purpose, we used the computational code goldfish (Shishkina et al. Reference Shishkina, Horn, Wagner and Ching2015; Reiter & Shishkina Reference Reiter and Shishkina2020; Reiter, Zhang & Shishkina Reference Reiter, Zhang and Shishkina2022) with sufficiently fine computational grids (Shishkina et al. Reference Shishkina, Stevens, Grossmann and Lohse2010; Kooij et al. Reference Kooij, Botchev, Frederix, Geurts, Horn, Lohse, van der Poel, Shishkina, Stevens and Verzicco2018) that have at least 14 grid points in each boundary layer and $128^3$, $192^3$, $256^3$, $360^3$ and $576^3$ grid points in total for ${{Ra}}=10^6$, $10^7$, $10^8$, $10^9$ and $10^{10}$, respectively.

We expected that the classical profiles of $c_p(T)$ would lead to results close to those from the OB approximation, but the ultimate profiles would lead to visible NOB effects, in particular, to clear differences in the heat transport. Surprisingly enough we failed to recognise any NOB effects or any tendency associated with the NOBness that increases with growing ${{Ra}}$; we found it neither for the classical profiles, nor for the ultimate ones, see table 1. Of course, different DNS show slightly different values of ${{Nu}}$, but there is no pattern that indicates the growth of the scaling exponent $\gamma$ for the ultimate profiles compared with the classical ones.

Table 1. Nusselt numbers $Nu$ obtained for different $Ra$ and fluids (pressurized SF$_6$ and cryogenic helium), and the classical or ultimate profiles of $c_p$, see the main text. Statistics is collected over a duration of at least 800 free-fall time units and at least 5 samples per time unit. The relative difference $d$ is defined as $d\equiv (Nu^{\text {`}ultimate\text {'}}- Nu_{OB})/Nu_{OB}$. The error $\sigma$ (standard deviation of $Nu(z)$ profile) is calculated as $\sigma ^2 =({1}/({N_z-1})) {\sum }_{n=1}^{N_z} (Nu -\langle Nu\rangle _{S_z,t}) ^2$, where $N_z$ is the total number of grid points in the vertical direction $z$.