Hostname: page-component-848d4c4894-ndmmz Total loading time: 0 Render date: 2024-05-19T22:13:14.670Z Has data issue: false hasContentIssue false

Towards the ultimate regime in Rayleigh–Darcy convection

Published online by Cambridge University Press:  02 February 2021

Sergio Pirozzoli
Affiliation:
Dipartimento di Ingegneria Meccanica e Aerospaziale, Sapienza Università di Roma, 00184Rome, Italy
Marco De Paoli
Affiliation:
Institute of Fluid Mechanics and Heat Transfer, TU Wien, 1060Vienna, Austria
Francesco Zonta
Affiliation:
Institute of Fluid Mechanics and Heat Transfer, TU Wien, 1060Vienna, Austria
Alfredo Soldati*
Affiliation:
Institute of Fluid Mechanics and Heat Transfer, TU Wien, 1060Vienna, Austria Polytechnic Department, Università degli Studi di Udine, 33100Udine, Italy
*
Email address for correspondence: alfredo.soldati@tuwien.ac.at

Abstract

Numerical simulations are used to probe Rayleigh–Darcy convection in fluid-saturated porous media towards the ultimate regime. The present three-dimensional dataset, up to Rayleigh–Darcy number $\textit {Ra}=8\times 10^4$, suggests that the appropriate scaling of the Nusselt number is $\textit {Nu}=0.0081\textit {Ra}+0.067\textit {Ra}^{0.61}$, fitting the computed data for $\textit {Ra}\gtrsim 10^3$. Extrapolation of current predictions to the ultimate linear regime yields the asymptotic law $\textit {Nu}=0.0081 \textit {Ra}$, about $16\,\%$ less than indicated in previous studies. Upon examination of the flow structures near the boundaries, we confirm previous indications of small flow cells hierarchically nesting into supercells, and we show evidence that the supercells at the boundary are the footprints of the megaplumes that dominate the interior part of the flow. The present findings pave the way for more accurate modelling of geophysical systems, with special reference to geological $\textrm {CO}_2$ sequestration.

Type
JFM Rapids
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© The Author(s), 2021. Published by Cambridge University Press.

1. Introduction

When a fluid-saturated porous layer undergoes bottom-up heating, it is prone to exhibit Rayleigh–Darcy convection (Horton & Rogers Reference Horton and Rogers1945; Lapwood Reference Lapwood1948; Graham & Steen Reference Graham and Steen1994). The main controlling parameter is the Rayleigh–Darcy number $\textit {Ra}$, expressing the ratio of buoyancy to dissipative forces, and the resulting heat flux is characterized in terms of the Nusselt number, $\textit {Nu}$. The scaling of $\textit {Nu}$ with $\textit {Ra}$ in natural convection is traditionally described via power-law dependency, $\textit {Nu}\sim \textit {Ra}^{\alpha }$, although $\alpha$ may change as a result of the change in the structure of the convective flow. For instance, in Rayleigh–Bénard convection (He et al. Reference He, Funfschilling, Nobach, Bodenschatz and Ahlers2012; Lepot, Aumaître & Gallet Reference Lepot, Aumaître and Gallet2018; Zhu et al. Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018; Wang, Zhou & Sun Reference Wang, Zhou and Sun2020), the general agreement is that $\alpha$ should gradually increase from about $1/3$ at moderate $\textit {Ra}$, to asymptotically approach a value of $1/2$, corresponding to the attainment of an ‘ultimate’ convection regime (Kraichnan Reference Kraichnan1962).

Classical theory suggests that the appropriate scaling in the ultimate regime of Rayleigh–Darcy convection should be linear (Malkus Reference Malkus1954; Howard Reference Howard1964). The main argument is that at high $\textit {Ra}$, the interior part of the flow is perfectly mixed, and temperature gradients are confined to near-wall boundary layers whose thickness is $\delta \propto \textit {Ra}^{-1}$ (whereas in Rayleigh–Bénard convection $\delta \propto \textit {Ra}^{-1/3}$); hence $\textit {Nu} \propto \delta ^{-1} \propto \textit {Ra}$. Other more rigorous arguments have also been provided based on variational calculus, which show that linear scaling yields maximum heat transfer among all suitably energy-constrained flows (Otero et al. Reference Otero, Dontcheva, Johnston, Worthing, Kurganov, Petrova and Doering2004; Wen et al. Reference Wen, Dianati, Lunasin, Chini and Doering2012; Hassanzadeh, Chini & Doering Reference Hassanzadeh, Chini and Doering2014). This prediction is supported by existing numerical studies, mainly based on two-dimensional computations (Otero et al. Reference Otero, Dontcheva, Johnston, Worthing, Kurganov, Petrova and Doering2004; Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2012; Wen et al. Reference Wen, Dianati, Lunasin, Chini and Doering2012; Wen, Corson & Chini Reference Wen, Corson and Chini2015; De Paoli, Zonta & Soldati Reference De Paoli, Zonta and Soldati2016) carried out in the range $\textit {Ra} \leqslant 10^5$. Three-dimensional simulations carried out in recent times up to $\textit {Ra}=2\times 10^4$ (Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2014) also seem to support the establishment of a shifted linear variation of $\textit {Nu}$ with $\textit {Ra}$. However, given the limited range of $\textit {Ra}$ in which data are available, uncertainties remain regarding the actual establishment of the expected ultimate convection regime, and especially how it is reached starting from finite $\textit {Ra}$ (Hewitt Reference Hewitt2020).

The quest for the $\textit {Nu}$ scaling in the high-$\textit {Ra}$ regime is crucial from a fundamental and also from an applied viewpoint, as the Rayleigh–Darcy model closely mimics flows in a number of geophysical applications, with special reference to geological $\textrm {CO}_2$ sequestration in deep saline aquifers (Hidalgo et al. Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012; Huppert & Neufeld Reference Huppert and Neufeld2014; Riaz & Cinar Reference Riaz and Cinar2014; Emami-Meybodi et al. Reference Emami-Meybodi, Hassanzadeh, Green and Ennis-King2015; De Paoli Reference De Paoli2021). Real-life instances may exhibit $\textit {Ra}$ up to ${\sim }{ {O}}(10^5 \div 10^6)$, and small deviations from the expected ultimate linear scaling can produce large differences in the overall predicted transfer flux and in the corresponding cumulative time integral (Slim Reference Slim2014; De Paoli et al. Reference De Paoli, Zonta and Soldati2016; De Paoli, Zonta & Soldati Reference De Paoli, Zonta and Soldati2017). The goal of the present work is to investigate the high-$\textit {Ra}$ range of Rayleigh–Darcy convection, using a database of three-dimensional numerical simulations up to $\textit {Ra}=8\times 10^4$. Our aim is to extract reliable trends for the Nusselt number in a wide range of $\textit {Ra}$, and possibly to infer robust asymptotic estimates for the ultimate linear regime.

