2. Boussinesq equations and numerical methods

Let us consider a fluid layer between two horizontal plates heated from below and cooled from above. We employ the Oberbeck–Boussinesq approximation, wherein the density variations are only significant in the buoyancy term. The time evolutions of velocity field ${\boldsymbol {u}}({\boldsymbol {x}},t)=u{\boldsymbol {e}}_{x}+v{\boldsymbol {e}}_{y}+w{\boldsymbol {e}}_{z}$ and temperature field $T({\boldsymbol {x}},t)$ are described by the Boussinesq equations

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

(2.2)\begin{gather}\frac{\partial {\boldsymbol{u}}}{\partial t}+({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}){\boldsymbol{u}}=-\frac{1}{\rho}\boldsymbol{\nabla} p+\nu\boldsymbol{\nabla}^{2}{\boldsymbol{u}}+g\alpha T{\boldsymbol{e}}_{z}, \end{gather}

(2.3)\begin{gather}\frac{\partial T}{\partial t}+({\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla})T=\kappa\boldsymbol{\nabla}^{2}T, \end{gather}

where $p({\boldsymbol {x}},t)$ is pressure and $\rho , \nu , g, \alpha$ and $\kappa$ are the mass density, kinematic viscosity, acceleration due to gravity, volumetric thermal expansivity and thermal diffusivity, respectively. Vectors ${\boldsymbol {e}}_{x}$ and ${\boldsymbol {e}}_{y}$ are mutually orthogonal unit vectors in the horizontal directions, while ${\boldsymbol {e}}_{z}$ is a unit vector in the vertical direction. The two no-slip and impermeable horizontal plates are positioned at $z=0$ and $z=H$, and the top (or bottom) wall surface is held at a lower (or higher) constant temperature:

(2.4a–c)\begin{equation} \boldsymbol{u}(z=0)=\boldsymbol{u}(z=H)=\boldsymbol{0};\quad T(z=0)=\Delta T>0,\quad T(z=H)=0. \end{equation}

The velocity and temperature fields are supposed to be periodic in the $x$ and $y$ directions with the same period $L_{x}=L_{y}=L$. The thermal convection is characterised by the Rayleigh number $Ra$ and the Prandtl number $Pr$:

(2.5a,b)\begin{equation} Ra=\frac{g\alpha\Delta T H^{3}}{\nu\kappa},\quad Pr=\frac{\nu}{\kappa}. \end{equation}

The vertical heat flux is quantified by the Nusselt number:

(2.6)\begin{equation} Nu=\frac{-\kappa{\langle \partial T/\partial z \rangle}_{xyt}+{\langle wT \rangle}_{xyt}}{\kappa\Delta T/H}=1+\frac{H}{\kappa\Delta T}{\left\langle wT \right\rangle}_{xyzt}, \end{equation}

where $\langle \cdot \rangle _{xyt}$ and $\langle \cdot \rangle _{xyzt}$ represent the horizontal and time average and the volume and time average, respectively. The second equality is given by the volume and time average of (2.3).

Equations (2.1)–(2.3) are discretised employing a spectral Galerkin method based on the Fourier series expansion in the periodic horizontal directions and the Chebyshev polynomial expansion in the vertical direction. The nonlinear terms are evaluated using a spectral collocation method. Aliasing errors are removed with the aid of the $2/3$ (or $1/2$) rule for the Fourier (or Chebyshev) transform. Time advancement is performed with the Crank–Nicholson scheme and the second-order Adams–Bashforth scheme for the diffusion terms and the rest, respectively. The nonlinear steady solutions are obtained using the Newton–Krylov iteration (for more details, see § 3 and appendix A in Motoki, Kawahara & Shimizu (Reference Motoki, Kawahara and Shimizu2018b)).

In this paper, we present the steady solution and the turbulent states in the horizontally square periodic domain with $L/H={\rm \pi} /2\approx 1.57$ for $Pr=1$. The domain is the same as that of the optimal states derived by Motoki et al. (Reference Motoki, Kawahara and Shimizu2018a). The numerical process is conducted on $128^{3}$ grid points for $Ra<10^7$ and $256^{3}$ grid points for $Ra\ge 10^7$. In the smaller (or larger) domain with $L/H=1$ (or $L/H=2{\rm \pi} /3.117\approx 2.02$), we have confirmed that the effects of the domain size on the heat flux and the spatial structures are insignificant, which is described in the following sections. The domain size dependence is shown in appendix B with the Prandtl number dependence and the adequacy of the spatial resolution.

