Multi-scale steady solution for Rayleigh-B\'enard convection

We have found a multi-scale steady solution of the Boussinesq equations for Rayleigh-B\'enard convection in a three-dimensional periodic domain between horizontal plates with a constant temperature difference by using a homotopy from the wall-to-wall optimal transport solution given by Motoki et al. (J. Fluid Mech., vol. 851, 2018, R4). The connected steady solution, which turns out to be a consequence of bifurcation from a thermal conduction state at the Rayleigh number $Ra\sim10^{3}$, is tracked up to $Ra\sim10^{7}$ by using a Newton-Krylov iteration. The exact coherent thermal convection exhibits scaling $Nu\sim Ra^{0.31}$ (where $Nu$ is the Nusselt number) as well as multi-scale thermal plume and vortex structures, which are quite similar to those in the turbulent Rayleigh-B\'enard convection. The mean temperature profiles and the root-mean-square of the temperature and velocity fluctuations are in good agreement with those of the turbulent states. Furthermore, the energy spectrum follows Kolmogorov's -5/3 scaling law with a consistent prefactor, and the energy transfer to smaller scales in the wavenumber space agrees with the turbulent energy transfer.


Introduction
Rayleigh-Bénard convection, the buoyancy-driven flow in a horizontal fluid layer heated from below and cooled from above, is one of the most canonical flows widely observed in nature, including engineering materials. The effect of buoyancy on a flow is characterised by the Rayleigh number Ra. When Ra exceeds a certain critical value Ra c , the thermal conduction state becomes unstable, and two-dimensional (2D) steady convection rolls appear (Drazin & Reid 1981). At higher Ra, the convection becomes time-dependent, and subsequently exhibits turbulent states with multi-scale thermal and vortex structures. One of the primary interests in the Rayleigh-Bénard problem is the scaling of turbulent heat transfer with Ra, i.e., the dependence of the Nusselt number N u on Ra. Over half a century ago, Malkus (1954) derived the scaling N u ∼ Ra 1/3 by a marginal stability argument, based on the assumption that the thermal boundary layer adapts its thickness δ as δ/H ≈ (Ra/Ra c ) −1/3 , where H is the height of the fluid layer, so that the local Ra in the boundary layer becomes marginally stable. Subsequently, based on the mixing-length theory, Kraichnan (1962) predicted a transition of the boundary layer from laminar to turbulent state, and derived the asymptotic scaling, N u ∼ Ra 1/2 , with a logarithmic correction for very high Ra. The scaling N u ∼ Ra 1/2 is known as ultimate scaling, and has been obtained as the rigorous upper bound on N u by variational approaches (Doering & Constantin 1996;Plasting & Kerswell 2003). In conventional turbulent Rayleigh-Bénard convection, however, the ultimate scaling has not been observed yet. A prominent experiment by Niemela et al. (2000) for very high Ra exhibits N u ∼ Ra 0.31 even at Ra ∼ 10 17 . Grossmann & Lohse (2000) proposed a unifying scaling theory (GL theory) of global properties for Ra and the Prandtl number P r, based on decomposing the total scalar and energy dissipation into contributions from the bulk region and the boundary layer. A lot of experiments and numerical simulations have demonstrated the validity of this theory (Ahlers et al. 2009;Stevens et al. 2013). Per the theory, the scaling N u ∼ Ra 1/3 is derived in the high Ra regime 10 8 Ra 10 14 for P r ∼ 1. The transition to the ultimate scaling is also predicted for Ra 10 14 ; however, the effective scaling is N u ∼ Ra 0.38 due to logarithmic corrections (Grossmann & Lohse 2011). Although some results have shown the transition to N u ∼ Ra 0.38 , the high-Ra scaling is still being discussed (Chillà & Schumacher 2012;Zhu et al. 2018). On the other hand, for 10 8 Ra 10 11 , a lot of turbulent data exhibit N u ∼ Ra 0.31 (see, e.g. Niemela & Sreenivasan 2006;He et al. 2012).
Recently, Waleffe et al. Sondak et al. 2015) found a scaling N u = 0.115Ra 0.31 , which is quite similar to the turbulent data fit N u = 0.105Ra 0.312 (He et al. 2012) in 2D steady Rayleigh-Bénard convection for 10 7 Ra 10 9 . They obtained optimal 2D steady solutions to maximise N u by changing the horizontal periods, and the scaling was achieved by a family of 2D solutions with the horizontal period that decreases with increasing Ra. Although the result suggests that simple and coherent structures can capture the essence of turbulent convection, it does not imply that just any single 2D steady solution with a fixed horizontal period (maximal wavelength) can do it. More recently, the wall-to-wall optimal transport problem, which is a variational problem of finding a divergence-free velocity field optimising scalar transport between two parallel plates, has been discussed (Hassanzadeh et al. 2014;Tobasco & Doering 2017), and Motoki et al. (2018a) found three-dimensional (3D) steady velocity fields to be the optimal states maximising heat transfer between two isothermal no-slip parallel plates under the constraint of fixed total enstrophy. The optimal states exhibit ultimate scaling, which is quite close to the rigorous upper bound N u − 1 = 0.02634Ra 1/2 (Plasting & Kerswell 2003), as well as hierarchical self-similar vortex structures. In 2D velocity fields, however, such multi-scale structures have not been observed to be the optimal state (Souza et al. 2020). Although a 3D optimal state needs an external body force other than buoyancy, we proved that the optimal state can be continuously connected to a steady solution of the Boussinesq equation by using the homotopy continuation method (see appendix A). In the Rayleigh-Bénard convection between horizontal boundaries, a 3D steady solution with convection cells also bifurcates from the conduction state at the same critical value of Ra as the 2D steady solution (see the upper-left inset in Figure 1), as convection rolls in any horizontal direction can exist simultaneously. The connected solution is the 3D steady solution. Although this solution is not stable, it exists even at high Ra. In this paper, we demonstrate the ability of the invariant solution to capture key statistical features, as well as coherent thermal and flow structures in the turbulent Rayleigh-Bénard convection. We then discuss the hierarchical multi-scale vortex structures and energy transfer.
The remainder of the paper is organised as follows. In §2 we introduce the governing equations, the boundary conditions and the dimensionless parameters to characterise the thermal convection, and describe the numerical procedures to obtain nonlinear solutions. The statistical properties and spatial structures of the 3D steady solution are presented in §3, and the hierarchical vortex structures are discussed in §4. Finally, summary and conclusions are presented in §5. In the appendix A, we present the homotopy continuation analysis from the optimal solution of the Euler-Lagrange equations for the wall-to-wall optimal transport problem, to the present 3D steady solution of the Boussinesq equations. The parameter dependence of the 3D steady solution and the adequacy of the spatial resolution are shown in appendix B.