2. Methodology

We consider a fluid-saturated porous medium in a three-dimensional domain (see figure 1) with uniform porosity $\phi$ and uniform permeability $\kappa$. The flow, which is incompressible and governed by the Darcy law, is characterized by an unstable density difference ($\Delta \rho ^{*}$) induced by heating the flow from the bottom and cooling it from the top. Indicating by $x^*,z^*$ the horizontal directions, by $y^*$ the vertical direction (along which gravity acceleration $\boldsymbol {g}$ is directed) and by $\theta ^{*}$ the fluid temperature, the physical instance we consider corresponds to a fluid-saturated cubic volume with side length $h^{*}$ heated at the bottom, $\theta ^*(y^{*}=0)=\theta ^{*}_{max}$, and cooled at the top, $\theta ^*(y^{*}=h^{*})=\theta ^{*}_{min}$. The evolution of the temperature field is controlled by the advection–diffusion equation

(2.1)\begin{equation} \phi\frac{\partial \theta^{*}}{\partial t^{*}}+\boldsymbol{\nabla} \boldsymbol{\cdot}(\boldsymbol{u}^{*}\theta^{*}-\phi D \boldsymbol{\nabla}\theta^{*})=0 { ,} \end{equation}

where $t^{*}$ is time, $\boldsymbol {u}^{*}=(u^{*},v^{*},w^{*})$ is the volume-averaged velocity field and $D$ is the thermal diffusivity, which is considered constant here. The superscript $^{*}$ is used to indicate dimensional variables. We assume that fluid density, $\rho ^{*}$, is a linear function of temperature,

(2.2)\begin{equation} \rho^{*}(\theta^{*})=\rho^{*}(\theta^{*}_{min})-\Delta\rho^{*}\frac{\theta^{*}-\theta^{*}_{min}}{\theta^{*}_{max}- \theta^{*}_{min}}{ ,} \end{equation}

with $\Delta \rho ^{*}=\rho ^{*}(\theta ^{*}_{min})-\rho ^{*}(\theta ^{*}_{max})$. Assuming validity of the Boussinesq approximation (Landman & Schotting Reference Landman and Schotting2007; Zonta & Soldati Reference Zonta and Soldati2018), the flow field is fully described by the continuity and Darcy equations

(2.3a,b)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^{*}=0,\quad \boldsymbol{u}^{*}=-\frac{\kappa}{\mu}\left(\boldsymbol{\nabla}P^{*}+\rho^{*}g \boldsymbol{j}\right) { ,} \end{equation}

with $\mu$ the fluid viscosity (constant), $P^{*}$ the pressure and $\boldsymbol {j}$ the vertical unit vector. The walls are assumed to be impermeable and isothermal, and periodicity is assumed in the wall-parallel directions.

Figure 1. Sketch of the computational domain – a cube with side length $h^*$ – used to study Rayleigh–Darcy convection. The flow is heated at the bottom, $\theta ^*(y^*=0)=\theta ^*_{{max}}$, and cooled at the top, $\theta ^*(y^*=h^*)=\theta ^*_{{min}}$, and boundaries in the $x^*$ and $z^*$ directions are assumed to be periodic. The gravity acceleration ($\boldsymbol {g}$) points downwards. The temperature distribution $\theta ^*$ for the case $\textit {Ra}=1\times 10^4$ is also shown for illustrative purposes on the side boundaries and in a plane very close to the top boundary (i.e. at a distance of $50h^*/\textit {Ra}$ from the top boundary).

2.1. Dimensionless equations

Natural velocity and length scales for the system are the buoyancy velocity, $V^{*}=g \Delta \rho ^{*} \kappa / \mu$, and the domain height, $h^{*}$, respectively. Using the set of dimensionless variables

(2.4ac)\begin{equation} \theta=\frac{\theta^{*}-\theta^{*}_{min}}{\theta^{*}_{max}-\theta^{*}_{min}} ,\quad t=\frac{t^{*}}{\phi h^{*}/V^{*}},\quad P=\frac{P^{*}}{\Delta \rho^{*}gh^{*}} , \end{equation}

and introducing the reduced pressure $p^{*}$, we obtain the dimensionless form of the governing equations (2.1), (2.3a,b):

(2.5)\begin{gather} \frac{\partial \theta}{\partial t}+\boldsymbol{\nabla} \boldsymbol{\cdot} \left( {\boldsymbol{u}} \theta -\frac{1}{\textit{Ra}} \boldsymbol{\nabla} \theta\right)=0 ,\end{gather}
(2.6a,b)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0,\quad {\boldsymbol{u}}=-\left( \boldsymbol{\nabla} p -\theta \boldsymbol{j}\right) ,\end{gather}

where $\textit {Ra}=g \Delta \rho ^{*} \kappa h^{*} / (\phi D \mu )=V^{*} h^{*} / (\phi D)$ is the Rayleigh–Darcy number. The wall boundary conditions for velocity and temperature then read as

(2.7a)\begin{gather} v(y=0)=0,\quad \theta(y=0)=1, \end{gather}
(2.7b)\begin{gather}v(y=1)=0,\quad \theta(y=1)=0. \end{gather}

As previously mentioned, for the physical system investigated here, the buoyancy velocity $V^{*}$ is a natural reference velocity scale (Fu, Cueto-Felgueroso & Juanes Reference Fu, Cueto-Felgueroso and Juanes2013; Wen et al. Reference Wen, Akhbari, Zhang and Hesse2018). At the same time, a possible reference length scale is the thickness of the porous layer, ${x}_{c}^{*}=h^*$ (convective scaling). However, an alternative choice for the reference length scale, which is perhaps more related to the physics of the phenomena under investigation, is ${x}_{d}^{*}=\phi D/V^{*}$ (diffusive–convective scaling). Note that ${x}^{*}_{d}$ represents the length over which advection and diffusion balance (Slim Reference Slim2014), and is independent of the physical domain thickness. When rescaled in the latter way, dimensions are bounded within the range ${x}^{*}/{x}_{d}^{*}\in [0,\textit {Ra}]$, and comparison between simulations at different $\textit {Ra}$ is easier. For this reason, lengths in this paper are rescaled with respect to $x_{d}^{*}$. Furthermore, introduction of this length scale also yields another interpretation of the Rayleigh–Darcy number, $\textit {Ra}={x}_{c}^{*}/{x}_{d}^{*}$, which may be regarded as the dimensionless height of the domain (Slim Reference Slim2014).