All the 3-D steady solutions presented in this paper satisfy the ${\rm \pi} /2$-rotation symmetry

(2.7)\begin{equation} [u,v,w,T](x,y,z)=[v,-u,w,T](\,y,-x,z) \end{equation}

as well as the reflection symmetry

(2.8)\begin{equation} [u,v,w,T](x,y,z)=[-u,v,w,T](-x,y,z) \nonumber\\ =[u,-v,w,T](x,-y,z) \end{equation}

and the shift-and-reflection symmetry

(2.9)\begin{equation} [u,v,w,T](x,y,z)=[u,v,-w,1-T](x+L/2,y+L/2,1-z). \end{equation}

These symmetries are not imposed on the solutions explicitly at $Ra\lesssim 10^6$. At $Ra\sim 10^{7}$, we impose the symmetries (2.7)–(2.9) in order to reduce the computational degrees of freedom.

3. Three-dimensional steady solution

3.1. The $Nu$–$Ra$ scaling

The 3-D and 2-D steady solutions in the domain with $L/H={\rm \pi} /2$ for $Pr=1$ bifurcate from the conduction state at $Ra\approx 1879$, as shown in figure 1. The 2-D roll solution is stable up to $Ra\approx 3.55\times 10^{4}$, and subsequently a time-periodic 3-D solution appears (the stable periodic solution has been obtained by the time evolution). Figure 2 presents $Nu$ in the 3-D and 2-D steady solutions as a function of $Ra$ at a higher $Ra$. The red line shows the 3-D steady solution, and the open and filled circles represent the present turbulent data in the horizontally square periodic domain and the experimental data in a cylindrical container (Niemela & Sreenivasan Reference Niemela and Sreenivasan2006), respectively. The blue circles represent the 2-D optimised steady solutions for $Pr=1$ obtained by Sondak et al. (Reference Sondak, Smith and Waleffe2015). Although the 3-D steady solution at high $Ra$ exhibits smaller $Nu$ than the 2-D optimised solutions, it maintains larger $Nu$ than the turbulent states even at $Ra\sim 10^{7}$. The inset shows $Nu$ compensated by $Ra^{\gamma }$. The scaling exponent $\gamma$ of turbulent states shows $2/7$ for $Ra\lesssim 10^{7}$; at higher $Ra$, it changes to $0.31$. Such a transition has been experimentally and numerically observed for $Pr\sim 1$ (Castaing et al. Reference Castaing, Gunaratne, Heslot, Kadanoff, Libchaber, Thomae, Wu, Zaleski and Zanetti1989; Silano, Sreenivasan & Verzicco Reference Silano, Sreenivasan and Verzicco2010). The exponent of the heat flux in the 3-D steady solution is greater than $2/7$ but less than $1/3$, with $0.31$ being the closest, whereas that in the 2-D steady solution with fixed $L/H={\rm \pi} /2$ (blue line) is smaller than $2/7$. The bumps in the 3-D steady solution correspond to the remarkable changes in the thermal and flow structures, which will be shown in § 3.3. At $Ra\sim 10^{4}$ the heat flux of the 3-D steady solution surpasses maxima in the 2-D steady solutions. As shown in appendix B, the 3-D steady solution with the smaller horizontal period ($L/H=1$) shows larger $Nu$ at $Ra\sim 10^{5}$. The results imply that a family of the 3-D steady solutions with optimised horizontal periods might achieve maximal heat transfer in the steady Rayleigh–Bénard convection. The orange dashed line indicates the optimal scaling $Nu-1=0.0236Ra^{1/2}$ (Motoki et al. Reference Motoki, Kawahara and Shimizu2018a) given by the 3-D optimal states maximising the wall-to-wall transport, and it is quite close to the rigorous upper bound, $Nu-1=0.02634Ra^{1/2}$ (Plasting & Kerswell Reference Plasting and Kerswell2003). Although the 3-D optimal states exhibiting significantly high heat flux require an external body force different from buoyancy, the optimal state can be continuously connected to the present 3-D steady solution of the full Boussinesq equations by a homotopy from the body force to the buoyancy, as shown in appendix A.

