Vibration-induced boundary-layer destabilization achieves massive heat-transport enhancement.

Thermal turbulence is well known as a potent means to convey heat across space by a moving fluid. The existence of the boundary layers near the plates, however, bottlenecks its heat-exchange capability. Here, we conceptualize a mechanism of thermal vibrational turbulence that breaks through the boundary-layer limitation and achieves massive heat-transport enhancement. When horizontal vibration is applied to the convection cell, a strong shear is induced to the body of fluid near the conducting plates, which destabilizes thermal boundary layers, vigorously triggers the eruptions of thermal plumes, and leads to a heat-transport enhancement by up to 600%. We further reveal that such a vibration-induced shear can very efficiently disrupt the boundary layers. The present findings open a new avenue for research into heat transport and will also bring profound changes in many industrial applications where thermal flux through a fluid is involved and the mechanical vibration is usually inevitable.

INTRODUCTION

Heat transported from a solid surface to a fluid is a common scene that one often confronts in many natural and industrial processes. Turbulent thermal convection is known for its ability to vigorously transport heat (, ), and it ubiquitously occurs in nature. For instance, in stars like the Sun and in planets like our Earth, convection can be found in the atmosphere () and oceans (), in the mantle, and also in its core, which is considered to be responsible for the geomagnetic field (). It can also be found in the daily life settings, such as the indoor air circulation in a house, and in many engineering applications, such as the production of semiconductor crystals. The canonical setup for the study of convective heat transport is the so-called Rayleigh-Bénard (RB) convection, i.e., a fluid layer heated from below by a horizontal hot plate and cooled from above by a cooling plate. It is one of the most classical paradigms in fluid mechanics ().

As efficient heat transport by passive means is usually essential for many thermal industrial applications, how to enormously promote the heat-transfer efficiency holds a central issue in the community of thermal turbulence. Recently, many strategies (, ) were proposed to achieve high heat flux. The most widely used method is to disturb the boundary-layer structures, such as by creating roughness on the conducting plates (–), which enables a notable increase of heat flux (, ) but may suppress the global heat transport in some parameter ranges (). Geometrical confinement provides an additional way to enhance heat transfer via plume condensation (, ), despite an intermediate range of the cell aspect ratio only. Very recently, the combination of inclination of the convection cell and confined geometries was found to be able to lead to an increase of heat transport for several times, especially for low Prandtl numbers and relatively low Rayleigh numbers (, ). Up to now, how to severely improve the efficiency of the global heat transport through the RB convection is still highly desirable, which is the main goal of the present investigation.

In the classical regime of thermal turbulence, the process of convective heat transport is tremendously restricted by thermal diffusion in the laminar boundary layers, which impedes the effective heat exchange between the turbulent bulk flow and boundary layers. The realization of efficient heat transfer thus relies on the minimization of the roles of boundary layers, but manipulating boundary layers is notoriously difficult in thermal turbulence, partly due to the weak shear generated by the large-scale mean flow near the conducting plates. As an example, the introduction of wall roughness has been widely adopted to disrupt boundary layers by enhancing the detachment of the thermal boundary layer from the tips of rough elements (–). As the analog to pressure drag is absent in the temperature advection equation (), however, the system settles back to the boundary-layer–controlled regime for large imposed temperature differences (, ). Another attempt is to perform the lateral confinement of turbulent flows (, ), which has been successful in manipulating the flow structures in the convective bulk, yet the laminar boundary layers still limit the global heat flux. Therefore, the only way for massively enhancing heat transport is based on the breakthrough of the boundary-layer limitation.