2.2. Computational details

The numerical simulations rely on the modified version of a second-order finite-difference incompressible flow solver, based on staggered arrangement of the flow variables (Orlandi Reference Orlandi2000), which has been extensively used for DNS of wall-bounded neutrally-buoyant and unstably-stratified turbulent flows (Pirozzoli Reference Pirozzoli2014; Pirozzoli et al. Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017). The temperature transport equation is advanced in time by means of a hybrid third-order low-storage Runge–Kutta algorithm, whereby the convective terms are handled explicitly and the diffusive terms are handled implicitly, limited to the wall-normal direction. This approach guarantees that the total temperature variance is discretely preserved in the limit of inviscid flow. The pressure field in the forced Darcy equation is determined by solving the Poisson equation resulting from the divergence-free constraint,

(2.8)\begin{equation} \nabla^2 p = \frac{\partial \theta}{\partial y}, \end{equation}

with $\partial p / \partial y = 0$ at walls to satisfy the impermeability condition, which is then fed into (2.6a,b). Fourier expansion along the horizontal periodic directions yields a system of tridiagonal equations in the wall-normal direction for each Fourier mode, which is then solved with standard highly efficient numerical techniques (Kim & Moin Reference Kim and Moin1985; Orlandi Reference Orlandi2000).

The mesh spacing in the wall-parallel directions was decided based on preliminary grid-resolution studies at low $\textit {Ra}$ and inspection of the temperature spectra, to prevent any energy pile-up at the smallest resolved flow scales. Regarding the resolution in the wall-normal direction, we have followed the criterion that twenty points should be placed within the thermal boundary layer edge, identified through the peak location of the temperature variance, and grid points are clustered towards the walls using a hyperbolic tangent stretching function. Given the expected linear growth of the temperature gradients, the number of points in each coordinate direction was increased proportionally to $\textit {Ra}$. The time step is selected so that the CFL number is approximately unity for all the simulations herein reported. Preliminary calculations, carried out at $\textit {Ra} \leqslant 5\times 10^3$, have shown excellent agreement with the results of Hewitt et al. (Reference Hewitt, Neufeld and Lister2014).

3. Results

3.1. Heat transfer

We have carried out numerical simulations in the Rayleigh number range $2.5\times 10^3<\textit {Ra}<8\times 10^4$. Note that, because of the time step restrictions, the simulation at $\textit {Ra}=8\times 10^4$ – using $6144\times 6144\times 2048$ grid points – is almost 4096 times more computationally expensive than the simulation at $\textit {Ra}=1\times 10^4$ – using $768\times 768\times 256$ grid points – hence the computational challenge is substantial. A summary of the main computational parameters is provided in table 1. The heat flux is quantified in terms of the Nusselt number $\textit {Nu}$, evaluated as the mean temperature gradient at the top and bottom boundaries,

(3.1)\begin{equation} \textit{Nu}=-\left \langle \frac{1}{2}\int_{0}^{1}\int_{0}^{1}\left.\frac{\partial \theta}{\partial y}\right\vert_{y=0}+\left.\frac{\partial \theta}{\partial y}\right\vert_{y=1} \,\text{d}x\,\text{d}z\right\rangle, \end{equation}

with angle brackets indicating time averaging.

Table 1. Parameters of the main three-dimensional simulations performed in the present study. For each simulation, we report Rayleigh number $\textit {Ra}$, grid resolution $N_{x}\times N_{z}\times N_{y}$, Nusselt number $\textit {Nu}$, and number of samples $N_{s}$ over which $\textit {Nu}$ is averaged (see appendix B for further details). Additional three-dimensional simulations (not listed here) have been performed at $\textit {Ra}=2.5\times 10^3, 5\times 10^3$.

In figure 2 we show the computed values of $\textit {Nu}$ as a function of $\textit {Ra}$ (filled circles), compared with data available in the literature (open triangles, Hewitt et al. Reference Hewitt, Neufeld and Lister2014). Also shown in light grey is the high-$\textit {Ra}$ portion of the $(\textit {Ra},\textit {Nu})$ parameter space covered in previous three-dimensional numerical simulations. The dash-dot line in figure 2 corresponds to the prediction of Hewitt et al. (Reference Hewitt, Neufeld and Lister2014) based on fitting of lower-$\textit {Ra}$ data ($\textit {Ra} \leqslant 2\times 10^4$), and amounting to a shifted linear variation, $\textit {Nu}=0.0096\textit {Ra} + 4.6$. As seen from figure 2, this fit overestimates our data by about $9\,\%$ at the highest $\textit {Ra}$; hence we believe that a different expression may be appropriate to more accurately characterize the trend towards the ultimate linear scaling.

Figure 2. Nusselt number as a function of the Rayleigh–Darcy number. Results obtained from the present numerical simulations are shown with filled circles, and the fitting curve, $\textit {Nu}=0.0081\textit {Ra}+0.067 \textit {Ra}^{0.61}$, is shown as a black solid line. Results obtained in previous simulations by Hewitt et al. (Reference Hewitt, Neufeld and Lister2014) (open triangles) and the extrapolated shifted linear scaling (blue dash-dotted line) are shown for comparison. The high-$\textit {Ra}$ portion of the $(\textit {Ra},\textit {Nu})$ parameter space covered by previous investigations is rendered by the light grey rectangle.

In figure 3 we thus show the compensated Nusselt number $\textit {Nu}/\textit {Ra}$, for a collection of current and previous numerical data obtained in three-dimensional (Hewitt et al. Reference Hewitt, Neufeld and Lister2014) and two-dimensional (Hewitt et al. Reference Hewitt, Neufeld and Lister2012; Wen et al. Reference Wen, Corson and Chini2015; De Paoli et al. Reference De Paoli, Zonta and Soldati2016) domains. According to the prediction of Hewitt et al. (Reference Hewitt, Neufeld and Lister2014), the asymptotic value $\textit {Nu}/\textit {Ra}=0.0096$ should be reached for $\textit {Ra} \gtrsim 10^5$. Our data (filled circles) for $\textit {Ra}>2\times 10^4$ fall well below that prediction, showing no clear evidence for inception of an ultimate regime. An improved prediction for the pre-asymptotic regime can be constructed by fitting the numerical data and relying on the compelling evidence that the asymptotic scaling should be linear (Hewitt Reference Hewitt2020). Data fitting then results in a functional representation of the type