Figure 2. Nusselt number $Nu-1$ as a function of Rayleigh number $Ra$. The red and blue solid lines represent the 3-D and 2-D steady solutions, respectively. The open black circles represent the present turbulent data obtained in the horizontally square periodic domain, and the filled ones are the experimental turbulent data in a cylindrical container (Niemela & Sreenivasan Reference Niemela and Sreenivasan2006). The blue filled circles indicate 2-D optimised steady solutions for $Pr=1$ obtained by Sondak et al. (Reference Sondak, Smith and Waleffe2015). The orange solid and dashed lines indicate the upper bound $Nu-1=0.02634Ra^{1/2}$ (Plasting & Kerswell Reference Plasting and Kerswell2003) and the optimal scaling $Nu-1=0.0236Ra^{1/2}$, respectively, evaluated from the wall-to-wall optimal transport states (Motoki et al. Reference Motoki, Kawahara and Shimizu2018a). The inset shows $Nu$ compensated by $Ra^{\gamma }$: $\gamma =2/7$ (plot A); $\gamma =0.31$ (plot B); $\gamma =1/3$ (plot C).

3.2. Mean temperature and root-mean-square profiles

Figure 3 compares the mean temperature ${\langle T \rangle }_{xyt}$, root-mean-square (RMS) temperature $T_{rms}={\langle {(T-{\langle T \rangle }_{xyt})}^{2} \rangle }_{xyt}^{1/2}$ and RMS vertical velocity $w_{rms}={\langle w^{2} \rangle }_{xyt}^{1/2}$ for the 3-D steady solution with those for the turbulent states at $10^{5}\le Ra\le 10^{7}$. In figure 3(b,d,f) the distance to the wall, $z$, is normalised by $\delta$, where $\delta$ is the thermal conduction layer thickness given by $\delta /H=1/(2Nu)$. The mean temperature in the 3-D steady solution well follows that in the turbulent states; furthermore, the corresponding RMS values are also in good agreement with each other. In the bulk region, all mean temperature profiles are flattened, as a result of the nearly complete mixing by large-scale circulation. The temperature difference $\Delta T/2$ exists only at the thermal conduction layer, $0\le z\lesssim 2\delta /H=1/Nu$, and $T_{rms}$ exhibits a peak at $z/\delta \approx 1$. If the advection, diffusion and buoyancy terms in the Navier–Stokes equation (2.2) at the conduction layer are balanced as

(3.1)\begin{equation} \frac{w'^{2}}{\delta}\sim\nu\frac{w'}{\delta^{2}}\sim g\alpha\Delta T \end{equation}

(the balance between the advection and diffusion terms is given by that in the energy equation (2.3) for $\nu \sim \kappa$), the near-wall vertical velocity would then be

(3.2)\begin{equation} w'\sim{(g\alpha\Delta T\delta)}^{1/2}\sim Ra^{1/3}\frac{\kappa}{H}, \end{equation}

yielding the scaling law

(3.3)\begin{equation} Nu\sim Ra^{1/3} \end{equation}

(Kawano et al. Reference Kawano, Motoki, Shimizu and Kawahara2021), which has been given by Malkus’ theory (Malkus Reference Malkus1954) and the argument of Grossmann & Lohse (Reference Grossmann and Lohse2000). As shown in figure 3(f), the RMS vertical velocity $w_{rms}$ scales as $Ra^{1/3}\kappa /H$ near the wall $z/\delta \sim 1$.

Figure 3. (a,b) Mean temperature, (c,d) RMS temperature and (e,f) RMS vertical velocity as a function of (a,c,e) $z/H$ and (b,d,f) $z/\delta$ in the 3-D steady solution and the turbulent state. Here $\delta$ is the thermal conduction layer thickness that scales as $\delta /H=1/(2Nu)$.

In the bulk region, where the effects of viscosity or thermal conduction are insignificant, the advection and buoyancy terms are balanced as

(3.4)\begin{equation} \frac{U_{b}^{2}}{H}\sim g\alpha\Delta T', \end{equation}

and the dominance of convective heat transfer in bulk and the scaling (3.3) give us

