Steady but generally unstable solutions of the 2D Boussinesq equations are obtained for no-slip boundary conditions and Prandtl number 7. The primary solution that bifurcates from the conduction state at Rayleigh number Ra ≈ 1708 has been calculated up to Ra ≈ 5.106 and its Nusselt number is Nu ∼ 0.143 Ra0.28 with a delicate spiral structure in the temperature field. Another solution that maximizes Nu over the horizontal wavenumber has been calculated up to Ra = 109 and scales as Nu ∼ 0.115 Ra0.31 for 107 < Ra ≤ 109, quite similar to 3D turbulent data that show Nu ∼ 0.105 Ra0.31 in that range. The optimum solution is a simple yet multi-scale coherent solution whose horizontal wavenumber scales as 0.133 Ra0.217. That solution is unstable to larger scale perturbations and in particular to mean shear flows, yet it appears to be relevant as a backbone for turbulent solutions, possibly setting the scale, strength, and spacing of elemental plumes.
Rayleigh-Bénard convection is the buoyancy-driven motion of a fluid contained between two horizontal plates and heated from below. It is a paradigmatic problem with a rich array of nonlinear physics from deterministic chaos1 and period doubling2 to pattern formation3 and turbulence.4 A fundamental issue is to determine how the heat flux scales with the temperature difference ΔT between the bottom and top plates. There is a critical ΔTc such that heat transport is by conduction with for ΔT < ΔTc but macroscopic fluid motion develops for ΔT > ΔTc and enhances heat transport. Bifurcations take place as ΔT is increased, leading to turbulent flow and and perhaps even as high as . Understanding these bifurcations and the actual asymptotic scaling of have been the subjects of many experimental, theoretical, and numerical studies focusing on turbulent convection. Here, we consider coherent convection—steady flows that may be the permanent form of the coherent structures (“plumes,” “thermals,” and large scale “winds”) permeating turbulent convection.
In Rayleigh-Bénard convection, the fluid is contained between two infinite horizontal plates and the Boussinesq approximation is made so the governing equations are
for the velocity v and the temperature T. Incompressibility ∇ ⋅ v = 0 is maintained by the kinematic pressure p and is the constant acceleration of gravity. The fluid density is ρ ≈ ρ0 (1 − αVT), with αV ≥ 0 the volumetric thermal expansion coefficient and ρ0 a constant reference density, ν > 0 is the kinematic viscosity, and κ > 0 is the thermal diffusivity. The bottom plate temperature is fixed at ΔT/2 and the upper plate is at −ΔT/2. The conduction state is the solution v = 0, T = − yΔT/H to (1) and (2), with −h ≤ y ≤ h = H/2. This is a steady solution for all values of the parameters, but Rayleigh showed that it is linearly unstable when the Rayleigh number Ra = gαVΔT H3/(νκ) is larger than a critical value Rac, independent of the Prandtl number Pr = ν/κ.
We consider no-slip boundary conditions for which Rac ≈ 1708 for horizontal wavenumber5 αc = 3.116/H. Convection develops and enhances heat transport for Ra > Rac. For statistically steady solutions, integration of (2) yields the heat flux,6
where is the heat flux in the absence of fluid motion and is the vertical velocity with and 〈∗〉 denoting horizontal and domain averages, respectively.
A classic scaling argument is that heat transport is determined by marginal stability of the thermal boundary layers and of the mean temperature gradient in the interior.6 A simple version of this argument assumes an isothermal T = 0 interior with boundary layers of thickness δ such that7 Ra(δ) ≡ gαV(ΔT/2) δ3/(νκ) ≈ 1708/16, yielding heat flux
independent of Pr, where the Nusselt number . This scaling also follows from assuming8 that the heat flux becomes independent of the height H as Ra → ∞.
Another simple argument imagines a fluid at temperature T = 0, with cold plumes at temperature −ΔT/2 free-falling a distance H = 2h from top to bottom plates (and hot ΔT/2 plumes “free-rising” from bottom to top) at average speed , where g′ = gαVΔT/2 is the reduced gravity. This yields the heat flux
This inertial scaling also follows from the standard turbulence assumption8 that heat transport becomes independent of ν and κ in the limit Vh/ν → ∞, Vh/κ → ∞.
Kraichnan’s mixing length theory9 yields various scalings in distinct regions of the (Ra, Pr) parameter space. He obtains Nu ∼ Ra1/3, for fixed Pr, but predicts a transition to Nu ∼ Ra1/2(lnRa)−3/2 for very large Ra ≳ 1012 when the shear boundary layers would be turbulent. The comprehensive scaling theory of Grossmann and Lohse10 incorporates classic Kolmogorov scaling of velocity and temperature dissipation rates in the bulk together with thermal and viscous boundary layers. Their theory fits many experimental data sets.11 They predict a transition to an “ultimate” regime similar to Kraichnan’s for Ra ≳ 1014 and argue that their log correction would yield an effective scaling12 of about Nu ∼ Ra0.38 in the range 1012 < Ra < 1015. Some experimental results show transition to an ultimate regime13–16 while others do not.17,18 However, all the data are well-fitted by Nu ∼ 0.105 Ra0.31 for Ra < 1011 and that is the regime for which we report on unstable coherent states, in particular optimum transport solutions with Nu ∼ 0.115 Ra0.31.
Rigorous upper bounds on heat transport for no-slip boundary conditions19–23 yield Nu − 1 ≲ 0.026 Ra1/2 as Ra → ∞, for any Pr, showing that free-fall scaling (5) cannot hold for large Pr. Rigorous bounds for free-slip24 yield Nu − 1 ≲ 0.2295 Ra5/12 that conflicts with (5) for any Pr (the bound has been improved to Nu − 1 ≲ 0.106 Ra5/12).29
We consider 2D flow with and use h = H/2, , τ = h/V, and ΔT/2 as our characteristic length, velocity, time, and temperature scales, respectively. Eliminating p by taking the component of the curl of (1) yields
where ν and κ are now non-dimensionalized by and relate to the Rayleigh and Prandtl numbers as
All results in this paper are for Pr = 7 (water). The horizontal velocity u is obtained from v and ∂xu + ∂yv = 0 except for its horizontal average that is determined by (8), where the overbar denotes a horizontal average. Equations (6)–(8) are considered with no-slip boundary conditions,
together with periodicity of period L = 2π/α in the x direction.
We look for steady solutions, ∂t = 0, that bifurcate from the conduction state at Ra ≈ 1708 for wavenumber α ≈ 3.116/2 corresponding to aspect ratio L/H ≈ 2.016. Those convective solutions obey mirror symmetry,
as well as shift-reflect symmetry,
Equations (6) and (7) are discretized using a Fourier expansion in x and Chebyshev integration25 in y with time-marching to steady states. Two distinct codes have been written; code 1 uses a 3rd order time-accurate scheme26 and code 2 uses a semi-implicit forward-backward Euler time discretization together with Chebyshev tau of the 2nd kind27 for the v equation. Code 2 typically uses inconsistent time integration, with smaller time steps for smaller wavenumbers and for the v equation than for the T equation, to speed up or enable convergence to steady state. Both codes are dealiased in x and y using the 2/3 rule. The results of both codes overlap or connect very well as seen in Figs. 1–3, discussed below. The results also match a 3rd code based on finite differences and a Newton-Krylov iteration to steady state (Sondak, Smith, Waleffe 2015, not shown here). Mirror symmetry (11) is imposed and eliminates mean flow (8). Shift-reflect symmetry (12) is not imposed explicitly but is satisfied by the steady solutions presented here.
We show two branches of nonlinear steady solutions that bifurcate from the conduction state at Ra ≈ 1708, α ≈ 1.558. The primary branch (PB) has fixed horizontal wavenumber α (rounded here to α = π/2⇔L/H = 2). The optimum branch (OB) adjusts α to maximize heat flux Nu. The (Ra, Nu) curves for both solutions are shown in Figures 1 and 2. The steady primary solution is stable up to Ra ≈ 53 000 where it spawns a time-periodic solution. The maximum and minimum Nu achieved by that periodic solution are shown in Fig. 1; that time periodic solution was computed with the time-accurate code 1. The unstable steady state was continued with code 2 up to Ra = 5. 106 (symbols ∘ in Fig. 2). Mirror symmetry (11) suffices to stabilize the optimum solution in its fundamental periodic domain. That optimum solution was calculated up to Ra = 1.4 106 with code 1 and Ra = 109 with code 2 (symbols * in Fig. 2). The optimum solution at Ra = 109 is well-resolved with 200 Chebyshev polynomials in y and 200 Fourier modes in x, after dealiasing.
The heat flux for the unstable primary solution scales as Nu − 1 ≈ 0.143 Ra0.28 according to a least square fit in 5.105 ≤ Ra ≤ 5.106. That solution develops a delicate spiral structure in the temperature field but not in the velocity (Fig. 4). The winding of these temperature spirals continuously increases with Rayleigh number. Convergence of our algorithm to that unstable steady solution becomes quite slow for increasing Ra, apparently because of the center region of these spiral structures where v, T ≈ 0.
The optimum solution increases wavenumber α with Ra as α ≈ 0.133 Ra0.217 (Fig. 3) to achieve an optimum heat flux Nu − 1 ≈ 0.115 Ra0.31 (Fig. 2), from a least square fit in 107 < Ra ≤ 109. This optimum heat flux scaling is quite similar to the heat flux observed in 3D turbulent convection experiments. Indeed, the optimum branch in Fig. 2 is only slightly above the cryogenic helium gas data fit18 Nu ≈ 0.088 Ra0.32 and the pressurized SF6 gas data fit16 Nu ≈ 0.105 Ra0.312 (both in a cylinder of aspect ratio 1/2 and both for Ra < 1011). The aspect ratio 4 data18 (Table 1) are well-fitted by Nu − 1 ≈ 0.102 Ra0.31 in 108 ≤ Ra ≤ 1010 and lie just below the optimum transport data with Nu − 1 ≈ 0.115 Ra0.31. This close agreement is remarkable given the differences in dimension, 2D optimum vs. 3D data, and Prandtl number, Pr = 7 for the 2D optimum vs. Pr ≈ 0.7 in the helium gas experiments.
The optimizing wavenumber α = αopt(Ra) (Fig. 3) shows an undulation between 1708 < Ra ≲ 106 that is linked to a temperature spiral structure in the optimum transport solution as well (Fig. 5). A temperature updraft of hot fluid (and downdraft of cold fluid) develops between Ra = 1708 and Ra ≈ 104, but a spiral structure with slight downdraft of warm fluid (and updraft of cool) develops in 104 ≲ Ra ≲ 105, corresponding to the bump in the curve α = αopt(Ra) in Fig. 3. The warm spiral begins winding back up at Ra ≈ 105. The spiral does not appear to ever wind back down for higher Ra, developing instead an increasingly jagged structure (Fig. 5, right), and αopt approaches the power law scaling αopt ≃ 0.133 Ra0.217. The structure of optimum Nu solutions depends on Prandtl number Pr. This is investigated in forthcoming work which also shows that the locally optimum steady solution discussed here is in fact the global optimum over α.
The close agreement between the optimum transport 2D steady state solutions and 3D turbulent data in Fig. 2 is intriguing. Although this agreement could be fortuitous, it strongly suggests that a single unstable steady solution may capture key statistical features of fully developed turbulent flows, such as the net heat flux and the mean temperature profile as well as the strength and scale of elemental plumes, as in shear flows.25,28 A search for such maximum momentum transport solutions in shear flows was initiated over 10 yr ago but not completed because of the higher computational complexity of those 3D solutions.30
Our results are connected with Malkus’ theory of turbulent convection6 and subsequent work on upper bounds.19–24 Two key ingredients of Malkus’ theory are maximum heat transport and marginal stability of the mean temperature profile and the smallest scales of motion. Both ingredients are included in our calculations that maximize heat transport over horizontal wavenumber α and track steady state solutions that bifurcate from the marginal stability critical point at Ra ≈ 1708, α ≈ 1.558. We conjecture that our 2D results are in fact the 3D optimum transport solutions (for infinite or periodic horizontal directions). The upper bound results also assume 2D optimizers.
Why the optimum transport solution should capture gross 3D turbulence characteristics might be understood as a “winner-take-all” effect, where the optimum solution consumes all available potential energy so no other flow can be sustained. The optimum solution appears to be stable when mirror symmetry (11) is imposed and length scales are restricted to be less or equal to the optimum wavelength. The optimum solution is unstable to larger scale perturbations, in particular to subharmonics where plumes merge and form bigger plumes in a cyclic or quasi-cyclic fashion. It is also unstable to mean shear flow (8) when mirror symmetry is allowed to be broken. These instabilities would be the source of the “turbulence” but the underlying unstable coherent solutions control the heat transport.
The authors would like to thank Charles Doering (University of Michigan) and David Sondak (University of Wisconsin-Madison) for helpful comments and discussions and Detlef Lohse (University of Twente) for providing many recent references on the ultimate regime and the Grossmann-Lohse scaling theory. This research was partially supported by NSF Grant Nos. DMS-0807349 (FW, AB) and DMS-1008396 (AB, LMS).
Estimated from linear stability about the conduction state but with free “outflow” boundary conditions at the top edge of the boundary layer, that is, .