To achieve this, we here put forward the concept of thermal vibrational turbulence via the introduction of horizontal vibration to thermal turbulence. Vibration is inevitable in almost all mechanical devices, and it can even be found in the form of Earth’s free oscillations after large earthquakes (). In the context of fluid motion, previous results have shown that under microgravity conditions, vibration offers a mechanism to develop mean flows in a nonuniformly heated fluid to transport heat and mass (, ), but the presence of gravity markedly breaks down the convective flow induced by external vibration (). The influence of vibration on the propagation of a temperature wave from a heat source in near-critical SF6 was also investigated under weightlessness (). On the other hand, under normal gravity conditions, Gershuni et al. () theoretically analyzed the linear stability of double diffusive convection in the presence of high-frequency vibration and determined the conditions of quasi-equilibrium. It was further found that vibration can generate or delay convective instabilities depending on the mutual orientation of the vibration axis and gravity (, ). Specifically, vertical vibration that is parallel to the gravity has been used to stabilize the unstable thermal gradients (), even in the turbulent regime (, ).

A central discovery of the present investigation is that when horizontal vibration that is perpendicular to the gravity is applied to thermal convection in its turbulent regime, high-frequency vibration vigorously destabilizes the boundary layers and thus can trigger massive emissions of thermal plumes that notably contribute to the vertical heat flux. The corresponding heat-transfer efficiency can be enhanced by up to 600%, indicating a promising technological potential once implemented in practice.

RESULTS

Heat-transport enhancement

We carry out three-dimensional (3D) direct numerical simulations (DNSs) of turbulent RB convection under the Oberbeck-Boussinesq approximation in a rectangular cell of height H, length L, and width W. As shown in Fig. 1F, horizontal vibration in the x direction with displacement amplitude δ and angular frequency ω, i.e., δcos (ωt), is applied to the convection cell. Further details about the simulations are given in Materials and Methods.

Fig. 1

Turbulent thermal vibrational convection.

(A) Instantaneous flow structures observed under different vibration frequencies with ω = 0 (left), ω = 700 (middle), and ω = 1400 (right). The 3D flow is illustrated by temperature isosurfaces with low (green) and high (brown) temperature, and the opacity is set to be 50% (see supplementary movies). For all the three cases, Ra = 108, Pr = 4.38, and δ = 0.1. Here, δ and ω are the dimensionless vibration amplitude and angular frequency with respect to the cell height H and the free-fall time scale , respectively. (B and C) Corresponding temperature contours (B) on the horizontal slice at z = δ(ω) and (C) on the vertical slice at y/W = 0.5 near the lower conducing plate. Here, δ(ω) is the thermal boundary layer thickness obtained using δ(ω) = 1/[2Nu(ω)]. Note that the colormap is the same for (B) and (C), but a different one is adopted for (A). In the left panel, the typical flow structures of classical thermal turbulence without any vibration are established. In the middle and right panels, with increasing ω, more and more plumes are generated and erupted due to the destabilization of thermal boundary layers. (D) Heat-transport enhancement expressed as the ratio of Nusselt numbers Nu(ω)/Nu(0) versus ω, where Nu(ω) is the Nusselt number of the vibrated convection measured at ω, and Nu(0) is the Nusselt number of classical thermal turbulence without any vibration. The horizontal dashed line marks Nu(ω)/Nu(0) = 1, and the yellow shaded area corresponds to the Nu-enhancement regime. (E) Normalized area A/LW of hot plumes as a function of ω obtained at z = δ(ω) near the lower plate. Inset shows the corresponding heat content Q/Q0 of hot plumes. Here, Q0 is the energy supplied to the system in one large-scale turnover time. (F) Sketch of the rectangular cell with the coordinate system. The horizontal vibration δcos(ωt) in the x direction is applied to the convection cell.

Figure 1 (A to C) shows the typical instantaneous flow structures observed in our simulations at Ra = 108 and Pr = 4.38. Here, the Rayleigh number and the Prandtl number are given byrespectively, where Δ is the temperature difference across the fluid layer, g is the acceleration due to gravity, α is the isobaric thermal expansion coefficient of the working fluid, v is the kinematic viscosity, and κ is the thermal diffusivity. The traditional flow structures of the RB flow are recovered at ω = 0, as shown in the left panel of Fig. 1 (A to C). Sheetlike plumes propagate horizontally near the upper and lower conducting plates, gather and cluster at the two ends of the plates, and then move upward (hot) or downward (cold) along the cell sidewalls. This forms a large-scale mean flow in the bulk and some secondary rolls at the corners.