(3.5)\begin{equation} \frac{U_{b}\Delta T'}{\kappa \Delta T/H}\sim Ra^{1/3} \end{equation}

(Kawano et al. Reference Kawano, Motoki, Shimizu and Kawahara2021). Equations (3.4) and (3.5) yield the bulk velocity scales as

(3.6)\begin{equation} U_{b}\sim Pr^{-2/3}Ra^{4/9}\nu/H. \end{equation}

The Reynolds number based on $U_{b}$ scales with $Pr^{-2/3}Ra^{4/9}$, exhibiting consistency with the scaling suggested by Grossmann & Lohse (Reference Grossmann and Lohse2000) when both energy and scalar dissipation are bulk-dominated. Figure 4 shows the Reynolds number $Re$ defined as

(3.7)\begin{equation} Re=\frac{{\langle w^{2} \rangle}_{xyzt}^{1/2}H}{\nu} \end{equation}

as a function of $Ra$. The Reynolds number $Re$ of the 3-D steady solution is comparable with that of the turbulent states exhibiting ${Re}\sim Ra^{4/9}$ at a high $Ra$. In contrast, the 2-D steady solution with fixed $L/H={\rm \pi} /2$ seems to approach $Re\sim Ra^{1/2}$, implying that $\langle w^{2}\rangle _{xyzt}^{1/2}$ scales with the buoyancy-induced terminal velocity $U={(g\alpha \Delta TH)}^{1/2}$.

Figure 4. Reynolds number $Re$ as a function of Rayleigh number $Ra$. The red and blue solid lines represent the 3-D and 2-D steady solutions, respectively. The open circles represent the present turbulent data obtained in the horizontally square periodic domain. The inset shows $Nu$ compensated by $Ra^{\gamma }$: $\gamma =4/9$ (plot A); $\gamma =1/2$ (plot B).

3.3. Thermal and flow structures

Figure 5 visualises the thermal and flow structures in the 3-D steady solution and the turbulent state. The yellow objects represent the isosurfaces of temperature $T/\Delta T=0.6$, representing high-temperature plumes; the grey objects represent the vortex structures visualised by the positive second invariant of the velocity gradient tensor

(3.8)\begin{equation} Q=-\frac{1}{2}\frac{\partial u_{i}}{\partial x_{j}}\frac{\partial u_{j}}{\partial x_{i}}. \end{equation}

At $Ra=10^{5}$ (figure 5a), the 3-D steady solution consists of large-scale thermal plumes as well as accompanying near-wall small-scale plume and vortex structures near the vertical planes $x/H=0$ and $y/H=0$. As $Ra$ for the 3-D steady solution increases, smaller thermal plume structures (and relevant smaller and stronger tube-like vortex structures) appear near the walls without affecting the already existing large-scale structures. As $Ra$ increases from $10^{5}$ to $10^{6}$, thermal and vortex structures newly emerge along the diagonals $y={\pm }x$ near the walls (figure 5b). As $Ra$ increases further from $10^{6}$ to $10^{7}$, smaller-scale structures appear along $x/H=\pm {\rm \pi}/8$ and $y/H={\pm }{\rm \pi} /8$ near the walls, and the large secondary plumes approach the opposite wall (figure 5c). The bumps of $Nu$ seen in figure 2 correspond to these remarkable changes in the thermal and vortex structures.

Figure 5. Thermal and flow structures in (a–c) the 3-D steady solution and (d) the turbulent state at (a) $Ra=10^5$, (b) $Ra=10^6$ and (c,d) $Ra=10^7$. The yellow and grey objects represent the isosurfaces of the temperature $T/\Delta T=0.6$ and the positive second invariant of the velocity gradient tensor, (a) $Q/(\kappa ^2/H^4)=1.28\times 10^5$, (b) $Q/(\kappa ^2/H^4)=1.28\times 10^6$ and (c,d) $Q/(\kappa ^2/H^4)=8\times 10^7$, respectively. The contours represent temperature $T$ in the plane $y/H={\rm \pi} /4 (=-{\rm \pi} /4)$, and the velocity vectors $(u,w)$ in the enlarged views in (c,d) are superposed.

At $Ra=10^7$ (figure 5c), we observe sheet-like thermal plumes with the smallest-scale vortices, which are quite similar to those observed in the snapshot of the turbulent state (figure 5d). The smallest-scale structures appear in the thermal conduction layer, and the sizes of the plumes and vortices scale with the thickness $\delta$. In 2-D optimised steady solutions (Sondak et al. Reference Sondak, Smith and Waleffe2015; Waleffe et al. Reference Waleffe, Boonkasame and Smith2015), the appearance of small-scale thermal and vortex structures has been observed with the shrinking horizontal period as $Ra$ increases, and the scaling $Nu\sim Ra^{0.31}$ is achieved by a family of solutions with smaller horizontal periods. It should be stressed that the single 3-D steady solution spontaneously reproduces the multi-scale coherent structures of convective turbulence.

4. Hierarchical vortices and energy transfer in wavenumber space

The 3-D steady solution is found to reproduce the near-wall coherent structures in convective turbulence; however, the flow in bulk should be addressed. The developed turbulence organises hierarchical coherent vortex structures of various scales (Goto, Saito & Kawahara Reference Goto, Saito and Kawahara2017; Motoori & Goto Reference Motoori and Goto2019). The smallest-scale vortex structures can still be extracted by employing the isosurface of $Q$, as shown in figure 5; however, it is difficult to identify large- and intermediate-scale structures. To examine the hierarchy of multi-scale vortices in the 3-D steady solution, we consider coarse-graining the velocity field ${\boldsymbol {u}}$. The coarse-grained velocity field ${\boldsymbol {u}}^{*}$ is obtained using the Gaussian low-pass filter (Lozano-Durán, Holzner & Jimenez Reference Lozano-Durán, Holzner and Jimenez2016; Motoori & Goto Reference Motoori and Goto2019) as follows:

(4.1)\begin{equation} {\boldsymbol{u}}^{*}({\boldsymbol{x}})=\int_{V} {a}{\boldsymbol{u}}({\boldsymbol{x}}')\exp{\left\{ -{\left( \frac{{\rm \pi}\Delta r}{\sigma} \right)}^2 \right\}}\,\textrm{d}\kern0.7pt{\boldsymbol{x}}', \end{equation}

where $\Delta r=|\boldsymbol {x}'-\boldsymbol {x}|$, $\sigma$ is the filter width and $a$ is a constant such that the integral of the kernel over the control volume $V$ is unity. In the wall-normal direction, the Gaussian filter is applied by reflecting it at the wall (Lozano-Durán et al. Reference Lozano-Durán, Holzner and Jimenez2016). Figure 6 shows hierarchical vortex structures in the 3-D steady solution at $Ra=2.6\times 10^7$. Non-filtered structures are shown in figure 6(a), and the isosurfaces of $Q$ of the filtered velocity ${{\boldsymbol {u}}}^{*}$ with $\sigma =H(=2L/{\rm \pi} ),L/2,L/4,L/8$ and $L/16$ are displayed in figure 6(b–f), respectively. The blue objects in figure 6(b) are the largest-scale structures corresponding to the large-scale circulation, whereas the red ones in figure 6(f) are the smallest-scale structures of size $\sigma /2=L/32\approx 2\delta$, which coincides with the size of vortices observed in the non-filtered field (figure 6a). The light blue, green and light red objects in figure 6(c,d,e) illustrate the intermediate-scale vortex structures of eight, four and two times the size of the smallest-scale vortices, respectively. The smaller-scale vortex structures exist closer to the wall, while the intermediate-scale ones are observed in the bulk region in figure 6(d,e).

Figure 6. Hierarchical vortex structures visualised by coarse-graining with Gaussian low-pass filter. (a) The yellow and red objects represent the isosurfaces of the non-filtered $T/\Delta T=0.6$ and $Q/(\kappa ^2/H^4)=2\times 10^8$, respectively. (b–h) The vortex structures are visualised by the isosurfaces of $Q/(\kappa ^2/H^4)$ of the filtered velocity field with filter widths of $\sigma =H(=2L/{\rm \pi} )$ (blue), $\sigma =L/4$ (light blue), $\sigma =L/8$ (green), $\sigma =L/16$ (light red) and $\sigma =L/32$ (red); they are superposed in (g,h). The isosurface levels are (blue) $5\times 10^5$, (light blue) $4\times 10^6$, (green) $1.2\times 10^7$, (light red) $3\times 10^7$ and (red) $1.6\times 10^8$.

Figure 6(g,h) shows the superposed structures, and from their spatial distribution, it is conjectured that the bulk flow is composed of multi-scale coherent structures. Figure 7(a) shows the energy spectrum $E(k,z)$ of the 3-D steady solution at the centre of the fluid layer, $z=H/2$, and the corresponding turbulence spectrum at $Ra=2.6\times 10^{7}$. The energy spectrum $E(k,z)$ is defined as

(4.2)\begin{equation} E(k,z)=\frac{L}{2{\rm \pi}}\sum_{k-{\Delta k}/{2}<|\boldsymbol{k}_{2D}|<k+{\Delta k}/{2}}\frac{1}{2}{\left\langle {|\tilde{\boldsymbol{u}}(\boldsymbol{k}_{2D},z)|}^{2} \right\rangle}_{t}, \end{equation}

where $\tilde {(\cdot )}$ indicates the Fourier coefficients in the periodic ($x$ and $y$) directions and ${\langle \cdot \rangle }_{t}$ represents the time average. Here $\boldsymbol {k}_{2D}=(k_{x},k_{y})$ and $k=|\boldsymbol {k}_{2D}|$ are the wavenumber vector and its magnitude, respectively, and $\Delta k=2{\rm \pi} /L$. The lateral and longitudinal axes are normalised by the kinematic viscosity $\nu$ and the energy dissipation rate

(4.3)\begin{equation} \varepsilon(z)=\frac{\nu}{2}{\left\langle {\left( \frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}} \right)}^{2} \right\rangle}_{xyt}, \end{equation}

