The effects of viscoelasticity, here caused by polymer additives, on Rayleigh Bénard convection flows are investigated via direct numerical simulations at a marginally turbulent Rayleigh number. Simulations with a range of polymer length and relaxation time scales show heat transfer enhancement (HTE) and reduction (HTR). The selection of HTE and HTR depends strongly on the maximum extensional viscosity of the solution, whereas the magnitude of heat transfer modification is a function of both the maximum extensional viscosity and relaxation time of the polymer solution. The underlying physics of HTE and HTR are explored, and a mechanism of the interaction between convection cells and polymers is proposed. The findings are extrapolated to high *Ra* to shed some new light onto experimental observations of HTR.

## I. INTRODUCTION

In the turbulence research community, high-molecular weight polymer additives are well known for their ability to reduce turbulent drag.^{1} Their effect on turbulent heat transport is far less documented, even though polymer solutions may be of interest in the management of heat transfer, for instance. Recent experimental studies^{2–5} demonstrated that dilute solutions of polymers may cause heat transfer reduction (HTR) or enhancement (HTE) in turbulent natural convection flow. The flow of interest here, specifically the Rayleigh–Bénard convection (RBC) flow, is generated by a vertical negative gradient of temperature between two horizontal plates. When the Rayleigh number *Ra*, ratio of buoyancy forces to viscous forces, is large enough, convection cells emerge and drive convection heat transfer. These convection cells possess two critical flow patterns. First, convection cells create boundary layer flows on the top and bottom walls. In the absence of side walls, e.g., as in simulations of infinite horizontal space, and at low *Ra*, the length of the boundary layers is determined by the distance between the upward and downward plumes, the second critical flow pattern, of convection cells. Plumes drive heat transfer through the vertical flux of temperature *wθ*, where *w* and *θ* are the vertical velocity component and the temperature, respectively.

As discussed in the original HTR experiment,^{2} the effects of polymers on a convection flow may not seem obvious from the knowledge of polymer drag reduction. A key property of a polymer flow is its relaxation time *λ*, which defines the time necessary for fully stretched polymers to return to the coiled configuration, minimum energy configuration of polymers. Lumley’s theory^{6} proposed that, for drag reduction to occur, the polymer relaxation time needs to be equal to or larger than the smallest time scale of the flow *λ*_{K} (Kolmogorov time scale), thus defining the critical Weissenberg number, ratio of relaxation time to flow time scale, as *Wi*_{c}4 = *λ*/*λ*_{K} ≥ 1. The typical relaxation time of drag-reducing polymer solutions may range from 1 ms to 1 s.^{7,8} Ahlers and Nikolaenko^{2} pointed out that the typical relaxation of the polyethylene oxide (PEO) solution they used (less than a second) is smaller than the typical time scale of RBC eddies and velocity fluctuations, and thus, the authors stated that their flow should not be affected by polymers according to Lumley’s criterion. The authors also remarked that the boundary layers in RBC are laminar,^{9} and polymers are known to reduce drag only in turbulent wall flows. They speculated that the heat transfer reduction may result from an increased stability of the boundary layers, which, in turn, would result in lesser plume emissions and HTR.

Cheng *et al*.^{10} showed in two-dimensional RBC simulations that the amount of HTR behaves non-monotonically with the Weissenberg number. They related this non-monotone behavior to modifications of the period of the large scale circulation and boundary layer thickness. Adding to the confusion, Benzi *et al*.^{11} simulated heat transfer enhancement (HTE) and reduction (HTR) in both the direct numerical simulation (DNS) and shell model simulation of homogeneous turbulent convection. The most puzzling feature is certainly HTE, obtained for *Wi* ≲ 1, which has been recently observed experimentally.^{3,4} HTE was also observed by Cheng *et al*.^{12,13} in simulations of a single plume in thermal convection for a very narrow range of the parameter space. They also highlighted the dependence of heat transfer modification (HTM) on the polymer solution elasticity.

The present article investigates the impact of polymers in marginally turbulent RBC flows. Direct numerical simulations using a common constitutive model for polymer flows attempt to capture the effects of viscoelasticity on the dynamics of natural convection flows with the highest fidelity allowed by available computing power. Consequently, the present study is limited to a low Rayleigh number on the spectrum of turbulent natural convection, and extrapolation of our results and conclusions to higher Rayleigh numbers of the order of the HTR and HTE experiments^{2–5} is speculative. Nonetheless, the detailed identification of interactions between polymers and natural convection cells at a low Rayleigh number allows for new insights into the mechanism driving HTM.

After a short introduction to the governing equations and numerical methods, the effects of polymers on two-dimensional convection cells are briefly discussed to highlight some of the key features of polymer/convection cell interactions in a simpler system than in the 3D system. This paper then focuses on three-dimensional simulations in a parameter space where the polymer length and relaxation time are varied. For all cases, heat transfer modification (HTM) is quantified as

Here, the Nusselt number *Nu* is the ratio of the convective to conductive heat fluxes and the subscripts *p* and *n* denote polymeric and Newtonian flows, respectively. The HTE and HTR regimes are then investigated using the two simulations, which achieve maximum HTE and HTR within our parameter space. This study focuses on the changes in convection cells by the addition of polymers and the energy transfers between polymers and flows and draws some conclusions on the mechanisms of HTE and HTR. This paper concludes with a discussion of HTR in experiments.

## II. METHODS

### A. Governing equations and parameters

The natural convection flows under investigation are simulated in a Cartesian domain defined by the orthonormal vector base (**e**_{x}, **e**_{y}, **e**_{z}), where *x*, *y*, and *z* are the two horizontal directions and the vertical direction, respectively. The components of the velocity vector **u** are *u*, *v*, and *w* and are normalized by the free-fall velocity or convection velocity $Uc=\alpha gH\Delta $, where *α*, *g*, *H*, and Δ are the fluid’s volumetric expansion coefficient, gravity, domain’s height bounded by two horizontal and isothermal walls, and temperature difference between these two walls. The upper wall is cold and the bottom is hot to drive natural convection.

The Rayleigh number, ratio of buoyancy forces to thermal and momentum diffusive forces, *Ra* = *αgH*^{3}Δ/(*νκ*) is low, *Ra* = 10^{5}, on the lower end of the turbulent convection regime.^{14} The kinematic viscosity and thermal diffusivity of the solvent are identified as *ν* and *κ*, respectively. The Prandtl number is also fixed, *Pr* = *ν*/*κ* = 7, which is realistic for water. The Reynolds number based on *U*_{c} and *H* is therefore *Re* = *Ra*^{1/2}*Pr*^{−1/2} = 119. The flow is incompressible (** ∇** ·

**u**= 0), and the buoyancy effect is modeled by the Boussinesq approximation. The direct numerical simulations (DNS) discussed in this paper solve the following transport equations for velocity and temperature