When horizontal vibration is applied to the system, a new kind of convective instability is born. At large frequency ω, the fast horizontal oscillations of the conducting plates instantaneously induce a strong shear to the fluids near them, which makes thermal boundary layers to be markedly destabilized. As a result, a huge number of thermal plumes are generated and emitted from the destabilized boundary layers, as shown in the middle and right panels of Fig. 1 (A to C). These plumes are much coherent and energetic; they erupt at random locations over the entire boundary layers and can efficiently bring the hot or cold fluid into the convective bulk, leading to more uniform and thinner thermal boundary layers (see Discussion). Therefore, one would expect an efficient heat exchange in the system, which is the central finding of our results.

The convective heat transport through the system is usually measured by the Nusselt number, defined aswhere w is the vertical velocity, θ is the temperature, ∂ is the vertical derivative, and 〈·〉 is the average over time and any horizontal plane. Figure 1D shows the ratio of Nusselt number Nu(ω)/Nu(0) as a function of vibration frequency ω, where Nu(ω) is normalized by Nu(ω = 0) to show the enhancement effects. At low vibration frequency (ω ≾ 200, the green shaded area of Fig. 1D), the measured Nu(ω) is essentially equivalent to the RB value of Nu(0). This is because when ω is small, the induced shear is too weak to effectively perturb the boundary layers and the system is still in a state of boundary-layer–controlled thermal turbulence. As ω increases, more and more plumes are generated due to the increased shear that destabilizes the boundary layers. Because thermal plumes are the primary heat carriers that are responsible for the coherent heat transport in thermal turbulence, the marked increase of their number corresponds to the remarkable enhancement of heat transfer. As shown in Fig. 1D, at high vibration frequency (ω > 200, the yellow shaded area), Nu(ω)/Nu(0) increases rapidly with ω and a fivefold enhancement in heat flux is achieved at ω = 1700.

To quantitatively characterize this process, we use the approach in (, ) to extract hot plumes near the lower plate as the set of coordinates, where θ − 〈θ〉 > θrms and ωθ/(κΔ/H) > Nu. This criterion is based on the fact that a hot plume is a localized thermal structure being hotter than the turbulent background and carrying an excess of heat in the vertical direction. With the extracted coordinates of hot plumes, we calculate their area as A and “heat” contained in the plumes as Q = ∑cρVgridθ, where c and ρ are the specific heat and density of the working fluid, and Vgrid is the volume of each grid point. In Fig. 1E, we plot the normalized A/LW as a function of ω obtained on the horizontal slices at z = δth(ω). In contrast to the low-ω cases, in which hot plumes occupy ∼22% of the area of the lower plate, the plume area reaches to ∼39% at the highest ω = 1700 studied. In the meantime, the heat content of thermal plumes raises twofold (Fig. 1E, inset). Together, these results quantitatively illustrate that with increasing ω, the fast horizontal vibration vigorously triggers more frequent and much stronger plume emissions that gives rise to massively enhanced heat transport through the convection system.

To systematically reveal the properties of the heat-transport enhancement induced by horizontal vibration, we carry out 3D DNS in a rectangular cell spanning the Rayleigh number range 106 ≤ Ra ≤ 108 and 2D counterparts in a square box of height H spanning 107 ≤ Ra ≤ 1010. All the simulations are performed at fixed Pr = 4.38, corresponding to the working fluid of water at 40°C, and fixed δ = 0.1. The results for different values of δ can be found in the Supplementary Materials (see fig. S7). Figure 2 (A and B) shows the variations of Nu normalized by the values of corresponding RB cases (i.e., ω = 0) for both 3D and 2D simulations. The enhancement of heat transfer is rather robust; it can be found for all the values of Ra studied. In general, this enhancement is stronger for higher ω or larger Ra, and specifically, a vast increase of nearly five (six) times that without any vibration is achieved for 3D (2D) realizations. In addition, the Nu-enhancement regime shifts toward smaller ω when Ra is increased, suggesting that the enhancement of Nu occurs more easily at higher Ra due to the relatively stronger shear generated by the more turbulent flow.