which yield the Kolmogorov micro-scale length $\eta=(\nu^3/\varepsilon)^{1/4}$. The spectra of the 3-D steady solution and turbulent state are in good agreement at high wavenumber $k\eta \gtrsim 10^{0}$. Furthermore, in the wavenumber band $2{\rm \pi} /(L/4)\lesssim k\eta \lesssim 2{\rm \pi} /(L/16)$, corresponding to the intermediate-scale range, the energy spectrum follows Kolmogorov's $-5/3$ power law, $E=C_{K}\varepsilon ^{2/3}k^{-5/3}$ (Kolmogorov Reference Kolmogorov1941), with the constant $C_{K}\approx 1.5$, which is consistent with that in the inertial subrange of high-Reynolds-number turbulence (Sreenivasan Reference Sreenivasan1995; Ishihara et al. Reference Ishihara, Morishita, Yokokawa, Uno and Kaneda2016).

Figure 7. (a) Energy spectrum $E$ and (b) energy flux $\varPi$ at the centre of the fluid layer, $z=H/2$, in the 3-D steady solution (circles) and the turbulent state (lines) at $Ra=2.6\times 10^7$. The lateral and longitudinal axes are normalised by the kinematic viscosity $\nu$ and the energy dissipation rate $\varepsilon$ at $z=H/2$, which yield the Kolmogorov micro-scale length $\eta=(\nu^3/\varepsilon)^{1/4}$. The red dashed lines represent $E=1.5\varepsilon ^{2/3}k^{-5/3}$ and $\varPi /\varepsilon =1$, respectively. The light blue, green and light red colours indicate $k=2{\rm \pi} /(L/4)$, $2{\rm \pi} /(L/8)$ and $2{\rm \pi} /(L/16)$, respectively, normalised with $\eta$ in the 3-D steady solution, corresponding to the intermediate-scale structures shown in figure 6.