*θ*:

where *p* is the pressure. The non-dimensional formalism employed here is common for simulations of natural convection.^{15–17}

The Nusselt number is defined as

where $a\u2032=a\u2212a\xaf$ is the fluctuation of *a* with respect to the average over time and horizontal directions $a\xaf$ and ⟨*a*⟩_{V} designates averaging over time and over the entire flow volume.

The parameter *β* is the ratio of solvent viscosity to the zero-shear viscosity of the polymer solution and affects both the viscous stress and polymer stress terms in Eq. (2). The polymer stress tensor **T** is computed using the FENE-P^{18} model,

where the tensor **C** is the local conformation tensor of the polymer solution (defined below) and **I** is the unit tensor. The maximum elongational viscosity,

is the increase in elongational viscosity under infinite uniaxial elongation rate. With the parameters chosen for our simulations (25 ≤ *L* ≤ 100 and *β* = 0.9), the maximum possible increase in elongational viscosity is three orders of magnitude. The properties of the polymer solution are *β*, the relaxation time *λ*, here based on the convection scales (*Wi* = *λU*_{c}/*H*), and the maximum polymer extension *L*. The FENE-P model assumes that polymers may be represented by a pair of beads connected by a non-linear spring and defined by the end-to-end vector **q**. The conformation tensor is the phase-average of the tensorial product of the end-to-end vector **q** with itself, **C** = ⟨**q** ⊗ **q**⟩, whose transport equation is

On the rhs of Eq. (7), the first two terms are responsible for polymer stretching by hydrodynamic forces, while the third term models the restoring entropic force that tends to bring stretched polymers to their least energetic state (coiled).

### B. Numerical methods

Equations (2)–(7) are solved using finite differences on a staggered grid, following the earlier description of our numerical code.^{19} The momentum equation uses non-dissipative, energy-conserving central schemes, whereas the advection term in the conformation tensor equation requires an upwind compact scheme with low small-scale dissipation supplemented by a local artificial diffusion (LAD).

The computational domain is bounded by two horizontal walls located at *z* = ±*H*/2. The walls are no-slip, no-through (*u* = *v* = *w* = 0), and isothermal [*θ*(−*H*/2) = 1, *θ*(+*H*/2) = 0]. Periodicity is imposed in the horizontal directions (*x*, *y*). At the walls, the boundary condition for *C*_{ij} is given by Eq. (7) without the advection term.

The code has been used to elucidate the mechanism of polymer drag reduction^{20,21} and for the discovery of Elasto-Inertial Turbulence (EIT).^{22–25} Simulations of EIT have also been reproduced by a spectral code,^{26} demonstrating that finite differences with carefully chosen numerical methods are perfectly adequate. The modification of the code for the simulation of natural convection has been validated against existing databases of natural convection simulations^{16} (not shown here).

### C. Resolution and parameter space

When temperature operates as a passive scalar in a turbulent flow, large Prandtl numbers cause the smallest scale of the turbulent temperature field, *η*_{B}, to be smaller than the smallest scale of turbulence, the Kolmogorov scale^{27} *η*_{K}. The Kolmogorov scale^{27} is predicted by the rate of dissipation of kinetic energy *ε* at which turbulent energy cascades down the scales to the viscous scale $\eta K=(\nu /\epsilon 3)1/4$. For large Prandtl numbers,^{28} $\eta B=\eta K/Pr$ is the Batchelor scale. The same theory obviously applies to the transport of any high-Schmidt number passive scalars. This discussion is relevant to the system of interest here, since both temperature and polymers have Prandtl and Schmidt numbers larger than unity. To the knowledge of the authors, it is not clear if the Batchelor scale directly applies in the case of active scalar, i.e., for two-way coupled scalar-flow dynamics. In the case of polymer flows, we^{25} established that an accurate representation of polymer dynamics requires the resolution of scales smaller than the Kolmogorov scales, but still not as small as the Batchelor scale predicted from the Schmidt number of polymers in water *Sc* = *O*(10^{6}). The polymer dynamics, more so than the temperature dynamics, poses a major constraint on the grid resolution. The balance between the resolution requirement and the need to probe a large parameter space of polymer solutions motivates the relatively low Rayleigh number of our study.

The adopted mesh size allows for the resolution of turbulence smallest scales, which was assessed by checking the spectral decay of kinetic energy and temperature variance, as well as doubling the resolution. The horizontal dimensions of the domain were doubled for the Newtonian simulation without any effect on the velocity and heat flux statistics.

All simulations use a smallest grid step in the wall normal direction of 2.5 × 10^{−5}*H*, as it was determined to be necessary for viscoelastic simulations. The present investigation generated two sets of production run. The first was dedicated to the mapping of the Nusselt number as a function of *L*, the maximum polymer extension, and the Weissenberg number *Wi*. This set of simulations was performed in a computational domain of dimensions 8*H* × 8*H* × *H* and using 128 × 128 × 129 mesh nodes. The second set is used to discuss in more details three cases: Newtonian and two significant representations of HTE and HTR in a larger domain, 16*H* × 16*H* × *H* (256 × 256 × 129 nodes). The larger domain is only needed for the Newtonian flow, whose horizontal integral scale is much larger than the polymeric RBC flows. Polymeric simulations were repeated in a 8*H* × 8*H* × *H* with 256 × 256 × 257 to ensure grid convergence. No significant change in heat transfer was observed between the three different sets of simulations.

As stated earlier, the Rayleigh and Prandtl numbers are fixed: *Ra* = 10^{5} and *Pr* = 7. For polymer parameters, the viscosity ratio *β* is kept constant, *β* = 0.9, matching previous viscoelastic channel flow simulations performed in previous studies.^{19,20,22} The parameter space under investigation is a matrix (*L*, *Wi*) with *L* = (10, 25, 50, 100) and *Wi* = (1, 5, 10, 15, 20, 25, 30, 35, 40, 45), with a few additional simulations for *L* = 15, 20, 30, 40 at *Wi* = 10 and *L* = 25, 40 at *Wi* = 40.

### D. Data processing

The reference simulation is the Newtonian case. To achieve rapid destabilization of the flow, the initial Newtonian flow field consists of a row of three-dimensional Taylor–Green vortices of diameter *H*/2 at the center of the channel. White noise with an rms of 10% of the buoyancy velocity is superimposed on the initial field.

Convergence to the statistical steady state is monitored through time series of the instantaneous Nusselt number and volume averaged turbulent kinetic energy (TKE). For production runs, statistical convergence is verified using the balance of volume averaged TKE, which will be derived below. Statistics presented here are averaged over 400 flow fields sampled at a non-dimensional frequency of 1. Variations of statistical quantities in the vertical direction are plotted over one half of the domain, leveraging the symmetry or antisymmetry of each quantity, as a function of the distance from the bottom wall,

Polymer flow simulations at a given *L* were started from a fully turbulent Newtonian flow field with the conformation tensor equal to unity at the lowest *Wi* = 1. The same protocol as described above is used to monitor statistical convergence. The last sample field in the simulation run at a given *Wi* is used as the initial condition for the next *Wi* run.

## III. EFFECTS OF POLYMER ADDITIVES ON 2D CONVECTION CELLS

The geometrical complexity of 3D convection cells hinders the investigation of the physics at play in plumes or boundary layers. For the same *Ra* and *Pr* as for the production runs, we turn to two-dimensional simulations to observe the interplay between polymers and convection cells. A distinct advantage of 2D natural convection flows between two horizontal plates (with the periodic condition in the horizontal direction) is the possibility to achieve quasi-stable flows. A similar reductionist investigation was recently performed in a single plume for a FENE-P fluid as well.^{29} In a computational domain with an aspect ratio of 2:1, flow systems with two pairs of convection cells may experience long durations of stability or, otherwise, may be quasi-perfect oscillators. Although the study of the stability of 2D natural convection flow is not the focus of this work, we performed a few simulations in a small domain (2*H* × *H*) with a resolution of 128 × 129 for various *L* and *Wi*.

This section discusses two simulations, one at low elasticity and the other at high elasticity. The former experiences a small HTE of the order of 5%, whereas the latter records an HTR of the order of 25%. As will be shown below, the 2D simulations discussed here are representative of the states observed in 3D.

For a Newtonian flow [Figs. 1(a) and 1(b)], the velocity profile of plumes exhibits the bell shape of a jet and the velocity distribution in boundary layers resembles that of a wall jet. The stagnation point of an impinging plume is naturally a high-pressure region and the region where the plume lifts off from the wall is one of low-pressure (not shown). Consequently, boundary layers are subjected to a non-zero pressure gradient and a mixture of shear and extensional strain. An interesting feature of this particular simulation is its quasi-periodic behavior, where the size of the two convection cells fluctuates resulting in one larger plume and one smaller plume, as seen in Fig. 1(a).

From the perspective of polymers, plumes and boundary layers create the most significant dynamics. In polymer drag reduction in 3D channel flows, polymers stretch in bi-axial extensional flows.^{21} Figure 1(b) shows the contours of the inverse of the kinematic vorticity number,^{30}

where *ω* = *∂w*/*∂x* − *∂u*/*∂z* is the vorticity and *S*_{ij} = (*∂u*_{i}/*∂x*_{j} + *∂u*_{j}/*∂x*_{i})/2 is the strain rate tensor. The kinematic vorticity number relates to the second invariant of the velocity gradient tensor,

which is commonly used to identify vortices (the *Q*-criterion^{31}). Here, $Nk\u22121$ is smaller than unity in rotation dominated flows and larger in extensional–compressional flows. The plumes and the upper edge of the boundary layers are regions where extension–compression dominates [Fig. 1(b)]. In the bulk of the boundary layers, $Nk\u22121$ is around unity, indicating mostly shear flow.

Figure 1(c) shows the effects of the polymer solution with low maximum extensional viscosity, *E* = 62.5, on the temperature field and velocity profiles in plumes and boundary layers. Whereas plumes are still regions of maximum or minimum temperature, the velocity profile across a plume exhibits two maxima surrounding a minimum at the centerline. The flow is quasi-steady, thanks to two convection cells, and thus plumes, of the same size.

The centerline of the plume is also where the maximum first principal normal stress difference occurs, as shown in Fig. 1(d). The first principal normal stress difference is defined as

where the parallel and orthogonal stresses are

with $T\u03031$ being the eigenvalue of the first principal axis of tensor **T** and components $T\u03032$ and $T\u03033$ being the eigenvalues with eigenvectors perpendicular to the first eigenvector. The first principal normal stress difference is very similar to the first normal stress difference at the origin of the famous rod climbing effect,^{18} or the Weissenberg effect, which causes a free-surface polymeric solution to climb a rotating rod. The shear generated by the rotation of the rod stretches polymers and the stretched polymers produce normal stress, or hoop stress, in the direction perpendicular to their stretching. The absence of the solid boundary at the free surface enables the normal stress to push the fluid upward. The definition of *N*_{p} is rendered necessary by the nature of flow generated by convection cells, which is a mixture of extensional, shear, and curved-streamline flows. In boundary layers, *T*_{xx} dominates while *T*_{zz} takes the lead in plumes. More importantly, *N*_{p} is also linked to the definition of extensional viscosity.^{32} The use of *N*_{p} is therefore motivated by its relation to two of the most critical physical phenomena of polymer flows: extensional viscosity and Weissenberg effect.

For low *E* [Fig. 1(d)], large values of the first principal normal stress difference are found between the two convection cells, in the first half of the plume from its inception where the flow is mostly extensional. As mentioned earlier, the location of velocity deficit in the plume coincides with the maximum of the first principal normal stress difference, which is also the location of the maximum polymer stretch. The normal stress is much smaller in boundary layers than in plumes. Large regions of *N*_{p} appear to be mostly confined to the extensional flow region close to the stagnation point, where the incoming plume impinges the wall. The velocity deficit in between convection cells may be caused by the increase in elongational viscosity and/or the normal stress exerted on the convection cell.

The high maximum elongational viscosity flow (*E* = 3 × 10^{4}) confirms the observation made with low-*E* fluid but in a more dramatic manner. The velocity deficit is now between a half to a third of the maximum upward velocity [Fig. 1(e)]. The thermal plumes are wider and more diffused horizontally than for the low-*E* fluid. Regions of large first principal normal stress difference [Fig. 1(f)] are extending higher in the plumes both vertically and horizontally. There are some clear signs of EIT, which will not be discussed here. Boundary layers also experience thin sheets of large first principal normal stress difference.

In wall-turbulent flows undergoing drag reduction, the regions of high polymer stretch are found in upward and downward flows caused by near-wall vortices and are the key mechanism of polymer drag reduction.^{20,21,33} A similar mechanism appears to be at play for both convection cells and quasi-streamwise vortices in the near-wall turbulent channel or boundary layer flow. For the flow with low maximum elongational viscosity [Fig. 1(d)], regions of large *N*_{p} are concentrated in the first half of plumes. The obvious correlation between the presence of large *N*_{p} and the observed velocity deficit on the plume’s centerline could be found in the Weissenberg effect described above. The first principal normal stress difference in plumes pushes adjacent convection cells apart, thereby creating the centerline velocity deficit. The presence of significant first principal normal stress difference in boundary layers appears to be a strong function of the maximum extensional viscosity of the solution, which is consistent with our observation in channel flows.

Another interesting observation is that the low maximum elongational viscosity solution is quasi-stable, whereas the other solution is chaotic, as can be inferred from the structure of *N*_{p} in plumes. Although not the objective of the present study, the presence of EIT in this 2D flow justifies the use of a local artificial diffusion method. Using a low-Schmidt number global diffusion approach should destroy EIT in a similar manner that was established for channel flow.^{25}

The conclusion of this brief 2D study is that either the polymer first principal normal stress difference mostly acts in plumes by pushing convection cells apart, resulting in a velocity deficit at the center of the plume, or the extensional viscosity causes the splitting of the cell, or a combination of both. This observation is consistent with the results of the investigation of a single buoyant plume modified by polymers.^{29} The most striking difference between HTE and HTR takes place in boundary layers. In the HTE flow, the first principal normal stress is visibly much less in boundary layers than in plumes. In the HTR flows, plumes and boundary layers appear to produce a similar magnitude of *N*_{p}.

## IV. EFFECTS OF POLYMER ADDITIVES ON 3D CONVECTION CELLS

### A. Heat transfer modification

HTE is observed at *Wi* ≤ 5 for nearly all polymer lengths except for *L* = 100, yet the present simulated enhancement (HTM = + 12% for *L* = 25 and *Wi* = 10) is much more modest than the maximum observed in the shell model of homogeneous convection,^{11} *Nu*_{p}/*Nu*_{n} ≈ 6, or in the DNS by the same authors where *Nu*_{p}/*Nu*_{n} ≈ 2.4 with *L* = 30. The absence of walls and possibly a lower Prandtl number, *Pr* ≤ 1, might therefore magnify the HTE ability of viscoelastic thermal convection flows. Figure 2(a) also highlights the critical role of the polymer length on the heat transfer performance. Long polymers, *L* ≳ 50 or maximum elongational viscosity larger than 250, show hardly any HTE. For the range of *Wi* considered, the short polymer simulations experience a decrease in the HTE effect at high *Wi* yet do not reach HTR. The two HTR flows exhibit a power-law behavior at high *Wi*, with *Nu* ∝ *Wi*^{−0.1} and *Wi*^{−0.2}, for *L* = 50 and 100, respectively. Although this result indicates a dependence of the power law exponent in *L*, the limited available data do not allow for a definitive conclusion. Combining the shell model and the scaling theory of Grossmann and L̂hose,^{34} Benzi *et al*.^{11} predicted the HTR behavior as *Nu* ∝ *Wi*^{−1}. Since the shell model predictions are performed at much higher Rayleigh numbers, it may be assumed that Benzi’s work may describe an asymptotic behavior of HTR when *Ra* → *∞*.

The only conclusion that can be drawn from our data is that the HTR behavior as a function of *Wi* is strongly dependent on *L*. Figure 2(b) provides some insight into the evolution of HTM as a function of *L* for two different Weissenberg numbers, *Wi* = 10 and 40. For both *Wi*, HTM first increases to the same maximum *L* ≈ 20 of HTE as *L* grows. The critical *L* at which the flow transitions from HTE to HTR seems to decrease with the increase in *Wi*. Similarly, Fig. 2(a) suggests that the critical *Wi* decreases with the increase in *L*. The current discretization of the parameter space does not allow for the derivation of a correlation for the transition boundary between HTE and HT in the (*Wi*–*L*) space. This will be investigated in the future.

### B. Velocity and temperature statistics

We now focus on two simulations representative of HTE (*L* = 25, *Wi* = 10, and HTM = + 11%) and HTR (*L* = 100, *Wi* = 40, and HTM = −31%). The maximum elongational viscosities [*E*, see Eq. (6)] are 62.5 for HTE and 1000 for HTR. Figure 3(a) shows the horizontal and vertical velocity rms as a function of the vertical distance from the wall in the lower half of the flow domain. The addition of polymers considerably reduces the thickness of the boundary layer, estimated from the maximum of *u*_{rms},^{16} as well as the intensity of horizontal velocity fluctuations. Despite the reduction of horizontal velocity fluctuations, the intensity of vertical velocity fluctuations at HTE remains very close to that of the Newtonian flow. On the other hand, the HTR flow shows significant reduction of velocity fluctuations in both directions. Interestingly, for both polymeric flows, the maximum of *w*_{rms} is roughly 1.5 times the maximum of *u*_{rms}, whereas the maxima of *u*_{rms} and *w*_{rms} are very close for the Newtonian flow.

Convection cells in HTR flow are clearly slowed down compared to Newtonian flow, but the HTE flow presents a more puzzling picture: how could convection cells maintain the same level of rms of velocity fluctuations in plumes as in Newtonian flow when the rms in boundary layers is nearly halved? For both flows, polymers affect heat transfer throughout the domain [Fig. 3(b)]. Adding to the confusion, temperature fluctuations in HTE are identical to the Newtonian flow in the wall region and show a slight decrease in the core, whereas, in HTR, these fluctuations experience a reduction in the near-wall region and an increase in the core.

### C. Topology and dimensions of convection cells in polymer flows

Isosurfaces of positive and negative temperature fluctuations in the entire domain (Fig. 4) show the dramatic impact of polymers on convection cells. In the Newtonian flow [Fig. 4(a)], the hot fluid flows upward at the center of the convection cells and cold fluid flows on the outside of the convection cell forming ridges. This pattern known as the l-hexagonal convection cells at lower *Ra* is to be expected for a liquid.^{14} At the present Rayleigh number, the organization has two scales. The first one from the center of the upward hot flow to the ridges of downward cold flow is roughly 3 − 4*H*. There is an apparent fragmentation within the large convection cells at a smaller scale of the order 1*H*.

This complex structure, composed of at least two scales, disappears in polymer flows for both HTE [Fig. 4(b)] and HTR [Fig. 4(c)]. The topology of cells shows a higher degree of organization, akin to the type of organization observed in lower *Ra* Newtonian flows, with only one observable scale, compared with that of the convection cells. Apart from reduction in upward and downward flow velocities, HTE and HTR flows are qualitatively similar.

Horizontal correlation of temperature fluctuations (Fig. 5) confirms the reduction in convection cell horizontal dimension for both HTE and HTR. Here, the correlation is defined as

where $r=x2+y2$, *r*_{∥} denotes an arbitrary horizontal distance, and ⟨⋅⟩_{r,t} denotes the averaging over the entire horizontal plane and time.

The significant increase in the number of convection cells observed in Figs. 4(b) and 4(c) for polymeric flows is quantified by Fig. 5. The location of the (negative) minimum of the temperature correlation is a measure of the convection cell horizontal length scale (anti-correlation between upward and downward plumes). Compared to the Newtonian flow, polymers reduce this length by more than a third in both HTE and HTR, which confirms the observations derived from Fig. 4. The convection cells in HTR are slightly larger than those in HTE (HTE ∼0.75*H*, HTR ∼0.81*H*, and Newtonian ∼3.1*H*). The positive extremum following the negative minimum is also more pronounced at HTR, indicating that the distribution of sizes of convection cells is narrower at HTR than at HTE.

With the knowledge that the horizontal length of convection cells is profoundly reduced in polymeric flows, we revisit the question that arose from the comparison of velocity fluctuations in plumes and boundary layers between the flows under consideration [Fig. 3(a)]. For the Newtonian flow, the peak rms of vertical and horizontal velocity fluctuations are almost identical, which was observed both numerically^{35} and experimentally^{36} at similar Rayleigh numbers. In HTE flow or low-*E* solutions, the peak rms of the horizontal velocity fluctuations for the simulation discussed here is about half that of the vertical fluctuations. The reduction is likely to result from the dramatic reduction in their development length or the horizontal length scale of the convection cells. The similarity in magnitude of the rms of vertical velocity fluctuations suggests that polymer stress does not impact the buoyancy force in HTE.

For high-*E* solutions, the horizontal length scale of convection cells is roughly the same size as that for HTE, yet the rms of horizontal and vertical velocity fluctuations in Fig. 3(a) is half that of HTE. The increase in the number of convection cells is now compensated by a reduction in momentum in both boundary layers and plumes. The HTR and HTE convection cells are visually similar, and only their boundary layer and plume velocities are much smaller in HTR than in HTE.

As discussed by Benzi *et al*.,^{37} the statistics suggest, at least for HTE, that the main polymer effect is in the boundary layer. However, our observation on the effects of the first principal normal stress difference on 2D convection cells forces us to consider other scenarios. At such a low Rayleigh number, it could be assumed that the action of *N*_{p} on the dynamics and/or structure of plumes may modify the dynamics of boundary layers, even if the effects of polymers in that region may be negligible. At this point, the only certainty remains that boundary layers become buoyant over a much shorter distance, which is a likely consequence of the slower velocity in this region of the convection cells. However, the data discussed so far does not allow us to understand why the convection cells are modified and how these modifications lead to HTE and HTR.

### D. Relation between 2D and 3D simulations

Figure 6 compares random snapshots of temperature contours for the three flows of interest in a vertical slice of the domain. This figure illustrates the reduction in the horizontal length of the convection cells found by the two-point correlation analysis of temperature fluctuations (Fig. 5). Plumes in the HTE flow appear thinner and straighter than in the HTR flow. As observed previously in the 2D simulation with long polymers [Fig. 1(e)], plumes in the HTR flow fan out more as they approach a wall.

The difference between HTE and HTR is more visible in the contours of the first principal normal stress difference shown in Fig. 7. As for the 2D simulation with short polymers [Fig. 1(d)], the HTE flow [Fig. 7(a)] does not produce significant first principal normal stress difference in the boundary layers of convection cells and regions of large normal stress are confined to the first half of plumes. In Fig. 7(b), the HTR flow shows similar patterns to the 2D simulation with long polymers [Fig. 1(f)]: the regions of large *N*_{p} extend further into the plumes and display a more chaotic behavior. Boundary layers experience levels of first principal normal stress difference comparable to the levels of *N*_{p} in plumes.

In conclusion, the similarity of first principal normal stress difference dynamics in 2D [Figs. 1(d) and 1(f)] and 3D (Figs. 6 and 7) suggests that the difference between HTE and HTR is likely governed by the dynamics of the first principal normal stress difference in the boundary layers of convection cells, since this is the only visual difference between the two flows.

## V. MECHANISMS OF HTE AND HTR

### A. Energy balance

The interactions between the dynamics of polymers and Bénard cells are investigated through the contributions of the term

to the heat flux $w\theta \u2032\xaf=\kappa \Delta /H(Nu\u22121)$. In Eq. (13), $\epsilon p\xaf$ can be interpreted as the combination of a transport and energy transfer term between polymers and flow^{24} or as the work or flux of the polymer “force” in Eq. (2), **f**_{p} = ((1 − *β*)/*Re*)** ∇T**. By integrating the transport equation of turbulent kinetic energy [derived from Eq. (2)] over the entire volume, one may derive an equation for the Nusselt number,

^{38}which in the case of polymeric flows is given by

where *ε*^{u} is the dissipation rate of turbulent kinetic energy. The distinction between the horizontal (subscript *bl* for the boundary layer) and vertical (*pl* for plumes) contributions of **u** · ($\u2207$ · **T**) to *ε*^{p} [see Eq. (13)] is inspired from the study of a unifying theory of turbulent RBC.^{34} These authors derived a decomposition from the volume integration over the boundary layer (stress carried by horizontal momentum) and plume (stress carried by vertical momentum) regions. Here, the proposed decomposition of $\u27e8\epsilon p\u27e9V$ is based on the horizontal and vertical contributions of **u** · ($\u2207$·**T**), i.e., (*u***e**_{x} + *v***e**_{y}) · ($\u2207$·**T**) and *w***e**_{z} · ($\u2207$ · **T**), respectively. This approach is motivated by the high level of organization of convection cells and the correlation between strong polymer-flow interaction (via first principal normal stress difference) and plumes and boundary layers.

Figure 8 shows the contribution of $\u2212\epsilon blp$ and $\u2212\epsilon plp$ to Eq. (14) in a semi-logarithmic scale. As observed from the contours of first principal normal stress difference in 2D and 3D (Figs. 1 and 7), the largest positive contribution of −*ε*^{p} comes from the plumes at HTE but from the boundary layers at HTR. Yet, both HTE and HTR experience negative contributions of $\u2212\epsilon blp$ in the very near-wall region. Above the near-wall, $\u2212\epsilon blp$ becomes positive before vanishing in the core, whereas the contribution of $\u2212\epsilon plp$ is first negative and then becomes positive in the core.

For the Newtonian flow, the volume average vertical heat flux is entirely balanced by *ε*^{u} [Eq. (14)]. Table I breaks down the contribution of the three terms involved in Eq. (14) for the two HTE and HTR flows under consideration. The Newtonian dissipation rate of TKE dominates the energy balance in the HTE case with 72% of the whole LHS of Eq. (14). The second contribution is the elastic energy in the plume (21%), three times that of the contribution in the boundary layers (7%). The energy balance is reversed in HTR flows with 64% from elastic energy: 38% from boundary layers and 26% from plumes.

Term . | Region . | HTE (%) . | HTR (%) . |
---|---|---|---|

$\u27e8\epsilon u\u27e9V/\u27e8w\theta \u2032\u27e9V$ | Whole domain | 72 | 36 |

$\u2212\u27e8\epsilon plp\u27e9V/\u27e8w\theta \u2032\u27e9V$ | Whole domain | 21 | 26 |

$\u2212\u27e8\epsilon blp\u27e9V/\u27e8w\theta \u2032\u27e9V$ | Whole domain | 7 | 38 |

Term . | Region . | HTE (%) . | HTR (%) . |
---|---|---|---|

$\u27e8\epsilon u\u27e9V/\u27e8w\theta \u2032\u27e9V$ | Whole domain | 72 | 36 |

$\u2212\u27e8\epsilon plp\u27e9V/\u27e8w\theta \u2032\u27e9V$ | Whole domain | 21 | 26 |

$\u2212\u27e8\epsilon blp\u27e9V/\u27e8w\theta \u2032\u27e9V$ | Whole domain | 7 | 38 |

### B. Work of polymer force

To gain a deeper understanding of the interactions between convection cells and polymer dynamics, we analyze the sign of $\epsilon p\xaf$, which also corresponds to the work of the polymer force [Eq. (13)]. In other words, when $\epsilon p\xaf$ is positive (negative), the polymer force works with (against) the momentum. As shown by Fig. 8, the polymer force exerts a favorable work on the momentum in the near-wall region, whereas the rest of the flow experiences an adverse work of the polymer force. The position at which $\epsilon p\xafbl$ ($\epsilon p\xafpl$) changes sign is shifted downward (upward) for HTR compared to HTE.

The relation between normal stress and the regions of positive and negative $\epsilon p\xaf$ is now analyzed in terms of the intensity of the first principal normal stress difference *N*_{p}. As shown in Fig. 7, the regions of high *N*_{p} are sheet-like. The intensity of these high-*N*_{p} sheets is further quantified based on conditional statistics, as commonly done for intermittent systems. Figure 9 shows the vertical profiles of the mean $N\xafp$ and conditional mean of *N*_{p}, the latter being defined as

where the sum is applied over all nodes of all sample fields for each horizontal plane. For both HTE and HTR cases, Fig. 9 surveys a range of threshold values *N*_{p,th} = [25, 50, 75, 100, 125, 150, 175, 200]. The range is centered around the median of the colormap used for the *N*_{p} contours in Fig. 7. The chosen range of thresholds allows us to verify that any observation made is not threshold dependent.

Figure 9 highlights a critical difference between HTE and HTR: the largest first principal normal stress difference occurs in plumes in the HTE case (*ζ* ≳ 0.01), whereas $NpNp,th$ is significantly larger in boundary layers than plumes in the HTR case for all thresholds. The conditional averaging of *N*_{p} confirms the relative importance of plumes and boundary layers surmised from the balance of energy depending on the regime of heat transfer modification and the structure of the flow.

As mentioned above, the first principal normal stress difference generates a force orthogonal to the principal direction of polymer stress, which is the main driving force underlying the Weissenberg effect. A similar mechanism might be at play here: the thin regions of high principal normal stress difference observed in Fig. 7 cause the streamlines of the convection cells to be pushed inward, i.e., toward the center of the cell, causing the flow to slow down at the center of the plume as shown in the 2D simulations (Fig. 1). On the other hand, HTE and HTR might also be caused by the large extensional viscosity. Finally, a combination of both mechanisms might best explain HTM.

In HTR, the magnitude of *N*_{p} is such that the plumes fan out considerably more than those in HTE (Fig. 6). This mechanism is similar to that of drag reduction where polymers apply a negative torque^{33} through their significant stretching in the extensional flows created by, and located on the edge of, quasi-streamwise vortices.^{20,21}

Results clearly demonstrate that, for moderate to large Weissenberg numbers, short polymers generate HTE and longer polymers result in HTR. The parameter *L* emerges as a, if not the, critical factor in HTM. This is not surprising since the extensional viscosity reaches its maximum for Deborah numbers $\lambda \epsilon \u0307\u22731$, where $\epsilon \u0307$ is the extensional rate. In our flows, it was verified that, in regions of large *N*_{p}, the local Deborah number is larger than unity (not shown). The maximum extensional viscosity for the FENE-P model asymptotes to a value proportional to *L*^{2} (see the work of Tamano *et al*.^{39}). Consequently, low to moderate extensional viscosity stabilizes the flow and reduces the size of convection cells but does not slow down the circulation of convection cells enough to reverse HTE. Longer polymers can create larger extensional viscosity and thus act as a break on the circulation of convection cells, thereby leading to HTR. Note that whether extensional viscosity, first normal stress difference (Weissenberg effect), or a combination of both is responsible for the observable repulsion between two adjacent convection cells remains an open question.

## VI. DISCUSSION OF THE REALISTIC VALUE OF OUR SIMULATION AND THE CONCERN OF AHLERS AND NIKOLAENKO

Despite clear experimental evidence of heat transfer reduction, Ahlers and Nikolaenko^{2} expressed concerns regarding the inconsistency between flow and polymer time scales. These authors argued that the turbulent time scale of RBC flow is too slow compared with the relaxation time scale of the polymer solution. In wall-bounded flows, Lumley^{6} defined that the necessary condition for polymers to interact with turbulence is that the relaxation time scale should be equal or larger than the time scale of near-wall turbulence, i.e., the inverse of the wall shear. Interestingly, this time scale is effectively the Kolmogorov time scale close to the wall.^{40}

Instead of using the large scale of RBC, as Ahlers and Nikolaenko did, we propose that the relevant time scale should be that of the dissipative structures of RBC turbulence in plumes and boundary layers. Our proposition certainly makes sense within the framework of our mechanism of heat transfer enhancement or reduction, since this mechanism is based on extensional flows arising from strong and local velocity gradients, not on the large scale motion of convection cells. Calling *σ*^{*} the magnitude of velocity gradients which significantly stretch polymers, we can estimate the correlation for the dissipation rate of turbulent kinetic energy in boundary layers, defined by^{34}

Using the appropriate scaling for natural convection in water for the Reynolds number *Re* = 0.039, *Ra*^{1/2}*Pr*^{−5/6}^{34} yields

For Ahlers and Nikolaenko’s experiment, using *Ra* = 7 × 10^{10}, *Pr* = 4.34, *ν* = 10^{−6} m^{2}/s, and *h* = 0.5 m, the estimated value of *σ*^{*} is 0.09 s^{−1}. For the highest relaxation time reported recently for drag reduction,^{7} the Weissenberg number based on *σ*^{*} is of the order of unity, which is sufficient to achieve a few orders in magnitude increase in extensional viscosity. Xie *et al*.^{4} estimated a Weissenberg of the order of 0.1 based on the averaged Kolmogorov time scale measured from particle image velocimetry but recognized that fluctuations of velocity gradients are likely to create shorter time scales, therefore larger *Wi*.

It is important to note that this analysis suffers from (a) the absence of scaling laws for extensional velocity gradients at the origin of boundary layers and (b) (deriving from a) the fact that the proposed scaling is based on an average of $(\u2202ui/\u2202xj)2$ in the boundary layer of convection cells. The statistical average is bound to underestimate the relevant value of *σ*^{*} in the region of maximum extension. When a similar analysis is applied to our simulation with *h* = 1 cm and a temperature gradient of 2.64 K (for *Ra* = 10^{5}), it is found that *σ*^{*} ≈ 0.2 s^{−1}. HTM should be observable for relaxation times of the order of unity or more, which is less than the Xanthan Gum solution used by Escudier *et al*.^{7} and comparable to typical PEO or polyacrylamide (PAM)^{41} solutions used in many experiments of polymer drag reduction.^{8}

## VII. CONCLUSION

At the low Rayleigh number of our simulation, polymers yield heat transfer enhancement (HTE) and reduction (HTR) in Rayleigh–Bénard convection flows through the same mechanism: large first principal normal stress difference concentrated in sheet-like regions at the center of plumes. This phenomenon causes the horizontal length of convection cells to be much shorter than that in the Newtonian flow. In the case of HTE, the plume velocity, measured by its rms, is similar to the Newtonian flow. It is the sheer increase in the number of cells, therefore in the volume occupied by plumes, that drives HTE. The key difference between HTE and HTR is in the boundary layers. At HTE, the effect of polymer stress as measured by its energetic contribution is threefold smaller than that of the polymer energetic contribution in plumes. At HTR, the boundary layer contribution is nearly 50% larger than the plume contribution. The net effect is a dramatic reduction in the circulation speed of the convection cells, resulting in heat transfer reduction as large as 30% in the parameter space investigated in this paper.

An interesting observation is the Weissenberg-like effect of polymers on the convection cells. Regions of first-normal stress *N*_{p} are concentrated in plumes for HTE and plumes and boundary layers for HTR. In plumes, *N*_{p} causes the plumes to tilt inward, i.e., toward the center of convection cells. Much like for polymer drag reduction where polymers apply a break via negative torque, polymers slow down the circulation of convection cells at HTR.

The present study calls for further investigation of the dynamics of natural convection in the presence of polymers. Whereas a direct numerical simulation of Ahlers and Nikolaenko^{2} or Xia’s group^{3,4} would require very large computing resources, exploration of the lower-Rayleigh number space should be encouraged. Such investigations may be confined to numerical simulations, as the relevant flow time scales may be too long for realistic polymer solution to achieve Weissenberg numbers of the order of unity. The topological change in the convection cells suggests that polymers may drive the flow toward topologies observed at lower Rayleigh numbers, such as the hexagonal cells in Fig. 4. At HTE, there seems to be a maximum of heat transfer at low Weissenberg numbers, which is a function of *L*. With the present range of *Wi*, it is not possible to discern the possible existence of asymptotic states at large *Wi*. The three main questions are as follows: (i) Does the HTM (Fig. 2) at different short *L* converge onto one curve? (ii) If so, does this curve asymptote to HTM = 0? (iii) Does HTM for large *L*, i.e., for HTR, asymptote to a plateau at very large *Wi*? The effect of shear thinning should also be investigated considering *β* → 1. The challenge with simulations at large *L* and/or *Wi* is that the numerical solution becomes stiffer, thus requiring higher resolution and lower time steps. Finally, the 2D simulations show clear evidence of chaos or elasto-inertial turbulence for large *L*. This matter is under investigation and will be reported in subsequent publications.

From the improved knowledge of the interactions between polymers and convection cells, the present study suggests a possible explanation for the heat transfer reduction and enhancement observed experimentally by Ahlers and Nikolaenko.^{2–4} Ahlers and Nikolaenko expressed concerns since the time scale of convection is much larger than the relaxation time scale of polymer solutions. Our proposition, aligned with the argument of Xie *et al*.,^{4} is to consider the time scale of velocity gradients obtained from scaling argument, which is of the same order of magnitude as the relaxation time scale of typical polymer solutions. Considering the challenge of viscoelastic simulations at high Reynolds (or Rayleigh) numbers, more experiments are very much needed to understand the range of Rayleigh numbers, Weissenberg numbers, and concentrations over which polymers modify heat transfer, as well as an investigation of possible asymptotic states.

## ACKNOWLEDGMENTS

Y. D. acknowledges the financial support of the National Science Foundation award CBET-Fluid Dynamics (Grant No. 1805636).

## REFERENCES

Polyethylene oxide and polyacrylamide