Boussinesq equations and numerical methods
Let us consider a fluid layer between two horizontal plates heated from below and cooled from above, and employ the Oberbeck-Boussinesq approximation, wherein the density variations are only significant in the buoyancy term. The time evolution of velocity field u(x, t) = ue x + ve y + we z and temperature field T (x, t) are described by the Boussinesq equations where p(x, t) is pressure, and ρ, ν, g, α and κ are the mass density, kinematic viscosity, acceleration due to gravity, volumetric thermal expansivity, and thermal diffusivity, respectively. e x and e y are mutually orthogonal unit vectors in the horizontal directions, while e z is a unit vector in the vertical direction. The two horizontal plates are positioned at z = 0 and z = H, and the top (or bottom) wall surface is no-slip and impermeable, and held at a lower (or higher) constant temperature: 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 P r, Ra = gα∆T H 3 νκ , P r = ν κ . (2.5) The vertical heat flux is quantified by the Nusselt number where · xyt and · 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 the equation (2.3). The equations (2.1)-(2.3) are discretised by 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. The aliasing errors are removed with the aid of the 2/3 rule and the 1/2 rule for the Fourier transform and the Chebyshev transform, respectively. 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 by the Newton-Krylov iteration (for more details, see §3 and appendix A in Motoki et al. 2018b).
In this paper, we present the steady solution and the turbulent states in the horizontally square periodic domain with L/H = π/2 ≈ 1.57 for P r = 1. The domain is the same as that of the optimal states derived by Motoki et al. (2018a). The numerical process is carried out on 128 3 grid points for Ra < 10 7 and 256 3 grid points for Ra 10 7 . In the domains with L/H = 2π/3.117 ≈ 2.02 and 1, and for P r = 7, we confirm that the effects of the domain size, P r and the spatial resolution on the heat flux at high Ra as well as the thermal and flow structures in 3D steady solutions, which will be discussed in the following sections, are insignificant in appendix B. All the 3D steady solutions presented in this paper satisfy the π/2-rotation symmetry as well as the mirror symmetry the shift-and-reflection symmetry Although these symmetries are not imposed explicitly, they are satisfied in the solutions at Ra 10 6 . At Ra ∼ 10 7 we directly impose the symmetries (2.7)-(2.9) in order to reduce the computational degrees of freedom.