(3.2)\begin{equation} \textit{Nu}/\textit{Ra}=0.0081+0.067\textit{Ra}^{-0.39} , \end{equation}

which amounts to linear dependence of $\textit {Nu}$ over $\textit {Ra}$ augmented with a sublinear corrective term, shown as a black solid line in figure 3. Compared with previous estimates, this correlation yields good prediction of $\textit {Nu}$ over a much larger range of $\textit {Ra}$, starting at $\textit {Ra} \sim 2.5\times 10^3$, with error no larger than $1.5\,\%$. Extrapolation to higher $\textit {Ra}$ than reached in the present computational campaign would yield the asymptotic scaling

(3.3)\begin{equation} \textit{Nu} = 0.0081 \textit{Ra} , \end{equation}

in the ultimate regime. Based on our data, we expect that to have deviations of less than $5\,\%$ from the ultimate regime, $\textit {Ra}$ in excess of $5\times 10^5$ is needed. At those huge Rayleigh numbers, which closely represent the case of geological $\textrm {CO}_{2}$ sequestration in deep saline aquifers, differences between the current prediction and previous predictions of Hewitt et al. (Reference Hewitt, Neufeld and Lister2014) would be as large as ${\sim }16\,\%$. These findings are contrasted in figure 3 with the case of two-dimensional flow, for which the data generally support $\textit {Nu} \approx 0.0069 \textit {Ra} + 2.75$, as suggested by Hewitt et al. (Reference Hewitt, Neufeld and Lister2012), and attainment of a linear trend at $\textit {Ra} \sim 3\times 10^4$. Hence, extrapolation of our results yields the prediction that in the ultimate regime the Nusselt number should be larger by about $18\,\%$ in the three-dimensional than in the two-dimensional case.

Figure 3. Compensated Nusselt number as a function of the Rayleigh–Darcy number, $\textit {Ra}$. Results obtained from the present numerical simulations, which we have run in three-dimensional and two-dimensional domains, are shown by filled circles ($\bullet$) and diamonds ($\blacklozenge$), respectively. For comparison purposes, a collection of previous data obtained in both two-dimensional domains (Hewitt et al. Reference Hewitt, Neufeld and Lister2012; De Paoli et al. Reference De Paoli, Zonta and Soldati2016; Wen et al. Reference Wen, Corson and Chini2015) ($\square$, $\triangledown$ and $\triangleright$, respectively) and three-dimensional domains (Hewitt et al. Reference Hewitt, Neufeld and Lister2014) ($\triangle$) is also shown, with open symbols. The black solid line denotes the best fit of our three-dimensional data, as from (3.2), and the blue solid line the shifted linear scaling $\textit {Nu}/\textit {Ra}=0.096+4.6/\textit {Ra}$ predicted by Hewitt et al. (Reference Hewitt, Neufeld and Lister2014). The scaling law $\textit {Nu}/\textit {Ra}=0.0069+2.75/\textit {Ra}$ proposed by Hewitt et al. (Reference Hewitt, Neufeld and Lister2012) for the two-dimensional case is shown with a solid red line. The asymptotic predictions for each scaling law are reported as dashed lines.

3.2. Flow structure

It is reasonable to expect that the heat transfer behaviour just discussed is reflected in the behaviour of the flow structure. We characterize the flow structure focusing on its signature in the near-boundary region where the contribution of both diffusion-driven and buoyancy-driven mechanisms is important. An estimate of the location at which diffusion and buoyancy are both important is given by the temperature boundary layer thickness, $\delta$, which, in Rayleigh–Darcy convection, is $\delta \sim 30/\textit {Ra}$ (see Otero et al. (Reference Otero, Dontcheva, Johnston, Worthing, Kurganov, Petrova and Doering2004), Hewitt et al. (Reference Hewitt, Neufeld and Lister2012), and also appendix A of the present manuscript). Accordingly, in figure 4 we show the temperature distributions in a horizontal plane at $y=50/\textit {Ra}$ from the bottom boundary. Indication of the domain size in diffusive–convective units, $x^*_d$, for the simulations at lower $\textit {Ra}$ is given by the white boxes in panel (a). The complex flow pattern near the boundary, which appears to be still developing at $\textit {Ra}<2\times 10^4$, seems to attain a roughly self-similar structure at larger $\textit {Ra}$. Hence, in the following we will mostly refer to the $S_8$ simulation, since the flow features we wish to discuss are emphasized at that value of $\textit {Ra}$.

Figure 4. Temperature distribution in a plane close to the bottom boundary $(x,y=50/\textit {Ra},z)$, for simulations $S_8$ (a), $S_4$ (b), $S_2$ (c) and $S_1$ (d). The domain size is also explicitly indicated in diffusive–convective units in each panel. Indication of the domain size – rescaled using the diffusive–convective length scale $x^*_d$ – for the simulations at lower $\textit {Ra}$ is given in panel (a) by the white boxes. The size of the superstructures, estimated as $\lambda _i=2{\rm \pi} \textit {Ra}/k_i$ with $k_i/\textit {Ra}$ defined as in figure 5, is also shown in each panel (black bars).