In figure 7(b), we show the energy flux in the wavenumber space, $\varPi (k,z)$ (Mizuno Reference Mizuno2016), defined as

(4.4)\begin{gather} \varPi(k,z)=\sum_{k'\ge k}\sum_{k-{\Delta k}/{2}<|{\boldsymbol{k}}_{2D}|<k+{\Delta k}/{2}}T^{s}({\boldsymbol{k}}_{2D},z), \end{gather}

(4.5)\begin{gather}T^{s}({\boldsymbol{k}}_{2D},z)=\mathrm{Re}\left[ {\left\langle \partial_{j}\tilde{u}_{i}{(\widetilde{u_{i}u_{j}})}^{{\dagger}ger} \right\rangle}_{t}-\frac{1}{2}\frac{\partial {\left\langle \tilde{u}_{j}{(\widetilde{u_{j}w})}^{{\dagger}ger} \right\rangle}_{t}}{\partial z} \right], \end{gather}

where $(\partial _{1},\partial _{2},\partial _{3})=(\textrm {i}k_{x},\textrm {i}k_{y},\partial /\partial z)$ and ${\dagger}ger$ denotes the complex conjugate. Here $T^{s}(\boldsymbol {k}_{2D},z)$ represents the energy transfer between the Fourier modes, and the sum of all spectral components does not contribute to the total energy budget, i.e. $\sum _{\boldsymbol {k}_{2D}}T^{s}({\boldsymbol {k}}_{2D},z)=0$. In the intermediate-scale range, the energy flux exhibits positive values, that is, the energy is transferred from large to small scale, and it scales with the same order of energy dissipation rate.

In the optimised 2-D steady solution (Sondak et al. Reference Sondak, Smith and Waleffe2015; Waleffe et al. Reference Waleffe, Boonkasame and Smith2015), as $Ra$ increases the horizontal scale of large-scale convective rolls reduces due to the shrinking of the horizontal period. The bulk flows at high $Ra$ diverge from those in the turbulent states, although it is intriguing that the optimised steady solutions exhibit the turbulent scaling $Nu\sim Ra^{0.31}$. Meanwhile, the 3-D steady solution reproduces not only the near-wall coherent structures but also the turbulence in the bulk of the Rayleigh–Bénard convection.

5. Summary and conclusions

We have discovered a 3-D steady solution to the Boussinesq equations, exhibiting a scaling of $Nu\sim Ra^{0.31}$ and leading to multi-scale coherent structures, which are similar to those observed in turbulent Rayleigh–Bénard convection. The invariant solution bifurcates from the conduction state at $Ra\sim 10^3$, and it has been tracked up to $Ra\sim 10^7$ using the Newton–Krylov iteration. The heat flux of the 3-D steady solution can surpass maxima in the 2-D steady solutions (Sondak et al. Reference Sondak, Smith and Waleffe2015; Waleffe et al. Reference Waleffe, Boonkasame and Smith2015), and thus a family of the 3-D steady solutions with optimised horizontal periods might achieve maximal heat transfer in steady Rayleigh–Bénard convection. In contrast, in the 3-D steady solution with fixed horizontal periods, the horizontal-averaged temperature and the RMS of the temperature and velocity fluctuations are in good agreement with the horizontal and temporal averages obtained for the turbulent states. In the near-wall region of both the 3-D steady solution and the turbulent states, smaller-scale thermal plumes emerge with an increase in $Ra$. The size of the coherent thermal structures and relevant vortices is comparable with the thermal conduction layer thickness $\delta /H=1/(2Nu)$, and the RMS vertical velocity at $z/\delta \sim 1$ scales with the velocity scale $Ra^{1/3}\kappa /H$, corresponding to $Nu\sim Ra^{1/3}$. The Reynolds number $Re$ based on the RMS vertical velocity scales as $Re\sim Ra^{4/9}$, which is related to the scaling $Nu\sim Ra^{1/3}$. The bulk flow in the 3-D steady solution comprises hierarchical multi-scale vortices. We have extracted the large- and intermediate-scale vortex structures by employing the coarse-graining method. The ratio of the largest to the smallest length scales in the 3-D steady solution at $Ra=2.6\times 10^{7}$ is approximately $20$. The energy spectrum at the centre of the fluid layer shows good agreement with that of the turbulent state. In the intermediate-scale range, the spectrum follows $E=1.5\varepsilon ^{2/3}k^{-5/3}$, which is commonly observed in the inertial subrange of the developed turbulence. Furthermore, energy is transferred from large to small scales in the wavenumber space, and the energy flux balances the energy dissipation rate, in accordance with the Kolmogorov–Obukhov energy cascade view.

Recently, van Veen, Vela-Martín & Kawahara (Reference van Veen, Vela-Martín and Kawahara2019) have found a time-periodic solution that reproduces inertial range dynamics in a triply periodic turbulence driven by a constant body force of the Taylor–Green type. They have obtained the invariant solution by applying large-eddy simulation based on the Smagorinsky-type eddy-viscosity model. Meanwhile, by introducing the buoyant force, we have succeeded in finding a multi-scale solution of the full incompressible Navier–Stokes equations without any empirical models. We believe that the current work and approaches based on multi-scale invariant solutions will trigger significant advances in the theoretical understanding and deductive modelling of coherent structures and energy transfer mechanisms in developed turbulence.