Three-dimensional steady solution
3.1. N u-Ra scaling Figure 1 presents N u as a function of Ra. The red line shows the 3D 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 2006), respectively. The 3D steady solution maintains slightly larger N u than the turbulent states even at Ra ∼ 10 7 . The bottom right inset shows N u compensated by Ra γ , and the scaling exponent γ of turbulent states shows 2/7 for Ra 10 7 , and at higher Ra it changes to 0.31. Such a transition has been experimentally and numerically observed for P r ∼ 1 (Castaing et al. 1989;Silano et al. 2010). Meanwhile, the exponent of heat flux in the 3D steady solution is greater than 2/7 but less than 1/3, and it can be approximated to N u − 1 = 0.115Ra 0.31 (blue dashed, Waleffe et al. 2015;Sondak et al. 2015), which is achieved by a family of 2D steady solutions with optimal horizontal periods. The orange dashed line indicates the optimal scaling N u − 1 = 0.0236Ra 1/2 (Motoki et al. 2018a) given by the 3D optimal states in the wall-to-wall optimal transport problem, and it is quite close to the rigorous upper bound (Plasting & Kerswell 2003). Although optimal states exhibiting significantly high heat flux have been achieved by external body force being different from buoyancy, the steady solution can be continuously connected to the present 3D steady solution of the full Boussinesq equations by a homotopy from the body force to the buoyancy, as shown in appendix A.

Mean temperature and root-mean-square profiles
The 3D steady solution reproduces the mean temperature of turbulent states in the whole region. Furthermore, the root-mean-square (RMS) values are also in good agreement with each other (figure 2). Note that the RMS values are obtained from the horizontal average for the steady solutions, and the time and horizontal averages for the turbulent states. In the bulk region, all mean temperature profiles are flattened,  (Niemela & Sreenivasan 2006). The blue dashed line indicates the optimal scaling in the 2D steady solutions, N u − 1 = 0.115Ra 0.31 Sondak et al. 2015). The orange solid and dashed lines indicate the upper bound N u − 1 = 0.02634Ra 1/2 (Plasting & Kerswell 2003) and the optimal scaling N u − 1 = 0.0236Ra 1/2 , respectively, evaluated from the wall-to-wall optimal transport states (Motoki et al. 2018a). The black curves in the top-left inset show the maximal and minimum values of N u in the 3D time-periodic solution. The bottom right inset shows N u compensated by Ra γ : γ = 2/7 (plot A); γ = 0.31 (plot B); γ = 1/3 (plot C). as a result of the nearly complete mixing by large-scale convection. The temperature difference ∆T /2 exists only at the thermal conduction layer, 0 z 2δ/H = 1/N u, and T ′ RMS has peaks at z/δ ≈ 1. If the advection, diffusion, and buoyancy terms in the Navier-Stokes equation (2.2) at the conduction layer are balanced as (the balance between the advection and diffusion terms is given by that in the energy equation (2.3) for ν ∼ κ) then the near-wall vertical velocity would be and yields the scaling law which has been given by Malkus' theory (Malkus 1954) and the GL theory (Grossmann & Lohse 2000). As shown in figure 2(f ) the RMS vertical velocity w ′ RMS scales as Ra 1/3 κ/H near the wall z/δ ∼ 1.   with the thickness δ. In 2D steady solutions Sondak et al. 2015), the appearance of such small-scale plume and vortex structures has not been observed for a fixed horizontal period, and the scaling N u ∼ Ra 0.31 is achieved by a family of solutions with smaller horizontal periods as the Ra increases. It should be stressed that the single 3D steady solution spontaneously reproduces the multi-scale coherent structures of convective turbulence.