Fig. 2

Heat-transport enhancement due to horizontal vibration.

(A and B) Ratio Nu(ω)/Nu(0) as a function of vibration frequency ω obtained at various Ra and fixed δ = 0.1 for (A) 3D and (B) 2D simulations. The dashed curves are the best fits of the crossover function y = log [101/4 + (ω/ω*)]4 to the respective data, from which the value of ω* is gained. The fitted values of ω* and a are given and discussed in the Supplementary Materials. (C) Ratio Nu(ω)/Nu(0) as a function of the normalized frequency ω/ω* for the same datasets as those in (A) and (B). Note that we also carried out simulations at different values of δ and plotted the ratio Nu(ω)/Nu(0) as a function of the vibrational Rayleigh number (also called Gershuni number) Ra. The same behaviors are obtained (see figs. S7 to S9).

We further note that the data shown in Fig. 2 (A and B) can be characterized by a crossover function y = log [101/4 + (ω/ω*)]4, with ω* and a being two fitted parameters. Here, ω* can be regarded as the critical vibration frequency, above which the Nu enhancement becomes notable. The dashed curves in Fig. 2 (A and B) are the best fits to the respective data, and the crossover function is seen to well describe the dependence between the normalized Nu and ω. To better compare the measured Nu(ω) at different Ra and different dimensionality, we adopt ω* to normalize the data and the results are plotted in Fig. 2C. It is seen that all data points collapse nearly on top of each other for all Ra and for both 3D and 2D cases. This signals that the Nu enhancements for all cases studied exhibit universal properties and are governed by the same vibrational convective mechanism.

Scaling of heat transport

In this section, we address the issue of how horizontal vibration modifies the constitutive scaling relation Nu ∼ Raβ of turbulent thermal convection. We show in Fig. 3 (A and B) the measured Nu versus Ra curves obtained at various ω for both 3D and 2D simulations. The traditional scaling of the RB flow in the classical regime, i.e., Nu ∼ Ra0.3 (the dashed lines in the figure), is gained at ω = 0, as expected, indicating that without any vibration the global heat transport is strongly dominated by boundary layers. As ω increases, all Nu-Ra curves display a clear power law and their scalings become steeper with increasing ω. At high vibration frequencies, the Nu-Ra data approach to the scaling of Nu ∼ Ra0.5 (the solid lines in the figure), i.e., the asymptotic ultimate scaling of heat transport in thermal turbulence. Figure 3C shows β in the best power-law fits of Nu ∼ Raβ to the respective data as a function of ω. One sees that β ≃ 0.3 at small ω and rises up with ω. At large ω, β attains a maximum value of 0.485 for three dimensions and 0.447 for two dimensions. Both of the maximum exponents are smaller than the ultimate value of 0.5, which may originate from the log-layer corrections due to the leftover boundary layers ().

Fig. 3

Heat-transport scaling of the vibrated convection.

(A and B) Log-log plots of the measured Nusselt number Nu as a function of the Rayleigh number Ra for various values of the vibration frequency ω for (A) 3D and (B) 2D simulations. From bottom to top, the symbols are △ : ω = 0; ☆ : ω = 200; ⊲ : ω = 400; ◊ : ω = 550; □ : ω = 700; ⊳ : ω = 850; ○ : ω = 1000; ▽ : ω = 1400; + : ω = 1700. The solid and dashed lines are eyeguides. (C) Fitted exponent in the power-law Nu ∼ Raβ as a function of vibration frequency ω. The solid and dashed lines mark the values of 0.5 and 0.3, respectively, for reference. Inset shows the shear Reynolds number Res as a function of ω at Ra = 108. The dash-dotted lines indicate the critical value of 200 taken by Grossmann and Lohse (, ) and that of 420 given by Landau and Lifshitz ().

