Tailoring boundary geometry to optimize heat transport in turbulent convection
By tailoring the geometry of the upper boundary in turbulent Rayleigh-Bénard convection we manipulate the boundary layer – interior flow interaction, and examine the heat transport using the Lattice Boltzmann method. For fixed amplitude and varying boundary wavelength , we find that the exponent in the Nusselt-Rayleigh scaling relation, , is maximized at , but decays to the planar value in both the large () and small () wavelength limits. The changes in the exponent originate in the nature of the coupling between the boundary layer and the interior flow. We present a simple scaling argument embodying this coupling, which describes the maximal convective heat flux.
Thermal and compositional convection underlie the behavior of a wide range of systems from planetary and stellar interiors and the motions of Earth’s atmosphere and oceans, to the solidification of multicomponent melts Kadanoff (2001); Worster (2000); Wettlaufer (2011). The simplest setting to study thermal convection is in a Rayleigh-Bénard cell Chandrasekhar (2013), wherein the flow is controlled by the Rayleigh number (), which describes the ratio of buoyancy to dissipative forces, and the Prandtl number (), which is the ratio of momentum to thermal diffusivities of the fluid, and the aspect ratio ().
The key quantity of interest is the vertical heat flux across the cell, expressed in non-dimensional form as the Nusselt number, , which describes the ratio of the total heat flux to the heat flux solely due to conduction. For , the function is usually sought in the form of a scaling law: . The value emerges from the classical argument that when the dimensional heat flux should become independent of the depth of cell, implying that the boundary layers (BLs) at the upper and lower surfaces do not interact Priestley (1954); Malkus (1954); Howard (1966). However, the exponents obtained from experiments and numerical simulations range from approximately Castaing et al. (1989); Kerr (1996); Johnston and Doering (2009); Urban et al. (2011, 2012) to Urban et al. (2011, 2012); Niemela et al. (2000); Verzicco and Camussi (2003); Niemela and Sreenivasan (2006). Theories with specific assumptions concerning the structure of the flow Castaing et al. (1989) and/or the nature of the BLs Shraiman and Siggia (1990) have been proposed to explain the scaling. For extremely large , however, is predicted to become independent of the molecular properties of the fluid, and hence the boundary layers, and heat transport is achieved solely by convection Kraichnan (1962); Spiegel (1971). In this so-called “ultimate regime”, Spiegel (1971). Finally, we note that a means of examining the various “crossovers” in the plane has been proposed Stevens et al. (2013).
Taking a different approach, Howard Howard (1963) sought to determine upper bounds on using a variational formulation with incompressibility as one of the constraints on the statistically stationary flow. When a single horizontal wavenumber dominates the flow, he found an upper bound of . Kerswell Kerswell (1998) and Hassanzadeh et al. Hassanzadeh et al. (2014) (and references therein) provide a detailed discussion of this approach. A recent variational study of the two-dimensional problem by Whitehead & Doering Whitehead and Doering (2011) has shown that the thermal BLs (TBLs) do play a role in limiting the heat flux in the “ultimate regime”, even when there are no momentum boundary layers (free-slip conditions were used). Hence, the nature of the interaction between the BLs and the core flow plays a central role in turbulent Rayleigh-Bénard convection.
This interaction can be probed either by manipulating the boundary geometry itself or by introducing inhomogeneous temperature boundary conditions Ripesi et al. (2014). The former can be achieved by corrugating one or both horizontal boundaries, although we argue here that an asymmetric geometry provides unique insights. Earlier studies on convection over rough surfaces were motivated by the need for a better understanding of the role of BLs in the high regime Shen et al. (1996); Du and Tong (2000). The geometry used by Shen et al. Shen et al. (1996) and Du & Tong Du and Tong (2000) consisted of a cylindrical cell with rough top and bottom boundaries made of pyramidal elements. The ratio of wavelength () to amplitude of roughness () was fixed at (). When the thickness of the thermal boundary layer was smaller than , the pre-factor of the scaling relation increased. There was also an increase in the plume production with an enhanced detachment near the tip of the pyramids. Similar observations were made by Ciliberto & Laroche Ciliberto and Laroche (1999) in which the rough surfaces were prepared by gluing glass spheres to the bottom plate and coating the surface with a thermally conductive paint.
The experiments of Qiu et al. Qiu et al. (2005) and the direct numerical simulations of Stringano et al. Stringano et al. (2006) have shown that changes for periodic roughness, with in experiments and in the simulations. Wei et al.’s Wei et al. (2014) experiments with different combinations of smooth and rough surfaces (with pyramidal elements of ) at the top and bottom of the cell revealed that for: (a) both surfaces rough: , (b) only the top surface rough: and (c) only the bottom surface rough: . However, Dirichlet (Neumann) conditions are applied on the top (bottom) surface.
Similar studies have been carried out using rough surfaces made of rectangular elements Shishkina and Wagner (2011); Tisserand et al. (2011); Salort et al. (2014). Tisserand et al. Tisserand et al. (2011) investigated the effects of a bottom rough boundary on the heat transport near a planar top boundary by analyzing separately at the top and bottom boundary. They found that the smooth top boundary is not influenced by the effects of the bottom rough boundary.
Surface roughness has also been used in attempts to reach the ultimate regime at smaller than predicted by the theory of Kraichnan Kraichnan (1962). Roche et al. Roche et al. (2001) used a cylindrical cell with an interior entirely covered with V-shaped grooves to trigger a transition to turbulence in the BLs, and they reported for . Detailed accounts of the developments spanning various periods with a variety of perspectives can be found in a number of reviews Siggia (1994); Ahlers et al. (2009); Chillà and Schumacher (2012).
Despite differences in characteristics of turbulent flows in two and three dimensions (3D) (e.g., Boffetta and Ecke (2012)), numerical studies of 2D Rayleigh-Bénard convection for large and have yielded surprisingly close to those from experiments, the differences being principally in the pre-factors Johnston and Doering (2009); DeLuca et al. (1990); Waleffe et al. (2015). Clearly, this correspondence provides an essential role for well resolved 2D simulations to probe the properties of the key components (BLs, plumes, core flow) of convective flow and to relate them to the 3D dynamics for Schmalzl et al. (2004). However, we note that the pre-factors in 3D are larger and hence so too are the values of van der Poel et al. (2013).
Here, we describe quantitative studies of the effects of sinusoidal roughness of the upper boundary, with fixed amplitude and varying wavelength, on in two-dimensional Rayleigh-Bénard convection using highly resolved Lattice Boltzmann method (LBM) numerical simulations. We use this roughness to systematically probe the coupling between the BLs and the core flow and hence the resulting changes in the heat transport. We find that the heat transport is maximized at a dimensionless wavelength , with , and decays to the planar value () in both the large () and small () wavelength limits. This maximum originates in the nature of the coupling between the boundary layer and the bulk flow.
I Governing Equations and Numerical Method
We describe thermal convection with the Oberbeck–Boussinesq equations Chandrasekhar (2013). We non-dimensionalize them by choosing the height of the cell, , as the length scale, the temperature difference across this height, , as the temperature scale, as the velocity scale, where is the acceleration due to gravity and is the thermal expansion coefficient of the fluid, and as the time scale. The dimensionless equations of motion and boundary conditions are shown in figure 1. Here, , and are the velocity, temperature and pressure fields respectively, is the unit vector along the vertical axis, , and . The heat transfer rate in terms of can be obtained from , where denotes a horizontal and temporal average taken after the statistically steady state has been reached.
The governing equations are solved using the LBM Benzi et al. (1992); Chen and Doolen (1998); Succi (2001) with separate distribution functions for the momentum and temperature fields Shan (1997); Guo et al. (2002). For planar geometries with horizontal periodicity, one can solve the macroscopic equations numerically using spectral methods Orszag (1971), which provide a natural way to cluster grid points near the boundaries where steep gradients in temperature result for and , and hence where resolution is important. However, employing spectral methods for rough geometries is both technically challenging and computationally demanding, whereas the LBM can handle rough geometries naturally Chen and Doolen (1998); Succi (2001). Given our focus on the interaction between the boundary layer, which is ‘perturbed’ by the imposed roughness, and the core flow, it is advantageous to use the LBM. For all the simulations reported here, a mid-grid bounce-back condition is used to enforce the no-slip condition at the solid boundaries, and Dirichlet conditions for temperature Succi (2001). Periodic boundary conditions are used along the horizontal.
The code developed has been parallelized using the Message Passing Interface system, and extensively tested by reproducing a range of classical results from different flows Lipps (1976); Clever and Busse (1974); Rozhdestvensky and Simakin (1984); Johnston and Doering (2009). The simulations of Johnston & Doering Johnston and Doering (2009) with planar upper and lower plates constitute the most relevant Rayleigh-Bénard test comparison. They used a Fourier-Chebyshev spectral method with at least grid points inside the thermal boundary layers (TBLs). The aspect ratio – ratio of cell width () to height () – is fixed at . We use a uniform grid with at least grid points in the TBL, ensuring very high resolution throughout the domain. Figure 2 shows the comparison of our with that of Johnston & Doering Johnston and Doering (2009) for . We recover their results for , and their fit of the highest eight data points with . We note that this is also remarkably close to obtained experimentally by Urban et al. Urban et al. (2011) using cryogenic \ce^4He gas in a cylindrical cell of for and .
We now replace the upper planar boundary by a sinusoidal surface of wavelength and amplitude (figure 1). For all the simulations discussed here, , , and . We consider: , and . The characteristic vertical length scale used to define , and is . Simulations were for performed for for all . To give an example of the resolutions used: For the smallest () and highest () the number of grid points we used along the horizontal and vertical directions was and . The wavelength and amplitude of each roughness element was resolved using and grid points respectively. With the turnover time, the total run time for the highest cases was no less than , and data for statistics were collected only after . The run times were longer for smaller . Grid independence was confirmed from simulations with and for two grids: (a) and (b) . The difference in for the two simulations was % and the averaged temperature profiles were indistinguishable.
Figure 3 shows temperature fields for and and . Analysis of the data from all runs reveals that the number of plumes produced at the rough surface is a maximum for , demonstrating an enhanced interaction between the boundary layer and the core flow relative to the long- and short-wavelength cases. Namely, as increases or decreases relative to , this interaction weakens leading to a lower production of plumes. This effect can also be seen in the average temperature field within the well-mixed core region . For and we have as in the planar case, but for we find . This clearly shows the effect of the roughness element wavelength on the dynamics of the cold plumes that are released from the element tips. While previous experiments Shen et al. (1996); Du and Tong (2000); Qiu et al. (2005); Wei et al. (2014) and numerical simulations Stringano et al. (2006) have reported enhanced plume production, here we find a wavelength dependence, which exhibits a maximum plume production at a particular wavelength .
Because plumes are the conveyors of heat from one BL to another Howard (1966), this change in the dynamics has a direct effect on . Figure 4 shows as a function of . For each , was determined from a linear least squares fit to the vs. data for . Clearly, is maximal () at , and when or we recover the planar boundary value of . A dimensional argument providing inscribes a negatively buoyant parcel of vertical () to horizontal () aspect ratio to the upper sinusoidal grooves. However, there is a more complex dependence of the flux as described next.
Although figure 4 most clearly demonstrates the main point, figure 5 shows that the higher flux is dominated by the large contributions, where the heat flux increases as wavelength decreases. However, at lower the heat flux decreases as wavelength decreases. This wavelength dependence is demonstrated further in figure 6 where we plot the results of the following analysis. For a given we determine the maximum value of among all of the considered and we define this as . As shown in figure 5 the heat fluxes are larger at small (large) for large (small) and hence averages over these two competing wavelength dependent high and low flux behaviors giving a . This is the average over the large and small values seen in figure 4.
Finally, we comment on two experiments. Firstly, as we discussed above, surface roughness has been used in hopes of reaching the ultimate regime at smaller than predicted by the theory of Kraichnan Kraichnan (1962). We note that the effect of roughness enhancing transport is seen in figure 6 where we find for less than , which is very nearly that found by Urban et al. Urban et al. (2011) for greater than for planar surfaces. Secondly, we highlight the results for (or ) because of the close correspondence to one of the geometries used by Wei et al. Wei et al. (2014) in their experiments. In this case we obtain and they obtain . The agreement is remarkable given that they used a cylindrical cell of with a rough top plate made of pyramidal elements of , and highlights the importance of a systematic experimental studies in which the geometry of one surface is changed.
ii.2 Scaling Arguments
We now propose a simple scaling argument for the maximal exponent that has been observed here. We assume that: () The flow field is dominated by plumes. () The cold plumes that are generated at the rough boundary are due to negative buoyancy. () The flow in the core region is dissipation free, and hence the energy is dissipated only in the boundary layers. By way of analogy, we appeal to the Kolmogorov picture of 3-D turbulence to understand the kinetic energy transfer from the cold plumes to the momentum boundary layer (MBL) on the other side of the cell, which acts as a sink. A similar analogy has been used in the study of wall-bounded turbulent flows Sreenivasan (1989). Negative buoyancy ‘injects’ energy into the cold plumes at the upper rough boundary, which then travel through the core region into the MBL on the opposite side, where they are dissipated. The velocity scale of the plumes is , and the velocity scale in the MBL is , where is the dimensional thickness of the MBL. The scale for implies that inertial and viscous effects are of similar order close to the boundaries Schlichting and Gersten (2000); in other words in the BLs Zocchi et al. (1990). The length scale associated with the plumes is . Hence, the rate of injection of energy per unit mass is , and the energy dissipation rate per unit mass is . In the statistically steady state we must have . After some manipulation we obtain:
and noting that , where is the dimensional thickness of the TBL, and Howard (1966), we arrive at
The theoretical exponent () is close to that obtained from our numerical simulations (), and our scaling argument embodies the dynamics of an optimal–dissipation free core–flow interacting with the BLs. Hence, the small difference in is due to (a) the assumption that the core region is dissipation free (an obvious oversimplification) and (b) finite size effects in the simulations. Finally, although we are aware that the relation between and could perhaps be refined, it is not an essential issue for our argument given the rather weak dependence found here in equation 2.
We have investigated the effects of an upper rough boundary with fixed amplitude and varying wavelength on turbulent convection in two dimensions over a wide range of Rayleigh number . In this manner, we systematically manipulated the interaction between the boundary layer and interior/core flows. We have shown that there is a wavelength that maximizes heat transport, whereas for small () and large ( ) wavelengths the planar case (lower flux) results are recovered.
This wavelength dependence of the heat flux is reflected in the non-monotonic behavior of in the Nusselt-Rayleigh scaling relation, shown in figure 4. For , the boundary-layer/core-flow interaction is enhanced by an increase in the number of plumes produced along the roughness elements and their direct injection from the tips into the core flow, circumventing an intermediate transition region. The effect of the enhanced upper surface plume injection is to decrease the average core temperature relative to and .
The relation obtained for () is in agreement with the recent experiments of Wei et al. Wei et al. (2014), highlighting the importance of a systematic experimental study in which the geometry of one surface is changed.
A simple scaling argument describing the dynamics of a maximal flux has been proposed, and it is in good agreement with the simulation results. This prompts us to speculate on the relation between and variational approaches that seek maximal fluxes in the planar case for single wavenumber flows, such as Howard’s Howard (1963) treatment of optimal heat flux and Doering and Constantin’s Doering and Constantin (1996) upper bound using background test fields. In both cases, for fixed , which is similar to obtained here numerically. Clearly, detailed theory, numerics and experimentation is necessary for a firm understanding of this correspondence. Furthermore, we note that the robust upper bound scaling holds not only for flat no-slip boundaries, but also for both one- and two smoothly modulated boundaries (Goluskin & Doering, pers. comm.). In our ongoing but preliminary simulations over for both upper and lower boundaries having (), we find a pre-factor of 0.0055 (0.0091) and ().
Finally, given the fact that outside of the laboratory setting the boundaries of convecting fluids are rarely uniform, the results presented here have important implications for turbulent transport in astrophysical Mitra et al. (2013), engineering Chillà and Schumacher (2012) and geophysical Bercovici (2007) settings.
Acknowledgements.The authors acknowledge the support of the University of Oxford and Yale University, and the facilities and staff of the Yale University Faculty of Arts and Sciences High Performance Computing Center. ST acknowledges a NASA Graduate Research Fellowship and Sumesh P. T. for helpful discussions. JSW acknowledges the Swedish Research Council and a Royal Society Wolfson Research Merit Award for support. This work was completed at the 2015 Geophysical Fluid Dynamics Summer Study Program at the Woods Hole Oceanographic Institution, which is supported by the National Science Foundation and the Office of Naval Research. We thank many of the staff for feedback, particularly C.R. Doering & D. Goluskin.
- L. P. Kadanoff, Phys. Today 54, 34 (2001).
- M. G. Worster, in Perspectives in Fluid Dynamics — a Collective Introduction to Current Research, edited by G. Batchelor, H. Moffatt, and M. Worster (Cambridge University Press, 2000) pp. 393 – 446.
- J. S. Wettlaufer, Phys. Today 64, 66 (2011).
- S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability (Dover Publications, 2013).
- C. Priestley, Austr. J. Phys. 7, 176 (1954).
- W. V. R. Malkus, Proc. R. Soc. Lond. A 225, 196 (1954).
- L. N. Howard, in Applied Mechanics, Proc. of the 11th Congr. of Appl. Mech. Munich (Germany), edited by H. Görtler (Springer, 1966) pp. 1109–1115.
- B. Castaing, G. Gunaratne, F. Heslot, L. Kadanoff, A. Libchaber, S. Thomae, X.-Z. Wu, S. Zaleski, and G. Zanetti, J. Fluid Mech. 204, 1 (1989).
- R. M. Kerr, J. Fluid Mech. 310, 139 (1996).
- H. Johnston and C. R. Doering, Phys. Rev. Lett. 102, 064501 (2009).
- P. Urban, V. Musilová, and L. Skrbek, Phys. Rev. Lett. 107, 014302 (2011).
- P. Urban, P. Hanzelka, T. Kralik, V. Musilova, A. Srnka, and L. Skrbek, Phys. Rev. Lett. 109, 154301 (2012).
- J. Niemela, L. Skrbek, K. R. Sreenivasan, and R. J. Donnelly, Nature 404, 837 (2000).
- R. Verzicco and R. Camussi, J. Fluid Mech. 477, 19 (2003).
- J. Niemela and K. R. Sreenivasan, J. Fluid Mech. 557, 411 (2006).
- B. I. Shraiman and E. D. Siggia, Phys. Rev. A 42, 3650 (1990).
- R. H. Kraichnan, Phys. Fluids 5, 1374 (1962).
- E. A. Spiegel, Annu. Rev. Astron. Astrophys. 9, 323 (1971).
- R. J. Stevens, E. P. van der Poel, S. Grossmann, and D. Lohse, J. Fluid Mech. 730, 295 (2013).
- L. N. Howard, J. Fluid Mech. 17, 405 (1963).
- R. R. Kerswell, Physica D 121, 175 (1998).
- P. Hassanzadeh, G. P. Chini, and C. R. Doering, J. Fluid Mech. 751, 627 (2014).
- J. P. Whitehead and C. R. Doering, Phys. Rev. Lett. 106, 244501 (2011).
- P. Ripesi, L. Biferale, M. Sbragaglia, and A. Wirth, J. Fluid Mech. 742, 636 (2014).
- Y. Shen, P. Tong, and K.-Q. Xia, Phys. Rev. Lett. 76, 908 (1996).
- Y.-B. Du and P. Tong, J. Fluid Mech. 407, 57 (2000).
- S. Ciliberto and C. Laroche, Phys. Rev. Lett. 82, 3998 (1999).
- X.-L. Qiu, K.-Q. Xia, and P. Tong, J. Turb. 6, 1 (2005).
- G. Stringano, G. Pascazio, and R. Verzicco, J. Fluid Mech. 557, 307 (2006).
- P. Wei, T.-S. Chan, R. Ni, X.-Z. Zhao, and K.-Q. Xia, J. Fluid Mech. 740, 28 (2014).
- O. Shishkina and C. Wagner, J. Fluid Mech. 686, 568 (2011).
- J.-C. Tisserand, M. Creyssels, Y. Gasteuil, H. Pabiou, M. Gibert, B. Castaing, and F. Chilla, Phys. Fluids 23, 015105 (2011).
- J. Salort, O. Liot, E. Rusaouen, F. Seychelles, J.-C. Tisserand, M. Creyssels, B. Castaing, and F. Chilla, Phys. Fluids 26, 015112 (2014).
- P.-E. Roche, B. Castaing, B. Chabaud, and B. Hébral, Phys. Rev. E 63, 045303 (2001).
- E. D. Siggia, Ann. Rev. Fluid Mech. 26, 137 (1994).
- G. Ahlers, S. Grossmann, and D. Lohse, Rev. Mod. Phys. 81, 503 (2009).
- F. Chillà and J. Schumacher, Eur. Phys. J. E 35, 58 (2012).
- G. Boffetta and R. E. Ecke, Ann. Rev. Fluid Mech. 44, 427 (2012).
- E. DeLuca, J. Werne, R. Rosner, and F. Cattaneo, Phys. Rev. Lett. 64, 2370 (1990).
- F. Waleffe, A. Boonkasame, and L. M. Smith, Phys. Fluids 27, 051702 (2015).
- J. Schmalzl, M. Breuer, and U. Hansen, Europhys. Lett. 67, 390 (2004).
- E. P. van der Poel, R. J. Stevens, and D. Lohse, J. Fluid Mech. 736, 177 (2013).
- R. Benzi, S. Succi, and M. Vergassola, Phys. Rep. 222, 145 (1992).
- S. Chen and G. D. Doolen, Ann. Rev. Fluid Mech. 30, 329 (1998).
- S. Succi, The Lattice-Boltzmann Equation (Oxford Univ. Press, Oxford, 2001).
- X. Shan, Phys. Rev. E 55, 2780 (1997).
- Z. Guo, C. Zheng, and B. Shi, Phys. Rev. E 65, 046308 (2002).
- S. A. Orszag, Phys. Rev. Lett. 26, 1100 (1971).
- F. B. Lipps, J. Fluid Mech. 75, 113 (1976).
- R. Clever and F. Busse, J. Fluid Mech. 65, 625 (1974).
- B. Rozhdestvensky and I. Simakin, J. Fluid Mech. 147, 261 (1984).
- K. R. Sreenivasan, in Frontiers in Experimental Fluid Mechanics (Springer, 1989) pp. 159–209.
- H. Schlichting and K. Gersten, Boundary-layer theory (Springer, 2000).
- G. Zocchi, E. Moses, and A. Libchaber, Physica A: Statistical Mechanics and its Applications 166, 387 (1990).
- C. R. Doering and P. Constantin, Phys. Rev. E 53, 5957 (1996).
- D. Mitra, J. S. Wettlaufer, and A. Brandenburg, Astrophys. J. 773, 120 (2013).
- D. Bercovici, in Treatise on Geophysics, Vol. 7, edited by D. Bercovici (Elsevier, 2007) 2nd ed., Chap. 1, pp. 1–30.