Heat transport by coherent Rayleigh-Bénard convection

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 chaos¹ and period doubling² to pattern formation³ and turbulence.⁴ A fundamental issue is to determine how the heat flux $H$ scales with the temperature difference ΔT between the bottom and top plates. There is a critical ΔT_c such that heat transport is by conduction with $H∼ΔT$ for ΔT < ΔT_c but macroscopic fluid motion develops for ΔT > ΔT_c and enhances heat transport. Bifurcations take place as ΔT is increased, leading to turbulent flow and $H∼ ( Δ T ) 4 / 3$ and perhaps even as high as $H∼ ( Δ T ) 3 / 2$⁠. Understanding these bifurcations and the actual asymptotic scaling of $H$ 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 $−g y ˆ$ is the constant acceleration of gravity. The fluid density is ρ ≈ ρ₀ (1 − α_VT), with α_V ≥ 0 the volumetric thermal expansion coefficient and ρ₀ 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 H³/(νκ) is larger than a critical value Ra_c, independent of the Prandtl number Pr = ν/κ.

We consider no-slip boundary conditions for which Ra_c ≈ 1708 for horizontal wavenumber⁵ α_c = 3.116/H. Convection develops and enhances heat transport for Ra > Ra_c. For statistically steady solutions, integration of (2) yields the heat flux,⁶ 

where $H 0 =κΔT/H$ is the heat flux in the absence of fluid motion and $v= y ˆ ⋅v$ 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.⁶ A simple version of this argument assumes an isothermal T = 0 interior with boundary layers of thickness δ such that⁷ Ra^(δ) ≡ gα_V(ΔT/2) δ³/(νκ) ≈ 1708/16, yielding heat flux

independent of Pr, where the Nusselt number $Nu=H/ H 0$⁠. This scaling also follows from assuming⁸ 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 $V= g ′ h$⁠, where g′ = gα_VΔT/2 is the reduced gravity. This yields the heat flux

This inertial scaling also follows from the standard turbulence assumption⁸ that heat transport becomes independent of ν and κ in the limit Vh/ν → ∞, Vh/κ → ∞.

Kraichnan’s mixing length theory⁹ yields various scalings in distinct regions of the (Ra, Pr) parameter space. He obtains Nu ∼ Ra^1/3, for fixed Pr, but predicts a transition to Nu ∼ Ra^1/2(lnRa)^−3/2 for very large Ra ≳ 10¹² when the shear boundary layers would be turbulent. The comprehensive scaling theory of Grossmann and Lohse¹⁰ 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.¹¹ They predict a transition to an “ultimate” regime similar to Kraichnan’s for Ra ≳ 10¹⁴ and argue that their log correction would yield an effective scaling¹² of about Nu ∼ Ra^0.38 in the range 10¹² < Ra < 10¹⁵. Some experimental results show transition to an ultimate regime^13–16 while others do not.^17,18 However, all the data are well-fitted by Nu ∼ 0.105 Ra^0.31 for Ra < 10¹¹ and that is the regime for which we report on unstable coherent states, in particular optimum transport solutions with Nu ∼ 0.115 Ra^0.31.

Rigorous upper bounds on heat transport for no-slip boundary conditions^19–23 yield Nu − 1 ≲ 0.026 Ra^1/2 as Ra → ∞, for any Pr, showing that free-fall scaling (5) cannot hold for large Pr. Rigorous bounds for free-slip²⁴ yield Nu − 1 ≲ 0.2295 Ra^5/12 that conflicts with (5) for any Pr (the bound has been improved to Nu − 1 ≲ 0.106 Ra^5/12).²⁹ 

We consider 2D flow with $v=u ( x , y , t ) x ˆ +v ( x , y , t ) y ˆ$ and use h = H/2, $V= g ′ h$⁠, τ = h/V, and ΔT/2 as our characteristic length, velocity, time, and temperature scales, respectively. Eliminating p by taking the $y ˆ$ component of the curl of (1) yields

where ν and κ are now non-dimensionalized by $Vh= g ′ h 3$ 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 $u ¯ ( y , t )$ 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 integration²⁵ in y with time-marching to steady states. Two distinct codes have been written; code 1 uses a 3rd order time-accurate scheme²⁶ and code 2 uses a semi-implicit forward-backward Euler time discretization together with Chebyshev tau of the 2nd kind²⁷ 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. 10⁶ (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 10⁶ with code 1 and Ra = 10⁹ with code 2 (symbols * in Fig. 2). The optimum solution at Ra = 10⁹ 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 Ra^0.28 according to a least square fit in 5.10⁵ ≤ Ra ≤ 5.10⁶. 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 Ra^0.217 (Fig. 3) to achieve an optimum heat flux Nu − 1 ≈ 0.115 Ra^0.31 (Fig. 2), from a least square fit in 10⁷ < Ra ≤ 10⁹. 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 fit¹⁸ Nu ≈ 0.088 Ra^0.32 and the pressurized SF₆ gas data fit¹⁶ Nu ≈ 0.105 Ra^0.312 (both in a cylinder of aspect ratio 1/2 and both for Ra < 10¹¹). The aspect ratio 4 data¹⁸ (Table 1) are well-fitted by Nu − 1 ≈ 0.102 Ra^0.31 in 10⁸ ≤ Ra ≤ 10¹⁰ and lie just below the optimum transport data with Nu − 1 ≈ 0.115 Ra^0.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 ≲ 10⁶ 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 ≈ 10⁴, but a spiral structure with slight downdraft of warm fluid (and updraft of cool) develops in 10⁴ ≲ Ra ≲ 10⁵, corresponding to the bump in the curve α = α^opt(Ra) in Fig. 3. The warm spiral begins winding back up at Ra ≈ 10⁵. 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 Ra^0.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.³⁰ 

Our results are connected with Malkus’ theory of turbulent convection⁶ 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).