The ultimate regime, first proposed by Kraichnan, has stimulated intense studies in the past two decades (, –). In this fully turbulent regime, the kinematic viscosity and thermal diffusivity of the working fluid become irrelevant and the scaling laws of heat transport can be extrapolated to arbitrarily high Rayleigh numbers. The observed asymptotic ultimate scaling of heat transfer, in the present relatively low-Ra regime of highly vibrated convection, is of great interest and fundamental importance. To physically understand this, we note that the fast horizontal vibration would induce a strong instantaneous shear to the body of fluid in the vicinity of the conducting plates. The typical length scale of such a process is the thickness of Stokes layer, . It is well known that the viscous boundary layer is much thicker than the thermal boundary layer for Pr > 1 in the unvibrational system. Whereas in thermal vibrational convection, the thickness of the Stokes layer induced by high-frequency vibration is much smaller than that of thermal boundary layer, i.e., the Stokes layer is nested deeply in the thermal boundary layer. (See also fig. S4, which displays that the strong velocity fluctuations occur inside thermal boundary layers.) Therefore, the vibration-induced shear can very efficiently disrupt thermal boundary layers and may, in turn, trigger the transition to ultimate regime.

To assess the magnitude of such a shear, we calculate the shear Reynolds number Re, based on the maximum vibration velocity δω and the Stokes layer thickness . According to Landau and Lifshitz (), the critical value for the laminar-turbulent transition of the boundary layer is Re = 420, and Grossmann and Lohse () further pointed out that such a transition can already start at a critical value of Re = 200 in turbulent RB convection, as a consequence of the experiments by Chavanne et al. (). In our present configuration, we find that the obtained Re exceeds the suggested critical values. As an example, we plot in the inset of Fig. 3C the shear Reynolds number Re as a function of ω computed at Ra = 108. It is found that Re = 196 for ω = 400 and Re = 404 for ω = 1700, which is qualitatively consistent with the fitted scaling exponents β as shown in Fig. 3C.

DISCUSSION

In the present study, we found a new kind of convective instability via the introduction of horizontal vibration to thermal turbulence. When horizontal vibration is applied to the convection system, a corresponding shear is generated for the fluids near the conducting plate. When the frequency of vibration exceeds a critical value, the induced shear becomes large enough to vigorously destabilize the boundary layers. As a consequence, massive eruptions of thermal plumes are triggered, which enables a very efficient heat exchange between the turbulent bulk flow and the boundary layers. This process overcomes the bottleneck effects imposed by boundary layers and realizes a remarkable enhancement of heat transport by up to 600%. Moreover, when the vibration-induced shear exceeds the critical value for laminar-turbulent transition, the boundary layers break down and the system starts to undergo a transition to ultimate regime with a steeper increase of Nu than the classical scaling. More generally, the concept in this work may serve as an example of how to intensify the turbulent mixing and heat and mass transfer in unstable stratification fluids by introducing translational vibration that are perpendicular to the density gradients.

For the present case of turbulent thermal vibrational convection, with increasing vibration frequency ω, more and more plumes, after detaching from thermal boundary layers, are entrained by the large-scale mean flow. This makes the mean flow more coherent and energetic, i.e., the flow moves much faster and contains the fluid that is much hotter or colder than the surrounding background (see fig. S6, A and B). Correspondingly, the flow can carry more heat. This can be quantitatively illustrated by the distributions of the local heat flux Nu (see fig. S6C), which reveals that very large values of positive Nu occur with much higher probabilities at large ω, implying the existence of more energetic events in the intensively vibrated convection. On the other hand, the vibration-induced destabilization occurs at random locations over the entire boundary layers, resulting in the eruption of thermal plumes from almost all horizontal positions. This process makes thermal boundary layers more uniform. Meanwhile, more plumes are generated by horizontal vibration, and a massive eruption of thermal plumes can efficiently bring the hot or cold fluid out of thermal boundary layers. This process highly contributes to the vertical heat flux, which, in turn, makes thermal boundary layers thinner. We find that the temperature gradient at the conducting plates becomes much steeper with increasing ω (see fig. S2). It is thus both of the above processes that minimize the roles of the boundary layers and correspondingly realize the efficient enhancement of the global heat flux through the system.