In figure 4 (and in particular in panel (a)) we clearly observe the presence of thin bright filaments which appear to be connected in a quasi-regular pattern of cells (Fu et al. Reference Fu, Cueto-Felgueroso and Juanes2013; Amooie, Soltanian & Moortgat Reference Amooie, Soltanian and Moortgat2018). These filaments correspond to regions where high-temperature fluid protrudes from the boundary, and are produced by tiny instabilities of the thermal boundary layer. Those are the three-dimensional counterpart of the protoplumes observed in previous two-dimensional studies (Hewitt et al. Reference Hewitt, Neufeld and Lister2012; Wen et al. Reference Wen, Corson and Chini2015; De Paoli et al. Reference De Paoli, Zonta and Soldati2016), and define the perimeter of rather homogeneous (dark) regions produced by downwelling motions (downwashes) of cold fluid towards the bottom boundary (see also the vertical slices shown in appendix A and comments therein). After impinging on the boundary, these downwashes spread horizontally and interact with the rising protoplumes, determining a dynamic pattern, in which some of the protoplumes eventually cluster into specific regions – thicker bright ridges (see figure 4a) – that define the boundaries of seemingly regular superstructures. These ‘supercells’, which encapsulate larger ensembles of smaller cells, have shape and dynamics recalling those of classical Rayleigh–Bénard turbulence (Stevens et al. Reference Stevens, Blass, Zhu, Verzicco and Lohse2018; Krug, Lohse & Stevens Reference Krug, Lohse and Stevens2020) or those of solar granules in the Sun's photosphere (Bahng & Schwarzschild Reference Bahng and Schwarzschild1961; Miesch Reference Miesch2005; Rieutord & Rincon Reference Rieutord and Rincon2010; Witze Reference Witze2020) (see also animations in the supplementary movies available at https://doi.org/10.1017/jfm.2020.1178 for a time-resolved rendering of the supercells). The size of the supercells $\lambda _i$, as inferred from the analysis of the wavenumber spectra (see figure 5 and the corresponding discussion), can be appreciated in figure 4. Since at large $\textit {Ra}$ the size of the supercells, which is approximately one-tenth of the box side length for simulation $S_8$ – see figure 4(a) – does not seem to evolve once the simulation has reached the steady state, we are keen to attribute to them an important role in the overall heat transfer mechanisms. At the same time, an important role is played by the smaller cells enclosed by supercells. Previous qualitative observations show a complex network of structures composed of small cells hierarchically nesting into supercells. To determine the size of such supercells, we use the time-averaged wavenumber power spectrum of the temperature signal, $P(k)$, which we obtain by taking the two-dimensional Fourier transform of the temperature field in a horizontal plane located near the boundary (such as the one sketched in figure 4a).

Figure 5. Power spectra $k_{r}P(k_{r})$ (solid lines and symbols) of temperature distribution near the boundary (sampled at the location $y=50/\textit {Ra}$, as in figure 4), shown as a function of the radial wavenumber $k_{r}=\sqrt {k_{x}^{2}+k_{z}^{2}}$, for all $\textit {Ra}$ considered here. Results of the $S_3$ simulation are omitted for ease of reading. In the inset of the figure we plot the power spectra of the temperature distribution at the centre (midplane) of the domain. Labels $k_i/\textit {Ra}$ and $k'_i/\textit {Ra}$ indicate the peaks of the spectra. The corresponding wavenumber is explicitly reported in the main panel.