Hierarchical vortices and energy transfer in wavenumber space
The where ∆r = |x ′ − x|, σ 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 ( Figure 4(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  5(a) shows the energy spectrum E(k, z) of the 3D steady solution at the centre of the fluid layer, z = H/2, and the corresponding turbulence spectrum at Ra = 2.6 × 10 7 . E(k, z) is defined as where (·) indicates the Fourier coefficients in the periodic (x-and y-) directions, and · t represents the time average. k 2D = (k x , k y ) and k = |k 2D | are the wavenumber vector and its magnitude, respectively, and ∆k = 2π/L. The lateral and longitudinal axes are normalised by the kinematic viscosity ν and the energy dissipation rate   (Kolmogorov 1941), with the constant, C K ≈ 1.5, which is consistent with that in the inertial subrange of high-Reynolds-number turbulence (Sreenivasan 1995;Ishihara et al. 2016).
In figure 5(b), we show the energy flux in the wavenumber space, Π(k, z) (Mizuno 2016), defined as where (∂ 1 , ∂ 2 , ∂ 3 ) = (ik x , ik y , ∂/∂z) and † denotes the complex conjugate. T s (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., k2D T s (k 2D , z) = 0. In the intermediate-scale range, the energy flux exhibits positive values, that is, the energy transfer from large to small scale, and it scales with the same order of energy dissipation rate.

Summary and conclusions
We have discovered a three-dimensional steady solution to the Boussinesq equations that exhibits scaling (N u ∼ Ra 0.31 ) and 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 ∼ 10 3 , and it has been tracked up to Ra ∼ 10 7 by using the Newton-Krylov iteration. The horizontal-averaged temperature and the RMS of the temperature and velocity fluctuations are in good agreement with the horizontal and temporal averages for the turbulent states. In the near-wall region, smaller-scale thermal plumes are generated with an increase in Ra. The size of the thermal coherent structures and relevant vortices is comparable with the thermal conduction layer thickness δ/H = 1/(2N u), and the RMS vertical velocity at z/δ ∼ 1 scales with the velocity scale Ra 1/3 κ/H, corresponding to N u ∼ Ra 1/3 . On the other hand, in the bulk region, the flow consists of 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 3D steady solution at Ra = 2.6 × 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ε 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 et al. (2019) 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. By introducing the buoyant force, meanwhile, we have succeeded in finding a multi-scale solution of the full incompressible Navier-Stokes equation 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.
where φ ′ (x ′ ), ψ ′ (x ′ ) and µ ′ are Lagrange multipliers, and prime (·) ′ represents a nondimensional variable based on H, ∆T , κ and ρ. The Euler-Lagrange equations are In our previous work (Motoki et al. 2018a), we obtained the optimal state so as to satisfy the equations (A 3)-(A 7). Thus, to fulfil the Boussinesq equations, the optimal velocity and temperature field (u ′ op , θ ′ op ) require an additional body force which is different from the buoyant force, where p ′ op is the pressure determined by the Poisson equation stemming from the Boussinesq equations. We consider homotopy from the Euler-Lagrange system to the steady Boussinesq system where α is a homotopy parameter. For a fixed Ra = 10 4 , P r = 1 and f ′ , we have tracked the solution from α = 1 to 0 by using the Newton-Krylov method (figure 6). The connected solution S 3D is the present three-dimensional steady solution shown in §3 and §4. Appendix B. Dependence of multi-scale steady solution on domain size, Prandtl number and spatial resolution Figure 7 shows the Nusselt number compensated by Ra 0.31 as a function of the Rayleigh number Ra in the three-dimensional steady solutions for the different horizontal period L and Prandtl number P r. The green, red and light blue symbols show the L/H = 2π/(k c H) ≈ 2.02, π/2 ≈ 1.57 and 1, respectively, for P r = 1, and the light red symbols represent L/H = π/2 for P r = 7, where k c = 3.117/H is the wavenumber corresponding to the minimal critical Ra c = 1708 (Drazin & Reid 1981). Although we observe the dependence of N u on Ra, at Ra 10 7 all the plots exhibit the values approximate to the optimal scaling N u − 1 = 0.115Ra 0.31 (blue dashed, Waleffe et al. 2015;Sondak et al. 2015) in the two-dimensional steady solutions. We expect that the variation in the domain size and P r (for 1 P r 10) of N u would not be significant at high Ra, since the emergence of the small-scale plume and vortex structures near the walls, which are robustly observed in different L and P r ( figure 3, 8 and 9), might be key ingredient in the vertical heat flux.
In figure 7 the symbols +, •, × and • show the results obtained on different grid points (N x , N y , N z ) = (64, 64, 64), (128,128,128), (192,192,128) and (256,256,256), respectively, and the effects of the spatial resolutions on the N u are minor. For our main results with L/H = π/2 and P r = 1, the grid points (N x , N y , N z ) = (128,128,128) are enough to evaluate the characteristics of the 3D steady solution at Ra 10 7 ; (N x , N y , N z ) = (256, 256, 256) are sufficient at Ra ∼ 10 7 . Furthermore, the Kolmogorov micro-scale length η and the thermal conduction layer thickness δ in the 3D steady solution and the turbulent states at Ra = 10 5 , 10 6 , 10 7 and 10 7.42 ≈ 2.6 × 10 7 are shown in table 1 together with the grid sizes. Since the energy dissipation rate is a function of the wall-normal coordinate z, ε(z) = (ν/2) (∂u i /∂x j + ∂u j /∂x i )  Table 1: Numerical details of the 3D steady solution and the turbulent states for P r = 1 and L/H = π/2. ∆x and ∆z are the spatial resolutions in the x-and z-directions. η z and η| c represent the Kolmogorov micro-scale length η = (ν 3 /ε) 1/4 based on the vertical averaged energy dissipation rate, ε z = (ν/2) (∂u i /∂x j + ∂u j /∂x i ) 2 xyzt , and that at the centre of the fluid layer, ε| z=H/2 , respectively. δ is the thermal conduction layer thickness, δ/H = 1/(2N u). τ is the integral time to obtain the statistics, and U = √ gα∆T H is the buoyancy-induced terminal velocity. and less than one third of δ. Therefore, the spatial resolution is sufficient to resolve the smallest-scale thermal and flow structures in the 3D steady solution and the turbulent states. The present turbulent DNS data is obtained by averaging time of more than 200 convective time units based on the buoyancy-induced terminal velocity U = √ gα∆T H.