MATERIALS AND METHODS

We consider an incompressible flow governed by the coupled equations of motion for velocity field u and temperature field T (dimensionless form is θ) in the Oberbeck-Boussinesq approximation in both two and three dimensions. Horizontal vibration of amplitude A and angular frequency Ω is applied to the convection cell in the x direction. In the coordinate system associated with the cells, the corresponding acceleration thus consists of a gravitational and a vibrational part, i.e., , where and are the unit vectors in the directions of the cell length (horizontal) and height (vertical), respectively. The governing equations readwhere p is the kinematic pressure field. All quantities studied above have been made dimensionless by the cell height H, the imposed temperature difference across the cell Δ, and the free-fall velocity and time . According to these choices, the dimensionless vibration amplitude δ and frequency ω are, respectively, given by

The dynamics of the system are then controlled by the Rayleigh number Ra, the Prandtl number Pr, the dimensionless vibration frequency ω, and the oscillatory Rayleigh number Raos. Here, the oscillatory Rayleigh number Raos is defined as

A combination of these parameters, named the Gershuni number Ra and given byis frequently used as a reasonable proxy of the ratio between the mean vibrational buoyancy force and the product of momentum and thermal diffusivities (, ).

The 3D simulations are carried out over the Rayleigh number range 106 ≤ Ra ≤ 108 and in a rectangular cell of height H, length L, width W, and aspect ratio H : L : W = 1 : 1 : 0.3. For 2D, the convection cell adopted is a square box of height H and length L, and the simulations cover 107 ≤ Ra ≤ 1010. For all cases, the Prandtl number is fixed at Pr = 4.38, corresponding to the working fluid of water at 40°C. No-slip boundary conditions for the fluid are applied at solid walls. The sidewalls are thermally insulated, and the upper and lower conducting plates are held at constant dimensionless temperatures θtop = −0.5 and θbot = 0.5, respectively.

Our simulations are conducted by using the open source spectral element code Nek5000. This code has been well validated in simulating turbulent RB convection (, ). In the present study, the number of spectral elements (N × N × N) increases with both Ra and ω. For example, it is increased from 64 × 20 × 72 to 128 × 40 × 160 as ω increases from 0 to 1700 at Ra = 108. The order of the Legendre polynomials is chosen to be 7 on each spectral element. To well resolve the boundary layers, the computational meshes are refined close to solid surfaces. For all runs, thermal boundary layers are resolved with at least 25 meshes, resulting from primary and secondary nodes, and the ratio of the maximum mesh width to the Batchelor length scale remains less than unity, implying that the smallest length scales are adequately resolved in our simulations. We have also used a finite-difference code to cross validate the 3D results obtained at Ra = 108. The code has been described in detail elsewhere (), and thus, we give only its main features here. The spatial derivative terms are approximated by a second-order central difference scheme with nonuniform staggered grids. A fractional step approach is used to solve momentum equations and a multi-grid strategy is implemented to accelerate the iteration process in the Poisson equation. The number of grid points increases from 512 × 256 × 512 at ω = 0 to 1024 × 384 × 1280 at ω = 1700.

For the temporal resolution, the time step Δ is chosen to fulfill the Courant-Friedrichs-Lewy (CFL) conditions at low ω (≤200), i.e., the CFL number is 0.2 or less, while Δ ≤ 2π/90ω at high ω (>200), i.e., at least 90 time steps are calculated during each vibration period. We note that Δ/τη < 0.01 for all runs with the Kolmogorov time scale τη, thus guaranteeing the adequate temporal resolution. All statistics are calculated over an averaging time of more than 500 dimensionless time units after the system has reached the steady state. The time convergence for the measured Nu is checked by comparing the time averages over the first and the last halves of the simulations (), and the resulting convergence is less than 2% for all the simulations. In addition, the discrepancy between the Nek5000 and finite-difference codes is less than 3% for all ω.