In figure 5, we show $k_rP(k_r)$ as a function of the normalized radial wavenumber $k_r/\textit {Ra}$, where $k_r=\sqrt {k_x^2+k_z^2}$, and $k_x$ and $k_z$ are the wavenumbers along $x$ and $z$, respectively. Results of the $S_3$ simulation are not shown, for better readability. For all considered values of $\textit {Ra}$, a sharp peak occurs at low wavenumbers, which clearly characterizes the size of the supercells. The wavenumber at which this peak occurs, labelled with $k_i/\textit {Ra}$ in figure 5 (where $i$ refers to the simulation number as in table 1, i.e. $i=1, 2, 4, 8$), decreases for increasing $\textit {Ra}$, thus indicating growth of the supercells with $\textit {Ra}$. Measurements of the temperature power spectrum at the midplane of the domain (inset of figure 5) is also characterized by the presence of sharp peaks (labelled as $k'_i/\textit {Ra}$), which in this case are due to large megaplumes developing in the interior part of the domain (Hewitt et al. Reference Hewitt, Neufeld and Lister2012; Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2013; De Paoli et al. Reference De Paoli, Zonta and Soldati2016; Hewitt & Lister Reference Hewitt and Lister2017).

Direct comparison with the spectra near the boundary shows close coincidence of the spectral peaks, thus clearly indicating that supercells at the boundary are the footprint of the megaplumes that dominate the interior part of the flow. At larger $k_{r}$, the spectrum near the boundary features a local minimum, followed by rapid increase. Such increase, which is not present in the spectrum at the domain centre, is naturally associated with the small-scale protoplumes populating the region near the boundary (see the small filaments in figure 4).

4. Conclusions

We have reported results of numerical simulation of Rayleigh–Darcy convection inside a cubic porous layer up to Rayleigh–Darcy number $\textit {Ra}=8\times 10^4$. Exploiting the numerical dataset we show that the Nusselt number variation at sufficiently high $\textit {Ra}$ can be approximately expressed as $\textit {Nu}=0.0081\textit {Ra}+0.067\textit {Ra}^{0.61}$, with error bar no larger than $1.5\,\%$. Extrapolation of this prediction to the ultimate regime, which we expect to set in at $\textit {Ra} \gtrsim 5\times 10^5$, yields the asymptotic scaling $\textit {Nu}=0.0081\textit {Ra}$ – almost $16\,\%$ less than inferred in previous studies, and exceeding by about $18\,\%$ the Nusselt number found in two-dimensional simulations. We have further characterized the complex flow structure near the boundaries, which consists of small flow cells hierarchically nesting into supercells, and we have shown evidence that the supercells at the boundary are the footprint of megaplumes dominating the interior part of the flow. We believe that it would be possible to conduct direct investigations of the ultimate convective regime to check our extrapolations, by running a simulation at $\textit {Ra} \simeq 5\times 10^5$. However, we estimate a cost for this simulation of approximately 2 billion CPU hours, well beyond present-day capabilities.

Supplementary movies

Supplementary movies are available at https://doi.org/10.1017/jfm.2020.1178.

Acknowledgements

We are grateful to the CINECA supercomputing centre (under the grants ISCRA IsB20–3DSIMCON and PRACE Pra19-4938) and the Vienna Scientific Cluster (VSC) for generous allowance of computational resources. We also acknowledge Dr D. Hewitt and Dr B. Wen for providing data from Hewitt et al. (Reference Hewitt, Neufeld and Lister2012, Reference Hewitt, Neufeld and Lister2014) and Wen et al. (Reference Wen, Corson and Chini2015), respectively.

Funding

S.P. and F.Z. gratefully acknowledge financial support from the ERASMUS+ project 29415-EPP–1–2014–1–IT–EPPKA3–ECHE between TU Wien and La Sapienza University. The authors acknowledge financial support from the TU Wien University Library through its Open Access Funding Program.

Declaration of interests

The authors report no conflict of interest.

Appendix A. Structure of the flow in a vertical plane

In figure 6 we show the temperature distribution in a vertical $(z,y)$ plane located at $x=1/2$. A close-up view of the temperature field near the boundary is also offered to exhibit the thickness of the thermal boundary layer, which, scaling as $\delta \sim 1/\textit {Ra}$, would otherwise be too small to be observed. Consistently with previous two- and three-dimensional studies (Hewitt et al. Reference Hewitt, Neufeld and Lister2012, Reference Hewitt, Neufeld and Lister2014; Wen et al. Reference Wen, Corson and Chini2015; De Paoli et al. Reference De Paoli, Zonta and Soldati2016), we observe that small fingers of light fluid emerge from the bottom boundary and move upwards, and correspondingly small fingers of heavy fluid emerge from the top boundary and move downwards. The small fingers then merge to form megaplumes (large columnar structures that dominate the core region of the flow) which, under the vigorous action of buoyancy, increase their vertical velocity and reach the opposite boundary. Upon impact with the boundary, the megaplumes are deflected and create the complex flow patterns analysed in detail in the main body of the present manuscript. As expected, the strength and persistence of the megaplumes increase with increasing $\textit {Ra}$.

Figure 6. Temperature distribution from the simulations $S_1$ (a,b) and $S_8$ (c,d) in a vertical slice located at $x=1/2$. The dimensionless domain size in diffusive–convective units is indicated. A close-up view of the near-wall region, indicated with black squares in panels (a,c), is shown in panels (b,d).

Appendix B. Evaluation of the Nusselt number

To estimate the averaged Nusselt number, we first run each simulation (starting from a still fluid, plus random perturbations) until the initial transient is finished. From that time, labelled as $t_{ss}$, the Nusselt number oscillates around a statistically steady value. In figure 7, we show the time evolution of the Nusselt number, $\textit {Nu}(t)$, over a specific time window after $t_{ss}$, for the main simulations herein performed. $\textit {Nu}(t)$ is observed to oscillate within $1\,\%$ to $3\,\%$ of the mean value (indicated by the dashed line). In particular, larger values of $\textit {Ra}$ correspond to lower relative amplitude of the fluctuations. The time-averaged Nusselt number (reported in table 1) is computed over a time window corresponding to 128, 64, 16, 8 and 4 dimensionless time units for $S_1$, $S_2$, $S_3$, $S_4$ and $S_8$, respectively. Naturally, the duration of the time window, $\Delta t=t-t_{ss}$, over which statistics are evaluated is less at the largest $\textit {Ra}$. However, it is important to point out that the Nusselt number is not only averaged in time, but also in space (along $x$ and $z$) over the entire wall surface, which increases the size of the ensemble over which statistics are evaluated. The effective number of samples, $N_{s}$, used to estimate the Nusselt number (and also reported in table 1) is here defined as $N_{s}=N_x\times N_z\times N_T$, where $N_T$ is the time interval (expressed in dimensionless convective units, i.e. turnover times) over which statistics are performed (starting from $t_{ss}$). Once the averaged value of the Nusselt number is estimated, we determine the best-fit scaling as in (3.2), yielding prediction error of no more than $1.5\,\%$, for $\textit {Ra} \geqslant 2.5 \times 10^3$.

Figure 7. Time evolution of instantaneous horizontally-averaged Nusselt number. $\textit {Nu}$ is shown for simulations $S_1$ (a), $S_2$ (b), $S_3$ (c), $S_4$ (d) and $S_8$ (e). The behaviour of $\textit {Nu}(t)$ is reported for a portion of the simulation, after the steady state is achieved (the beginning of the steady-state regime is indicated by $t=t_{ss}$). We observe that $\textit {Nu}(t)$ (solid line) oscillates within $1\,\%$ to $3\,\%$ of the mean value (dashed line). For all simulations, $\textit {Nu}(t)$ is shown in the range corresponding to ${\pm }10\,\%$ of the time-averaged value. Since each point of the $\textit {Nu}(t)$ curve is the result of a space average (over the entire top and bottom wall), the larger $\textit {Ra}$ is, the larger the number of points over which the space average is evaluated, and so the lower the amplitude of the fluctuations.

References

REFERENCES

Amooie, M.A., Soltanian, M.R. & Moortgat, J. 2018 Solutal convection in porous media: comparison between boundary conditions of constant concentration and constant flux. Phys. Rev. E 98, 033118.CrossRefGoogle Scholar
Bahng, J. & Schwarzschild, M. 1961 Lifetime of solar granules. Astrophys. J. 134, 312.CrossRefGoogle Scholar
De Paoli, M. 2021 Influence of reservoir properties on the dynamics of a migrating current of carbon dioxide. Phys. Fluids 33 (1), 016602.CrossRefGoogle Scholar
De Paoli, M., Zonta, F. & Soldati, A. 2016 Influence of anisotropic permeability on convection in porous media: implications for geological $\textrm {CO}_{2}$ sequestration. Phys. Fluids 28 (5), 056601.CrossRefGoogle Scholar
De Paoli, M., Zonta, F. & Soldati, A. 2017 Dissolution in anisotropic porous media: modelling convection regimes from onset to shutdown. Phys. Fluids 29 (2), 026601.CrossRefGoogle Scholar
Emami-Meybodi, H., Hassanzadeh, H., Green, C.P. & Ennis-King, J. 2015 Convective dissolution of $\textrm {CO}_2$ in saline aquifers: progress in modeling and experiments. Intl J. Greenh. Gas Control 40, 238266.CrossRefGoogle Scholar
Fu, X., Cueto-Felgueroso, L. & Juanes, R. 2013 Pattern formation and coarsening dynamics in three-dimensional convective mixing in porous media. Phil. Trans. R. Soc. Lond. A 371 (2004), 20120355.Google ScholarPubMed
Graham, M.D. & Steen, P.H. 1994 Plume formation and resonant bifurcations in porous-media convection. J. Fluid Mech. 272, 6790.CrossRefGoogle Scholar
Hassanzadeh, P., Chini, G.P. & Doering, C.R. 2014 Wall to wall optimal transport. J. Fluid Mech. 751, 627662.CrossRefGoogle Scholar
He, X., Funfschilling, D., Nobach, H., Bodenschatz, E. & Ahlers, G. 2012 Transition to the ultimate state of turbulent Rayleigh–Bénard convection. Phys. Rev. Lett. 108 (2), 024502.CrossRefGoogle ScholarPubMed
Hewitt, D.R. 2020 Vigorous convection in porous media. Proc. R. Soc. Lond. A 476 (2239), 20200111.Google ScholarPubMed
Hewitt, D.R. & Lister, J.R. 2017 Stability of three-dimensional columnar convection in a porous medium. J. Fluid Mech. 829, 89111.CrossRefGoogle Scholar
Hewitt, D.R., Neufeld, J.A. & Lister, J.R. 2012 Ultimate regime of high Rayleigh number convection in a porous medium. Phys. Rev. Lett. 108 (22), 224503.CrossRefGoogle Scholar
Hewitt, D.R., Neufeld, J.A. & Lister, J.R. 2013 Stability of columnar convection in a porous medium. J. Fluid Mech. 737, 205231.CrossRefGoogle Scholar
Hewitt, D.R., Neufeld, J.A. & Lister, J.R. 2014 High Rayleigh number convection in a three-dimensional porous medium. J. Fluid Mech. 748, 879895.CrossRefGoogle Scholar
Hidalgo, J.J., Fe, J., Cueto-Felgueroso, L. & Juanes, R. 2012 Scaling of convective mixing in porous media. Phys. Rev. Lett. 109 (26), 264503.CrossRefGoogle ScholarPubMed
Horton, C.W. & Rogers, F.T. 1945 Convection currents in a porous medium. J. Appl. Phys. 16 (6), 367370.CrossRefGoogle Scholar
Howard, L. 1964 Convection at high Rayleigh number. In Proceedings of the 11th International Congress of Applied Mechanics, pp. 1109–1115. Springer.CrossRefGoogle Scholar
Huppert, H.E. & Neufeld, J.A. 2014 The fluid mechanics of carbon dioxide sequestration. Annu. Rev. Fluid Mech. 46, 255272.CrossRefGoogle Scholar
Kim, J. & Moin, P. 1985 Application of a fractional-step method to incompressible Navier–Stokes equations. J. Comput. Phys. 59, 308323.CrossRefGoogle Scholar
Kraichnan, R.H. 1962 Turbulent thermal convection at arbitrary Prandtl number. Phys. Fluids 5 (11), 13741389.CrossRefGoogle Scholar
Krug, D., Lohse, D. & Stevens, R.J.A.M. 2020 Coherence of temperature and velocity superstructures in turbulent Rayleigh–Bénard flow. J. Fluid Mech. 887, A2.CrossRefGoogle Scholar
Landman, A.J. & Schotting, R.J. 2007 Heat and brine transport in porous media: the Oberbeck–Boussinesq approximation revisited. Trans. Porous Med. 70 (3), 355373.CrossRefGoogle Scholar
Lapwood, E.R. 1948 Convection of a fluid in a porous medium. Math. Proc. Cambridge 44, 508521.CrossRefGoogle Scholar
Lepot, S., Aumaître, S. & Gallet, B. 2018 Radiative heating achieves the ultimate regime of thermal convection. Proc. Natl Acad. Sci. USA 115 (36), 89378941.CrossRefGoogle ScholarPubMed
Malkus, W. 1954 Discrete transitions in turbulent convection. Proc. R. Soc. Lond. A 225, 185195.Google Scholar
Miesch, M.S. 2005 Large-scale dynamics of the convection zone and tachocline. Living Rev. Sol. Phys. 2 (1), 1.CrossRefGoogle Scholar
Orlandi, P. 2000 Fluid Flow Phenomena: A Numerical Toolkit. Kluwer.CrossRefGoogle Scholar
Otero, J., Dontcheva, L.A., Johnston, H., Worthing, R.A., Kurganov, A., Petrova, G. & Doering, C.R. 2004 High-Rayleigh number convection in a fluid-saturated porous layer. J. Fluid Mech. 500, 263281.CrossRefGoogle Scholar
Pirozzoli, S. 2014 Revisiting the mixing-length hypothesis in the outer part of turbulent wall layers: mean flow and wall friction. J. Fluid Mech. 745, 378397.CrossRefGoogle Scholar
Pirozzoli, S., Bernardini, M., Verzicco, R. & Orlandi, P. 2017 Mixed convection in turbulent channels with unstable stratification. J. Fluid Mech. 821, 482516.CrossRefGoogle Scholar
Riaz, A. & Cinar, Y. 2014 Carbon dioxide sequestration in saline formations: part I – review of the modeling of solubility trapping. J. Petrol. Sci. Engng 124, 367380.CrossRefGoogle Scholar
Rieutord, M. & Rincon, F. 2010 The Sun's supergranulation. Living Rev. Sol. Phys. 7 (1), 2.CrossRefGoogle Scholar
Slim, A.C. 2014 Solutal-convection regimes in a two-dimensional porous medium. J. Fluid Mech. 741, 461491.CrossRefGoogle Scholar
Stevens, R.J.A.M., Blass, A., Zhu, X., Verzicco, R. & Lohse, D. 2018 Turbulent thermal superstructures in Rayleigh–Bénard convection. Phys. Rev. Fluids 3 (4), 041501.CrossRefGoogle Scholar
Wang, B.F., Zhou, Q. & Sun, C. 2020 Vibration-induced boundary-layer destabilization achieves massive heat-transport enhancement. Sci. Adv. 6 (21), eaaz8239.CrossRefGoogle ScholarPubMed
Wen, B., Akhbari, D., Zhang, L. & Hesse, M.A. 2018 Convective carbon dioxide dissolution in a closed porous medium at low pressure. J. Fluid Mech. 854, 5687.CrossRefGoogle Scholar
Wen, B., Corson, L.T. & Chini, G.P. 2015 Structure and stability of steady porous medium convection at large Rayleigh number. J. Fluid Mech. 772, 197224.CrossRefGoogle Scholar
Wen, B., Dianati, N., Lunasin, E., Chini, G.P. & Doering, C.R. 2012 New upper bounds and reduced dynamical modeling for Rayleigh–Bénard convection in a fluid saturated porous layer. Commun. Nonlinear Sci. Numer. Simul. 17 (5), 21912199.CrossRefGoogle Scholar
Witze, A. 2020 World's most powerful solar telescope is up and running. https://www.nature.com/articles/d41586-020-00224-z (accessed 19 July 2020).CrossRefGoogle Scholar
Zhu, X., Mathai, V., Stevens, R.J.A.M., Verzicco, R. & Lohse, D. 2018 Transition to the ultimate regime in two-dimensional Rayleigh–Bénard convection. Phys. Rev. Lett. 120 (14), 144502.CrossRefGoogle ScholarPubMed
Zonta, F. & Soldati, A. 2018 Stably stratified wall-bounded turbulence. Appl. Mech. Rev. 70, 040801.CrossRefGoogle Scholar
Figure 0

Figure 1. Sketch of the computational domain – a cube with side length $h^*$ – used to study Rayleigh–Darcy convection. The flow is heated at the bottom, $\theta ^*(y^*=0)=\theta ^*_{{max}}$, and cooled at the top, $\theta ^*(y^*=h^*)=\theta ^*_{{min}}$, and boundaries in the $x^*$ and $z^*$ directions are assumed to be periodic. The gravity acceleration ($\boldsymbol {g}$) points downwards. The temperature distribution $\theta ^*$ for the case $\textit {Ra}=1\times 10^4$ is also shown for illustrative purposes on the side boundaries and in a plane very close to the top boundary (i.e. at a distance of $50h^*/\textit {Ra}$ from the top boundary).

Figure 1

Table 1. Parameters of the main three-dimensional simulations performed in the present study. For each simulation, we report Rayleigh number $\textit {Ra}$, grid resolution $N_{x}\times N_{z}\times N_{y}$, Nusselt number $\textit {Nu}$, and number of samples $N_{s}$ over which $\textit {Nu}$ is averaged (see appendix B for further details). Additional three-dimensional simulations (not listed here) have been performed at $\textit {Ra}=2.5\times 10^3, 5\times 10^3$.

Figure 2

Figure 2. Nusselt number as a function of the Rayleigh–Darcy number. Results obtained from the present numerical simulations are shown with filled circles, and the fitting curve, $\textit {Nu}=0.0081\textit {Ra}+0.067 \textit {Ra}^{0.61}$, is shown as a black solid line. Results obtained in previous simulations by Hewitt et al. (2014) (open triangles) and the extrapolated shifted linear scaling (blue dash-dotted line) are shown for comparison. The high-$\textit {Ra}$ portion of the $(\textit {Ra},\textit {Nu})$ parameter space covered by previous investigations is rendered by the light grey rectangle.

Figure 3

Figure 3. Compensated Nusselt number as a function of the Rayleigh–Darcy number, $\textit {Ra}$. Results obtained from the present numerical simulations, which we have run in three-dimensional and two-dimensional domains, are shown by filled circles ($\bullet$) and diamonds ($\blacklozenge$), respectively. For comparison purposes, a collection of previous data obtained in both two-dimensional domains (Hewitt et al.2012; De Paoli et al.2016; Wen et al.2015) ($\square$, $\triangledown$ and $\triangleright$, respectively) and three-dimensional domains (Hewitt et al.2014) ($\triangle$) is also shown, with open symbols. The black solid line denotes the best fit of our three-dimensional data, as from (3.2), and the blue solid line the shifted linear scaling $\textit {Nu}/\textit {Ra}=0.096+4.6/\textit {Ra}$ predicted by Hewitt et al. (2014). The scaling law $\textit {Nu}/\textit {Ra}=0.0069+2.75/\textit {Ra}$ proposed by Hewitt et al. (2012) for the two-dimensional case is shown with a solid red line. The asymptotic predictions for each scaling law are reported as dashed lines.

Figure 4

Figure 4. Temperature distribution in a plane close to the bottom boundary $(x,y=50/\textit {Ra},z)$, for simulations $S_8$ (a), $S_4$ (b), $S_2$ (c) and $S_1$ (d). The domain size is also explicitly indicated in diffusive–convective units in each panel. Indication of the domain size – rescaled using the diffusive–convective length scale $x^*_d$ – for the simulations at lower $\textit {Ra}$ is given in panel (a) by the white boxes. The size of the superstructures, estimated as $\lambda _i=2{\rm \pi} \textit {Ra}/k_i$ with $k_i/\textit {Ra}$ defined as in figure 5, is also shown in each panel (black bars).

Figure 5

Figure 5. Power spectra $k_{r}P(k_{r})$ (solid lines and symbols) of temperature distribution near the boundary (sampled at the location $y=50/\textit {Ra}$, as in figure 4), shown as a function of the radial wavenumber $k_{r}=\sqrt {k_{x}^{2}+k_{z}^{2}}$, for all $\textit {Ra}$ considered here. Results of the $S_3$ simulation are omitted for ease of reading. In the inset of the figure we plot the power spectra of the temperature distribution at the centre (midplane) of the domain. Labels $k_i/\textit {Ra}$ and $k'_i/\textit {Ra}$ indicate the peaks of the spectra. The corresponding wavenumber is explicitly reported in the main panel.

Figure 6

Figure 6. Temperature distribution from the simulations $S_1$ (a,b) and $S_8$ (c,d) in a vertical slice located at $x=1/2$. The dimensionless domain size in diffusive–convective units is indicated. A close-up view of the near-wall region, indicated with black squares in panels (a,c), is shown in panels (b,d).

Figure 7

Figure 7. Time evolution of instantaneous horizontally-averaged Nusselt number. $\textit {Nu}$ is shown for simulations $S_1$ (a), $S_2$ (b), $S_3$ (c), $S_4$ (d) and $S_8$ (e). The behaviour of $\textit {Nu}(t)$ is reported for a portion of the simulation, after the steady state is achieved (the beginning of the steady-state regime is indicated by $t=t_{ss}$). We observe that $\textit {Nu}(t)$ (solid line) oscillates within $1\,\%$ to $3\,\%$ of the mean value (dashed line). For all simulations, $\textit {Nu}(t)$ is shown in the range corresponding to ${\pm }10\,\%$ of the time-averaged value. Since each point of the $\textit {Nu}(t)$ curve is the result of a space average (over the entire top and bottom wall), the larger $\textit {Ra}$ is, the larger the number of points over which the space average is evaluated, and so the lower the amplitude of the fluctuations.

Pirozzoli et al. supplementary movie 1

Temperature distribution at the distance 50/Ra from the bottom wall (S1, Ra=10000).

Download Pirozzoli et al. supplementary movie 1(Video)
Video 9.4 MB

Pirozzoli et al. supplementary movie 2

Temperature distribution at the distance 50/Ra from the bottom wall (S8, Ra=80000).

Download Pirozzoli et al. supplementary movie 2(Video)
Video 50 MB