<figure><img src="image_1.png" /></figure>

HAL
open science

# Oscillating onset of the Rayleigh-Bénard convection with viscoelastic fluids in a slightly tilted cavity

Xin Zheng, Shihe Xin, M'Hamed Boutaous, Chao Wang, Dennis A Siginer,

Wei-Hua Cai

► To cite this version:

Xin Zheng, Shihe Xin, M'Hamed Boutaous, Chao Wang, Dennis A Siginer, et al. Oscillating onset of the Rayleigh-Bénard convection with viscoelastic fluids in a slightly tilted cavity. Physics of Fluids, 2023, 35 (2), pp.023107. 10.1063/5.0137501. hal-03992564

HAL Id: hal-03992564

https://cnrs.hal.science/hal-03992564

Submitted on 16 Feb 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Oscillating onset of the Rayleigh-Bénard convection with viscoelastic fluids in a slightly tilted cavity

Xin Zheng (郑鑫),<fnref n="1" />,<fnref n="2" /> Shihe XIN,<fnref n="2" /> M'hamed BOUTAOUS,<fnref n="2" /> Chao WANG* (王超),<fnref n="1" /> Dennis A. SIGINER,<fnref n="3" /> and Wei-Hua Cai* (蔡伟华)<fnref n="4" />

The oscillating onset of the Rayleigh-Bénard convection (RBC) with viscoelastic fluids in a slightly tilted 2-dimension (2D) rectangular cavity with aspect ratio $Γ = 2$ was investigated for the first time via direct numerical simulation. A series of simulations were conducted in the plane of the Rayleigh number ($Ra$) and the tilt angle ($α \in [0°, 5°]$) with three Weissenberg numbers ($Wi = (0.1, 0.15, 0.2)$) at a fixed Prandtl number $Pr = 7.0$. The evolutionary path of the oscillating convection onset in the ($Wi, α$)-plane was determined and corresponding complex flow structures were observed. The inclination of the box delays the onset of the oscillations and the corresponding Rayleigh number $Rac$ as compared to the horizontal configuration. Oscillating flow structures acquire the attributes of a traveling wave. A specific feature of the oscillating convection in the case of the horizontal cavity, the periodicity in space and time exists in the inclined box case as well. But, the evolution of the oscillatory flow structure is very different from the horizontal case in that the counter-clockwise cell assimilates the clockwise cell [Physical Review Fluids 7, 023301 (2022)].

## I. INTRODUCTION

Rayleigh-Bénard convection (RBC) is a simplified non-linear dynamic model of a slew of physical processes including thermal storage and management of electrical equipment, atmospheric physics, chemical processes and mantle flow among many others. An exhaustive review of the reported research in the literature can be found in references$^{1,2}$. Most studies about RBC were carried out with Newtonian fluids, and comparatively limited work was done with non-Newtonian fluids. However, this situation has been changed in recent decades.

The first study on the Rayleigh-Bénard convection with viscoelastic fluids (VRBC) can be traced back to Green III$^{3}$, who experimentally investigated heat transfer characteristics in VRBC and found the convective oscillating phenomena. His pioneering work was followed by a series of studies to gain a deeper understanding of several aspects of the VRBC$^{4–7}$ together with investigations on the flow reversal in turbulent flow$^{8}$, overstability$^{9–11}$, chaos$^{10}$, properties of the RBC through linear and nonlinear stability analyses$^{12,13}$, the effect of surfactants on the RBC$^{14}$, experimental findings on the turbulent$^{15}$ and yield stress fluids VRBC$^{16,17}$ and on heat transfer properties in turbulent VRBC$^{18}$. Direct numerical simulation (DNS) added much later to the tool kit for the study of both laminar and turbulent VRBC made possible remarkable contributions towards a thorough mapping of this flow$^{19–21}$.

It should be noted that all the progress made in the complete mapping out of the VRBC is related to the horizontal flow systems. However, horizontally placed flow systems represent

ideal situations and seldom occur in nature and industry. They are difficult to realize even in laboratory settings. The available literature on the RBC in tilted flow systems with Newtonian fluids shows that any angle of inclination will bring large departures from all aspects of the RBC determined for the horizontally positioned systems including pattern selection, heat transfer enhancement and turbulent flow reversal$^{22–25}$. Torres et al.$^{25}$ determine the evolution of the flow pattern with increasing angle of inclination from horizontal to vertical in a 3-dimension (3D) cube heated and cooled from the side walls with $Pr = 0.71$. In their study, the region of flow bifurcation under changing inclination is near the convection onset. Wang et al.$^{26}$ studied the coexistence of the flow states in 2-dimension (2D) RBC in tilted boxes with Newtonian fluids (NRBC) and $Pr = 0.71$ at large aspect ratios and angles of inclination $α \in [0, π]$. Shishkina and Horn$^{23}$, Zwirner and Shishkina$^{24}$, Wei and Xia$^{27}$, Vasil’ev et al.$^{28}$ reported on the heat transfer characteristics in tilted NRBC with rich parameters variations for Ra, Pr, $Γ$ and $α$. Wang et al.$^{29}$ investigated using DNS the effect of $α$ on flow reversal in turbulent RBC in 2D cavity with two aspect ratios $Γ = (1, 2)$, and mapped out the flow states in $α – Pr$ space for $Ra = 10^7$ and $α – Ra$ space for $Pr = 0.3$. They pointed out that the inclination has the opposite effects on triggering flow reversal for the cases with $Γ = 1, 2$. It promotes and suppresses it for $Γ = 1$ and for $Γ = 2$, respectively.

Compared to the Newtonian counterpart, the reports about the tilted Rayleigh-Bénard convection with non-Newtonian fluids are very few. To our knowledge the only available studies are by Vinogradov, Khezzar, and Signiger$^{30}$ and Khezzar, Signiger, and Vinogradov$^{31}$. They numerically simulated 2D

<footnote n="1">¹) College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, China</footnote>

<footnote n="2">²) Univ Lyon, CNRS, INSA-Lyon, Université Claude Bernard Lyon 1, CETHIL UMR5008, F-69621, Villeurbanne Cedex, France</footnote>

<footnote n="3">³) Departamento de Ingeniería Mecánica, Universidad de Santiago de Chile, Avenida Libertador Bernardo O'Higgins No. 3363, Estación Central, Santiago, Chile</footnote>

<footnote n="4">⁴) Laboratory of Thermo-Fluid Science and Nuclear Engineering, Northeast Electric Power University, Jilin 132012, China</footnote>

<footnote>(*Electronic mail: Co-corresponding author: caiwh@neepu.edu.cn)</footnote>

<footnote>(*Electronic mail: Corresponding author: wangchao0104@hrbeu.edu.cn)</footnote>
tilted RBC with power-law fluids, and systematically examined the influences of $\alpha$, $Ra$, $Pr$, and $\Gamma$ on the heat transfer capacity of RBC. Beyond that, there is no literature on the tilted viscoelastic RBC (TVRBC) with viscoelastic fluids with a nonlinear constitution. It should also be pointed out that from a strictly mathematical point of view the existence and uniqueness of the tilted flow configuration is an open question which has not been explored in depth except for the recent publications of Neustupa and Sigier$^{32,33}$.

There is clearly a distinction between weakly elastic VRBC which shows similarities to the NRBC and strongly elastic VRBC with much more interesting behavior, regular flow reversal for instance. The distinction extends to the convection onset and the corresponding pattern selection which depends on the rheological parameters $Wi$ and $\beta^{34,35}$. Similar dependency is also found by Park and Ryu$^{36}$. Park, Shin, and Sohn$^{37}$ also arrived at similar findings through linear and nonlinear stability analysis. We investigate in this paper the effect of the inclination angle on the VRBC with an emphasis on the strongly elastic fluids. The influence of the parameters ($\alpha$, $Wi$) on TVRBC in the vicinity of the convection onset is explored in detail. We work with an aspect ratio of $\Gamma = 2$ because flow patterns in a 2 : 1 cavity with strongly elastic fluids show a rich variety compared with aspect ratios $\Gamma = 1$ and $\Gamma > 3^{35,36}$ about the critical Rayleigh number (Rac).

The rest of the paper is organized as follows. In Section II, we briefly introduce the governing equations and numerical schemes used in the present work. In Section III, the numerical results are presented. Section III D studies the flow pattern selection in the plane of ($\alpha$, $Wi$); Section III B details two types of spatio-temporal periodic oscillating convection in TVRBC; Section III C analyses the budget of viscoelastic kinetic energy for the two types of oscillating flows observed in TVRBC in order to understand the corresponding physical mechanisms. Finally, a short summary and future directions are presented in Section IV.

## II. MATHEMATICAL FORMULATION

The computational domain inclined by an angle $\alpha$ over the horizontal is shown in Fig. 1. The width and height of the cavity are respectively $L$ and $H$, and the aspect ratio of the domain, $\Gamma = L/H$, is fixed to 2. A temperature difference $\Delta T = (T_2 - T_1)$ is imposed between the bottom and top boundaries with $T_2 > T_1$. No heat flux or the adiabatic condition is applied to the vertical cavity walls ($\partial T / \partial x = 0$ at $x = 0$ and $x = L$). No-slip boundary condition is imposed for the velocity on all the boundaries.

### A. Governing equations

The gravitational acceleration can be decomposed into two parts, one component parallel to the bottom heated boundary and the other perpendicular to the bottom,

$$g \tilde{e}_g = g(\sin(\alpha), \cos(\alpha))e_g \quad (1)$$

<figure><img src="image_3.png" /><figcaption>FIG. 1: Schematic diagram of the Rayleigh-Bénard convection with the viscoelastic fluids in an inclined cavity.</figcaption></figure>

where $e_g$ is the unit vector in the gravitational direction and $\tilde{e}_g = (\sin(\alpha), \cos(\alpha))e_g$ consists of two components of $e_g$ in the parallel and perpendicular directions to the y axis. Using Eq. (1) in the momentum balance equation, the field equations with the Oberbeck-Boussinesq approximation read as,

$$\nabla \cdot \mathbf{u} = 0 \quad (2)$$

$$\rho_0 \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla) \mathbf{u} \right) = - \nabla \bar{p} + \nabla \cdot \boldsymbol{\sigma} + \rho_0 \alpha (T - T_0) g \tilde{e}_g \quad (3)$$

$$\frac{\partial T}{\partial t} + (\mathbf{u} \cdot \nabla)T = \frac{k}{\rho_0 C_p} \Delta T + 2\mu_s D : \nabla \mathbf{u} + \boldsymbol{\tau}_p : \nabla \mathbf{u} \quad (4)$$

where $\mathbf{u} = (u_1, u_2)$ stands for the velocity vector with $u_1$ and $u_2$ representing the velocity components in the x and y directions, respectively. $\bar{p} = p - \rho_0 g \tilde{e}_g$, $T$, $k$, $C_p$ and $\rho_0$ stand for the pressure, the temperature, the thermal conductivity, the specific heat and the liquid density at the temperature $T_0 = (T_2 + T_1)/2$, respectively. The last two terms in Eq. (4) represent viscous dissipation and elastic dissipation, further discussed in detail in Section II D. Total stress $\boldsymbol{\sigma}$ in the momentum balance Eq. 8 is decomposed into its components $\boldsymbol{\sigma} = \boldsymbol{\tau}_s + \boldsymbol{\tau}_p$ the Newtonian solvent contribution $\boldsymbol{\tau}_s = 2\mu_s D$ and the viscoelastic polymer contribution $\boldsymbol{\tau}_p$ with $\mu_s$ the solvent viscosity. $D = \frac{1}{2}(\nabla u + \nabla u^T)$ is the deformation rate tensor and $\boldsymbol{\tau}_p$ is the extra-stress tensor. In present work, the Phan-Thierry-Tanner constitutive model (PTT)$^{38,39}$ is chosen to describe $\boldsymbol{\tau}_p$:

$$\nabla \cdot \boldsymbol{\tau}_p = -\frac{1}{\lambda} \boldsymbol{\tau}_p + 2 \frac{\mu_p}{\lambda} D - \frac{\varepsilon}{\mu_p} r(\boldsymbol{\tau}_p) \boldsymbol{\tau}_p - \xi (D \boldsymbol{\tau}_p + \boldsymbol{\tau}_p D) \quad (5)$$

where $\lambda$ is the relaxation time, $\mu_p$ is the polymer viscosity at zero shear rate, rheological parameters $\varepsilon$ and $\xi$ represent the elongational and slippage behaviors of the macromolecules. $\nabla \cdot \boldsymbol{\tau}_p$ indicates the contravariant convected time derivative of the elastic stress,

$$\nabla \cdot \boldsymbol{\tau}_p = \frac{\partial \boldsymbol{\tau}_p}{\partial t} + (\mathbf{u} \cdot \nabla) \boldsymbol{\tau}_p - \nabla u^T \cdot \boldsymbol{\tau}_p - \boldsymbol{\tau}_p \cdot \nabla u \quad (6)$$
The following scale factors are used to non-dimensionalize the governing equations:

$$x^* = \frac{x}{H}, \quad t^* = \frac{U_c}{H}, \quad \mathbf{u}^* = \frac{\mathbf{u}}{U_c}, \quad T^* = \frac{T - T_0}{T_2 - T_1}, \\ p^* = \frac{p}{\rho_0 U_c^2} \quad \text{and} \quad \tau_p^* = \frac{\tau_p}{\rho_0 U_c^2}$$

where the characteristic velocity is defined as $U_c = \frac{\kappa}{H} \sqrt{Ra}$. In order to simplify the notation, we drop hereafter (*) from all the dimensionless variables. Hence, the dimensionless equations governing 2D tilted Rayleigh-Bénard convection with viscoelastic fluids and Oberbeck-Boussinesq approximation are rewritten as:

$$\nabla \cdot \mathbf{u} = 0 \qquad (7)$$

$$\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\nabla p + \beta \frac{Pr}{\sqrt{Ra}} \Delta \mathbf{u} + \nabla \cdot \tau_p + (\sin(\alpha), \cos(\alpha)) Pr Te_g \qquad (8)$$

$$\frac{\nabla}{\tau_p} + \frac{\tau_p}{Wi\sqrt{Ra}} - 2\frac{1-\beta}{Ma^2}D = -\epsilon\frac{\sqrt{Ra}}{(1-\beta)Pr}tr(\tau_p)\tau_p - \xi(D\tau_p + \tau_pD) \qquad (9)$$

$$\frac{\partial T}{\partial t} + (\mathbf{u} \cdot \nabla)T = \frac{1}{\sqrt{Ra}}\Delta T + 2\beta \frac{EcPr}{\sqrt{Ra}} D: \nabla \mathbf{u} + Ec\tau_p: \nabla \mathbf{u} \qquad (10)$$

where $\beta = \mu_s/\mu_0$ is the ratio of the solvent viscosity to the total viscosity $\mu_0 = \mu_s + \mu_p$. The boundary conditions read as,

* At the top boundary $y = 1$: $u_1 = u_2 = 0$; $T = 0$

* At the bottom boundary $y = 0$: $u_1 = u_2 = 0$; $T = 1$

* At the side boundaries $x = 0, 2$: $u_1 = u_2 = 0$; $\frac{\partial T}{\partial x} = 0$

The dimensionless parameters include: the Rayleigh number $Ra = \frac{\alpha g \Delta T H^3}{\kappa}$, the Prandtl number $Pr = \mu_0 C_p / k$, the Weissenberg number $Wi = \frac{\lambda \kappa}{H^2}$ and the Mach number $Ma = \frac{\sqrt{Ra Wi}}{Pr}$ describes the ratio of the inertia force and the elastic force. The intensity of the elasticity is described by the elasticity number $E = Re / Wi = \sqrt{Ra / Pr Wi}$. In Eq. (10), the dimensionless Eckert number $Ec = \frac{U_c^2}{C_p \Delta T}$ quantifies the relationship between the kinetic energy and the enthalpy, and is used to characterize the extent of the influence of the heat dissipation. For thermal convection, the Nusselt number is used to show the intensity of convective heat transfer, $Nu = - \frac{\partial T}{\partial y}|_{y=0}$. The corresponding spatially averaged Nusselt number $Nu_s$ is defined as:

$$Nu_s = - \frac{1}{L} \int_{0}^{L} \left. \frac{\partial T}{\partial y} \right|_{y=0} dx$$

Tab. I summarises the dimensionless parameters characterizing the TVRBC.

<table><caption>TABLE I: Summary of numbers characterizing the TVRBC</caption><thead><tr><th>Ra = αgΔTH<sup>3</sup>/√κ</th><th>Pr = μ<sub>0</sub>C<sub>p</sub>/k</th><th>Re = √Ra/Pr</th></tr></thead><tbody><tr><td>β = μ<sub>s</sub>/μ<sub>0</sub></td><td>Wi = λκ/H<sup>2</sup></td><td>E = √Ra/PrWi</td></tr><tr><td>Nu = -∂T/∂y|<sub>y=0</sub></td><td>Ec = U<sub>c</sub><sup>2</sup>/C<sub>p</sub>ΔT</td><td></td></tr></tbody></table>

## B. Working fluids

The Phan-Thien and Tanner (PTT) model is widely used to describe the response of viscoelastic fluids under stress. There are two special constitutive parameters in the PTT model, $\epsilon$ that governs the elongational response of the fluid and $\xi$ that controls the non-affine movement of the molecular long chain lattice relative to the flow of the solvent and quantifies the shear thinning capability of the fluid in the absence of the former parameter $\epsilon$. If both parameters are present in the constitutive structure as is the case of the PTT fluid shear thinning is controlled by both proportionally to the magnitude of the parameters. The numerical values of the constitutive parameters $\xi$ and $\epsilon$ in the PTT model obtained by experimental measurements are usually in the range of $\epsilon \in [0, 0.5]$ and $\xi \in [0, 0.7]$ for different concentrated industrial materials such as polyethylene (LDPE), molten polyethylene (HDPE) and polyisobutylene (PIB)$^{40-44}$. Moderate values of the parameters are adopted in the present work $\epsilon = 0.1$ and $\xi = 0.05$, which is close to the solutions of 2.5% polyisobutylene$^{40}$. Other dimensionless parameters are set at: $Pr = 7.0$, $\beta = 0.2$, $Wi = (0.1, 0.15, 0.2)$ and $\alpha \in [0^\circ, 5^\circ]$.

The definition of the Weissenberg number $Wi$ used in this paper is different than the definition used by other researchers such as Cheng et al.$^{20}$ because the reference velocity $U_c = \frac{\kappa}{H} \sqrt{Ra}$ adopted in this paper depends on $Ra$. This definition implies that any change in $Ra$ leads to a change of $Wi$. The conventional definition of $Wi$ is independent of $Ra$, $Wi = \kappa \lambda / H^2$. For example, for $Ra = 1600$ and $Wi = 0.1$ in the present work, the equivalent $Wi$ in the work by Cheng et al.$^{20}$ (in their work $U_c = \sqrt{\alpha g H \Delta T}$) is $Wi = 10.58$.

## C. Numerical procedures

A recently developed in-house solver was adopted to solve the governing equations. The capability of this solver to simulate VRBC has been shown in our previous work$^{21,34,35,45}$. Various validation tests have been done. The highlight of this solver is the quasi-linear treatment of the hyperbolic terms in the momentum equation and the viscoelastic constitutive equation. As we know, extreme numerical instability can be brought by the hyperbolicity of the governing equation of the viscoelastic constitutive equation, which easily leads to numerical divergence. The quasi-linear treatment is used to reorganize the hyperbolic terms in the momentum equation and the viscoelastic constitutive equation, which makes them can be discretized by the HOUC to concentrate numerical dissipation toward high wavenumbers. A second-order semi-implicit time scheme is used. The implicit terms are mass conservation, pressure gradient, molecular diffusion, relaxation term,
and thermal diffusion; other terms including the quasi-linear forms are explicit. The time integration was carried out using a second-order backward differential formula (BDF2). In space, second-order central differencing is applied to most of the terms except for the quasi-linear parts which are treated in the eigenspace of $A_i$ by a third-order High-Order Upstream Central scheme (HOUC3)46,47 according to the sign of each eigenvalue. Finally, the velocity-pressure coupling is treated by the projection method.

### D. Viscous and elastic dissipation

The magnitudes of the viscous and elastic dissipations in the VRBC are evaluated through the last two terms on the RHS of Eq. (10). Fig. 2 shows the spatially-averaged dimensionless viscous and elastic dissipation for the case with $\alpha = 0^{\circ}$, $Pr = 7.0$, $Ra = 1600$, $\beta = 0.2$ and $Wi = 0.1$. Viscous dissipation remains always positive and is thus a positive source term for the temperature field whereas elastic dissipation is positive for most of the time and negative for a short time during one period. Therefore elastic dissipation could lead to decreasing temperatures. $Ec$ can be expressed as

$$Ec = \frac{U_c^2}{C_p \Delta T} = \frac{(\frac{\tilde{\alpha}}{H} \sqrt{Ra})^2}{C_p \Delta T} = \frac{\tilde{\alpha}^2 Ra}{H^2 C_p \Delta T} \quad (11)$$

Substituting the numerical values of the parameters for the working fluid into Eq. (11) the magnitudes of $Ec$ and $2\beta \frac{EcPr}{\sqrt{Ra}}$ in Eq. (10) determined to be very small $O(10^{-10})$ and $O(10^{-11})$, respectively. $2\beta \frac{EcPr}{\sqrt{Ra}} D: \nabla u$ and $Ec \tau_p: \nabla u$ are also very small. All that means that the effects of the viscous and elastic dissipations can be neglected.

### III. RESULTS AND DISCUSSION

#### A. Pattern selection and bifurcation diagrams

In Newtonian RBC it is well known that any inclination of the cavity cancels the pure conduction regime and results in a weak 1-cell convective flow for small inclination angles. With increasing Rayleigh number, the 1-cell convective pattern is transformed continuously into 2-cell convective pattern through an imperfect pitchfork bifurcation. The same situation is observed in TVRBC for big $\beta^{48}$ because polymeric contribution to the viscosity is small and Newtonian solvent plays an overwhelmingly important role leading to behavior similar to those of the Newtonian and weakly viscoelastic fluids$^{25}$.

In the present work we are interested in the cases with small $\beta$. Steady-state convection with a single cell ($T1S$) rotating in the inclination direction always exists when $Ra$ is above zero for small values of $\beta = 0.2$ and $\alpha = 2^{\circ}$ and the remaining parameters set at $Pr = 7.0$, $Wi = 0.2$, $\epsilon = 0.1$ and $\xi = 0.05$. Fig. 3 displays streamlines of the 1-cell convection obtained for $Ra = 800, 960$ and $1120$: the weak 1-cell convection flows in the counter-clockwise direction and the temperature field

is mainly conductive. With further increase of the Rayleigh number, oscillating convective flows are observed as will be detailed below.

To clearly understand the flow pattern transition in TVRBC, we plotted the bifurcation diagram along $Ra$ at fixed $\alpha = 2^{\circ}$ in detail. In the tested cases, rheological parameters are $Pr = 7.0$, $\beta = 0.2$, $Wi = 0.1$, $\epsilon = 0.1$ and $\xi = 0.05$. Fig. 4 shows the bifurcation diagrams in the range of ($600 < Ra < 1800$) for TVRBC ($\alpha = 2^{\circ}$) and HVRBC ($\alpha = 0^{\circ}$). The bifurcation diagram in HVRBC follows grey lines for easy comparison. The curves with other colors represent different flow patterns in tilted cases, (—) corresponds to stable convection with counterclockwise single-cell ($T1S$), (—) corresponds to time-dependent oscillating convection with three cells ($T3R$), and (—) to time-dependent oscillating convection with two cells with horizontal movement ($T2R$). The flow pattern branches and bifurcations are also annotated in the figure in the corresponding colors.

There are two solution branches (positive and negative) for each pattern of $H2R$ and $H3R$. Taking $H2R$ as an example, the positive and negative values of the branches are the maximum positive and negative velocities respectively in one oscillation period. The time-dependent $H2R$ and $H3R$ have both branches starting from the pure conduction state distinguished by the Hopf bifurcations $P1$ and $P2$. The colored

<figure><img src="image_4.png" /><figcaption>(a) Time evolution of $D: \nabla u$ integrated in space</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) Time evolution of $\tau_p: \nabla u$ integrated in space</figcaption></figure>

<caption>FIG. 2: (a) Viscous dissipation without the coefficient ($2\beta \frac{EcPr}{\sqrt{Ra}}$) and (b) elastic dissipation without the coefficient ($Ec$) when $\alpha = 0^{\circ}$, $Pr = 7.0$, $Ra = 1600$, $\beta = 0.1$, $Wi = 0.1$.</caption>
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0137501

<figure><img src="image_6.png" /><figcaption>FIG. 3: Streamlines for steady-state convection with 1 cell with $\alpha = 2^{\circ}$, $Pr = 7.0$, $\beta = 0.2$, $Wi = 0.1$, $\epsilon = 0.1$ and $\xi = 0.05$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 4: The bifurcation diagrams in HVRBC (grey) and TVRBC (colors, $\alpha = 2^{\circ}$), where the y-axis represents $u_2$ at the monitoring point $(x,y) = (7/8, 1/2)$. The simulated cases correspond to $Pr = 7.0$, $\beta = 0.2$, $Wi = 0.1$, $\epsilon = 0.1$ and $\xi = 0.05$.</figcaption></figure>

curves (TVRBC) in Fig. 4, demonstrate that a slight inclination breaks the pure conduction state found in HVRBC before $P1$ and $P2$, and forms $T1S$, which rotates in the same direction as the inclination. After that, the branches of $T3R$ grow from the bifurcation $H1$. Considering the oscillating flow configuration in HVRBC³⁵, we guess that there may be a certain $Ra$ range where $H2R$ and $H3R$ may coexist. To obtain another pattern, we take the stable solution $H2R$ of HVRBC as the

initial value and raise the inclination to $\alpha = 2^{\circ}$. The branches of $T2R$ were found, growing from the bifurcation $H2$. Due to $F_x$ the component of the buoyancy force parallel to the bottom of the enclosure, the cells in $T2R$ have the characteristic of left-right movement (detailed flow structure is in Fig. 5) compared to $H2R$. Compared to the bifurcation points $P1$ and $P2$ in HVRBC, $H1$ and $H2$ are slightly delayed. It should be noted that, due to the inclination, the branches of $T2R$ and
<figure><img src="image_8.png" /><figcaption>FIG. 5: Evolution of the periodic flow at $\alpha = 2^\circ$, $Ra = 1120$, when $Pr = 7.0$, $\beta = 0.2$, $Wi = 0.1$, $\epsilon = 0.1$ and $\xi = 0.05$. The time period $\lambda$ is equal to 13.76. The green arrows show the growth direction of the main counterclockwise cell.</figcaption></figure>

$T3R$ are not absolutely symmetric about the $x$ axis, and the maximum value of the downward velocity is slightly larger than that of the upward velocity.

## B. Oscillating convection structures

In the work of Zheng et al.³⁵, Park and Ryu³⁶, Park and Lee⁴⁹ on HVRBC, an oscillating convection, which exhibits strong spatio-temporal periodicity, takes place with $\Gamma = 2$ and a certain range of $\beta$ and $Wi$. The time-dependent convective cells show standing wave characteristics in space in the pattern of $H2R$ (two cells) meaning that there is no mass exchange between cells. However, the standing wave characteristics vanish in the pattern of $H3R$. The main reason for the periodic oscillating process is the alternating growth (phase difference in time and space) of the velocity field and the elastic stress field. The generation and growth of the new vortex is driven by the elastic stress and buoyancy flux³⁴. The cells in the oscillating convection are rotating in the clockwise and counterclockwise directions, and each circulation is of equal importance during the whole process in HVRBC.

Interestingly, the introduction of an inclination breaks the symmetry between the cells and a new 2-cell periodic oscillating convection ($T2R$) is formed. Take the case with $\alpha = 2^\circ$,

$Wi = 0.1$ and $Ra = 1120$ as an example. The evolution of the velocity vector of $T2R$ in one period in TVRBC is described in Fig. 5. The time stages ($t1 - t5$) and ($t5 - t9$) present respectively the two half-periods, in which the velocity field is completely reversed twice. From Fig. 5a to Fig. 5c: the counterclockwise and clockwise cells are of about the same size at the initial time $t1$, after that the counter-clockwise cell (the left cell) grows and moves to the right, and almost occupies the whole cavity at $t3$. From Fig. 5c to Fig. 5d: the cell keeps moving to the right and frees the left part for another secondary cell which is completely developed in Fig. 5e. Then the left cell is weakened again (Fig. 5f) and the main cell (the right cell) grows and moves from the right (Fig. 5g) to the left (Fig. 5h). Finally, the right part of the cavity is freed and another cell is formed again in Figs. 5h and 5i. The repeated occurrence of the above process forms the time-dependent oscillating convection in the tilted cavity ($T2R$). The main counter-clockwise cell is always stronger in terms of kinetic energy than the clockwise one. In the reverse process, it becomes bigger and oscillates alternatively from the left to the right side and from the right to the left side. We call the cell in counter-clockwise rotation the main cell and the one in clockwise rotation the secondary cell.

The evolution of the tilted oscillating convection $T3R$ is shown in Fig. 6. The main cells (the two counter-clockwise
<figure><img src="image_9.png" /><figcaption>FIG. 6: Periodic evolution of the periodic flow at $\alpha = 2^\circ$, and $Ra = 2080$, when $Pr = 7.0$, $\beta = 0.2$, $Wi = 0.1$, $\epsilon = 0.1$ and $\xi = 0.05$. The time period $\lambda$ is equal to 12.01. The green arrows show the growth direction of the main counterclockwise cells.</figcaption></figure>

cells) near the cavity ends grow and merge to form one large main cell approximately centered on the cavity, which then is broken down to form two main cells and initiate a reversal. In Fig. 6a, the main cells are located near the left and the right ends of the cavity. The secondary cell in the central part is gradually weakened and absorbed into the main cells followed by the cell on the right overtaking and merging with the main cell on the left (Figs. 6b, 6c, 6d). The main cell occupies the full cavity at t4 (Fig. 6d). The secondary cells form near the right and left ends of the enclosure at t5 or at half period (Fig. 6e) They start growing weaker ceding space to the growing central main cell at t6 (Fig. 6f). The main cell continues to grow and occupies the full cavity at t7 (Fig. 6g). It is broken into two cells near the left and right ends at t8 (Fig. 6h) with the secondary cell forming again in the central part of the cavity. The flow configuration returns to its original setting at t9 and the oscillation cycle continues.

In the reversal process of T3R, the two main cells located near the cavity ends grow and approach each other and merge in the central part of the cavity and two secondary cells are formed near the cavity ends on the one hand and on the other hand the large main cell grows to fill the space close to the cavity and is broken into two cells filling the cavity ends and a secondary cell is formed in the central part of the cavity.

Thus, viscoelasticity drives the left-right movement of the

oscillating convection and the periodic assimilation of the secondary cell by the primary cell.

### C. Kinetic energy transfer

The flow pattern and structure of the tilted oscillating convection have been discussed in the previous section. The viscoelastic kinetic energy budget (VKE)$^{20,50}$, similar to the well-known turbulent kinetic energy budget (TKE), is a powerful tool at our disposal to gain a deeper understanding of the oscillating convection. The actions of different forms of viscoelastic kinetic energy in the system can be clearly observed through VKE, to gain a deeper understanding of the oscillations in TVRBC. The equation for the instantaneous VKE is derived from the momentum balance Eq. (8) by integrating the
components over the computational domain $\omega$.

$$ \oint_{\omega} (\frac{\partial \epsilon}{\partial t}) d\omega = \underbrace{\oint_{\omega} \left( -\frac{\partial (pu_j)}{\partial x_j} + \beta \frac{Pr}{\sqrt{Ra}} \frac{\partial^2 e_{ij}}{\partial x_j^2} \right) d\omega}_{\Phi_D} + \underbrace{\oint_{\omega} \left( -\beta \frac{Pr}{\sqrt{Ra}} \frac{\partial u_i}{\partial x_j} \frac{\partial u_i}{\partial x_j} \right) d\omega}_{\Phi_V} + \underbrace{\oint_{\omega} \left( \tau_{ij} \frac{\partial u_i}{\partial x_j} - \frac{\partial (u_i \tau_{ij})}{\partial x_j} \right) d\omega}_{\Phi_G} + \underbrace{\oint_{\omega} \left( (\sin(\alpha), \cos(\alpha)) Pr T e_g u_i \right) d\omega}_{\Phi_F} \quad (12) $$

where, $e = \frac{1}{2} u_i u_i$ denotes the kinetic energy. In the RHS of Eq. (12), the first term ($\Phi_D$) represents the pressure diffusion and molecular viscous transport; the second term ($\Phi_V$) represents the viscous dissipation, the work done by the Newtonian solvent fluid due to the action of shear forces and transformed into heat; the third term ($\Phi_G$) denotes the energy transfer between the flow and the macromolecules due to the stretching and relaxation of the polymeric chains. $\frac{\partial(u_i \tau_{ij})}{\partial x_j}$ and $\tau_{ij} \frac{\partial u_i}{\partial x_j}$ denote the elastic dissipation and the interaction between the flow and polymeric chains, respectively. The last term ($\Phi_F$) is the input from the work done by the buoyancy using the Boussinesq approximation. The kinetic energy map of a TVRBC system is shown in Fig. 7. The kinetic energy reservoir $E_K$ exchanges with the elastic energy reservoir $E_E$ through the term $\Phi_G$ and with the buoyancy potential energy reservoir $E_F$ through the term $\Phi_F$. $\Phi_D$ and $\Phi_V$ exchange directly with $E_K$. According to the sign of each term, the kinetic energy reservoir $E_K$ can be enhanced or weakened and the elastic energy reservoir $E_E$ and the buoyancy potential energy reservoir $E_F$ can act as source terms or dissipation terms. The evolution of the energy transfers can explain very well the physics of the oscillating convection, a simplified way based on the temporal and spatial phase differences/delays among different energy forms, especially, $\Phi_F$, $\Phi_G$, and $E$. More details can be found in the work of Zheng et al. 34,35.

Fig. 8 plots the energy transfer process in the oscillating convection at fixed values of the parameters $\beta = 0.2$, $W_i = 0.1$, $Pr = 7.0$, $\epsilon = 0.1$ and $\xi = 0.05$. Figs. 8a and 8c concern the tilted cases ($\alpha = 2^\circ$), Figs. 8b and 8d concern the horizontal cases ($\alpha = 0^\circ$). The horizontal cases are plotted for comparison and to put an emphasis on the effects of the inclination on energy transfers.

Compared with the horizontal cases, the amplitude of all energy forms decreases, but the frequency increases, especially in the 2-cell cases. Also the magnitude of energies in each half-period in every convection pattern is different. Each half period in an oscillation period is characterized by momentum energy flowing from one maximum to another maximum. The energy flow/transfer in each half-period is almost of the same magnitude, for horizontal cases. The reason is that $H2R$ or $H3R$ always show spatial symmetry or temporal symmetry. Specifically, $H2R$ is symmetric about the central x axis at all times and $H3R$ is axisymmetric about the center point of the domain at all times. These symmetrical flow structure charac-

teristics lead to the symmetry of the energy transport and the time period in two half-periods.

As aforementioned in Sec. III B, there are main and secondary circulations in the domain, the main circulation is always stronger than the secondary one. The oscillating process shown in Fig. 6 clearly describes the two half-period processes including merging of the two main cells (side cells) with the secondary cell (middle cell) (Fig. 6a to 6e) and one main cell (middle cell) with two secondary cells (side cells) (Fig. 6e to 6h). In the tilted configuration, there are strong half-periods (which contain two main circulations) and weak half-periods (which contain a single secondary cell) in one period. Each energy type formed in the strong half-period is always greater than that in the weak half-period. The thermal buoyancy is the only energy input for the global convection system, Fig. 8c has two different magnitudes of $\Phi_F$ in two half-periods.

The repeated energy transfers in two half-periods also occurs in $T2R$ mode (Fig. 8a), but the physics is not the same as that in horizontal cases. There is only one main cell and one secondary cell during the entire oscillation period. The evolutionary process takes place only between these two cells (Fig. 5). From the perspective of energy transport, the energy evolution in these two half periods is unchanged.

Details of the oscillation, $u_2$ at the monitoring point $(x, y) = (7/8, 1/2)$ and the variation of $Nu_s$ versus dimensionless time for TVRBC and HVRBC (for comparison) are plotted in Fig. 9 at the fixed values of the parameters $Ra = 1760$, $\beta = 0.2$, $W_i = 0.1$, $Pr = 7.0$, $\epsilon = 0.1$ and $\xi = 0.05$. The horizontal cases were discussed in detail in the authors' previous works$^{35}$. It is observed that the time period for the velocity is twice that for $Nu_s$. The time evolution of $Nu_s$ also shows a strong-weak period in $T3R$. Globally, the inclinations tremendously weaken the heat transfer capacity in viscoelastic oscillating convection, especially in 2-cell cases. The peak of $Nu_s$ in $T2R$ and $H2R$ is basically the same as the second peak of $Nu_s$ in $T2R$ and $H2R$, respectively.

<figure><img src="image_10.png" /><figcaption>FIG. 7: Energy conversion framework for TVRBC. Three energy reservoirs exist in the convection system: the kinetic energy reservoir ($E_K$), the elastic energy reservoir ($E_E$) and the buoyancy potential energy reservoir ($E_F$). Kinetic energy transport and dissipation within the flow and the macromolecules are resulted from the interaction between buoyancy flux input ($\Phi_F$), kinetic diffusion ($\Phi_D$), viscous dissipation of kinetic energy ($\Phi_V$) and energy exchange ($\Phi_G$).</figcaption></figure>
PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0137501

<figure><img src="image_11.png" /><figcaption>FIG. 8: Comparison of the time evolution of the integrated quantities in tilted (left) and horizontal (right) VRBC. The parameters studied are Ra = 1760, β = 0.2, Wi = 0.1, Pr = 7.0, ε = 0.1 and ξ = 0.05.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 9: (a) and (b) show the variation of $u_2$ at the monitoring point $(x, y) = (7/8, 1/2)$ and the variation of $Nu_s$ versus dimensionless time in TVRBC and HVRBC, respectively. The periodic averaged values of $Nu_s$ are 1.05(T2R), 1.07(T3R), 1.13(H2R) and 1.10(H3R), respectively. In the figure, (—) and (—) relate to 2R and 3R, respectively. Ra = 1760, β = 0.2, Wi = 0.1, Pr = 7.0, ε = 0.1 and ξ = 0.05.</figcaption></figure>

D. Weissenberg number and inclination dependence

The influence of $Wi$ and $\alpha$ on the flow pattern selection in TVRBC with $Ra$ is presented in the vicinity of the imperfect bifurcation, depicting the first transition from the steady convection to the oscillating convection in the tilted configuration.

The relevant parameters $Wi = (0.1, 0.15, 0.2)$ and $\alpha \in [0^{\circ}, 5^{\circ}]$ and the corresponding flow states are summarised in Fig. 10. It should be noted that large enough $\alpha$ and $Wi$ will increase the local intensity of the elasticity (Ma) and thus may force the violation of the criterion for the convergence of this solver, namely $\tau_{it} \ge \frac{1-\beta}{Ma^2}$ depending on the local $Re$ and $Wi$.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

Δ Pure conduction
○ Reversal convection with 2 cells
+ Reversal convection with 3 cells
X Steady convection with 2 cells
☐ Steady convection with 1 cell

<caption>FIG. 10: Phase diagrams in the Ra – α plane for $Wi = (0.1, 0.15, 0.2)$. Green triangles denote pure conduction, blue circles correspond to time-dependent reversal convection with 2 cells, red crosses represent time-dependent reversal convection with 3 cells, pink crosses indicate steady convection with 2 cells and brown squares represent steady convection with a single cell.</caption>

The details are available in our previous work$^{21,35}$. Therefore, Fig. 10 does not show the cases with big $α$ and $Wi$, but, that does not affect the analysis of this paper. As a reminder, we again point out here that the values of the Weissenberg number $Wi$ as defined in this paper are equivalent to much higher values of the $Wi$ defined in the traditional way. For instance, as we already pointed out in Sec. II B the value of $Wi = 0.1$ in this paper corresponds to $Wi = 10.58$ in the work of Cheng et al.$^{20}$. Fig. 10 clarifies that there is no pure conduction state in the tilted cases. The original pure conduction state in horizontal cases breaks down immediately and convection starts as soon as a small inclination is introduced. As an example see the flow configuration changing from ($α = 0°$, $Wi = 0.1$, $Ra = 800$) to ($α = 1°$, $Wi = 0.1$, $Ra = 800$).

As $Ra$ further increases, take the case with $Wi = 0.1$ and $α = 1°$ as an example, the flow state quickly crosses the region of $T1S$, and changes to the time-dependent reversal convection with two or three rolls ($T2R$ or $T3R$). In fact, the time-dependent oscillating convection in TVRBC presents different features compared with that in HVRBC. The cells in TVRBC are characterized by left-right movement under the effect of

gravity described in detail in Sec. III B. In horizontal cases, the flow structure of the time-dependent oscillating convection with multi-cells shows similarities to the standing wave meaning no mass transfer takes place between cells and cells oscillate within their own space. Details are available in$^{35}$. The bifurcation causing the transition of $T1S$ to $T2R/T3R$ is what's called an imperfect bifurcation. Which one of $T2R$ or $T3R$ takes place first depends on the elasticity of the working fluids. For instance, $T2R$ takes place earlier when $Wi$ is larger, as shown in cases with $Wi = (1.25, 0.2)$, Figs. 10a and 10b.

With further increase of the angle of inclination, the unstable oscillating convection disappears, and only the flow pattern $T1S$ is observed, as is the case with $α = 3°$ when $Wi = 0.2$. The critical tilt angle for turbulent Newtonian Rayleigh-Bénard convection with $Ra = 4 × 10^7$, $Pr = 2$ is about $7°^{29}$. That means the horizontal component $F_x$ of the buoyancy force plays a greater role in initiating the oscillating mode with an increasing inclination angle. We observe that in Fig. 10 the flow transition from stable convection $T1S$ to time-dependent oscillating convection $T2R$ will occur first as $Wi$ increases. In addition, $T3R$ is delayed as $Wi$ increases.
In general, the increase in inclination makes the convection system more stable, but elasticity works against this effect.

## IV. CONCLUDING REMARKS

Effects of the inclination of the flow configuration on the Rayleigh-Bénard convection with strongly elastic viscoelastic fluids, highly relevant to the real-life situations are investigated numerically for the first time through an in-house solver. The computational domain is a 2D rectangular cavity with an aspect ratio of 2. The magnitude of the tilted angle is limited to a small range to observe the related complex flow bifurcations when Ra is in the vicinity of *Rac*. Large inclinations beget steady convection especially for large aspect ratios with Newtonian fluids as well<fnref n="29" />. The major findings for the tilted Rayleigh-Bénard convection with viscoelastic fluids center on the delayed imperfect bifurcation and the oscillating convection acquiring the feature of a left-right movement.

The combined impact of $(\alpha, Wi)$ on the flow pattern selection in the tilted Rayleigh-Bénard convection with viscoelastic fluids is investigated in detail. In the vicinity of *Rac* for the convection onset, three types of flow patterns, one steady mode and two unstable modes with different properties, are found for the Rayleigh-Bénard convection with strongly elastic viscoelastic fluids. The flow structures of different modes are depicted and explained in Sec. III B. The influence of the inclination on the convection structure stability is studied. The viscoelastic kinetic energy budget is very helpful in understanding the energy transfer process in the Rayleigh-Bénard convection in the tilted configuration with viscoelastic fluids. The results are important to help us understand the onset and the evolution of the oscillatory flow in the Rayleigh-Bénard convection with viscoelastic fluids in a tilted 2D box.

**Acknowledgement:** One of the authors, X.Z., thanks the Chinese Scholarship Council (CSC) for the financial support. Part of the results was obtained in the scope of PHC CAIYUANPEI project 44033SA.

<footnote n="1">¹E. Bodenschatz, W. Pesch, and G. Ahlers, "Recent developments in Rayleigh-Bénard convection," Annual review of fluid mechanics 32, 709-778 (2000).</footnote>

<footnote n="2">²P. Manneville, "Rayleigh-Bénard convection: thirty years of experimental, theoretical, and modeling work," Dynamics of spatio-temporal cellular structures, 41-65 (2006).</footnote>

<footnote n="3">³T. Green III, "Oscillating convection in an elasticviscous liquid," The Physics of Fluids 11, 1410-1412 (1968).</footnote>

<footnote n="4">⁴H. Hiromitsu, "Overstability of a viscoelastic liquid layer with internal heat generation," International journal of heat and mass transfer 29, 645-647 (1986).</footnote>

<footnote n="5">⁵R. W. Kolkka and G. R. Lerley, "On the convected linear stability of a viscoelastic Oldroyd-B fluid heated from below," Journal of non-Newtonian fluid mechanics 25, 209-237 (1987).</footnote>

<footnote n="6">⁶H. Park and H. Lee, "Hopf bifurcations of viscoelastic fluids heated from below," Journal of non-Newtonian fluid mechanics 66, 1-34 (1996).</footnote>

<footnote n="7">⁷P. Siddheshwar, G. Sekhar, and G. Jayalatha, "Effect of time-periodic vertical oscillations of the Rayleigh-Bénard system on nonlinear convection in viscoelastic liquids," Journal of non-Newtonian fluid mechanics 165, 1412-1418 (2010).</footnote>

<footnote n="8">⁸A. Castillo-Castellanos, A. Sergent, and M. Rossi, "Reversal cycle in square Rayleigh-Bénard cells in turbulent regime," Journal of Fluid Mechanics 808, 614-640 (2016).</footnote>

<footnote n="9">⁹C. M. Vest and V. S. Arpaci, "Overstability of a viscoelastic fluid layer heated from below," Journal of fluid Mechanics 36, 613-623 (1969).</footnote>

<footnote n="10">¹⁰R. E. Khayat, "Chaos and overstability in the thermal convection of viscoelastic fluids," Journal of non-Newtonian fluid mechanics 53, 227-255 (1994).</footnote>

<footnote n="11">¹¹R. E. Khayat, "Non-linear overstability in the thermal convection of viscoelastic fluids," Journal of non-Newtonian fluid mechanics 58, 331-356 (1995).</footnote>

<footnote n="12">¹²I. Eltayeb, "Nonlinear thermal convection in an elasticviscous layer heated from below," Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 356, 161-176 (1977).</footnote>

<footnote n="13">¹³H. Park and D. Ryu, "Nonlinear convective stability problems of viscoelastic fluids in finite domains," Rheologica acta 41, 427-440 (2002).</footnote>

<footnote n="14">¹⁴T. Wei, W. Cai, C. Huang, H. Zhang, W. Su, and F. Li, "The effect of surfactant solutions on flow structures in turbulent Rayleigh-Bénard convection," Thermal Science 22, 507-515 (2018).</footnote>

<footnote n="15">¹⁵W. Cai, T. Wei, X. Tang, Y. Liu, B. Li, and F. Li, "The polymer effect on turbulent Rayleigh-Bénard convection based on p/v experiments," Experimental Thermal and Fluid Science 103, 214-221 (2019).</footnote>

<footnote n="16">¹⁶C. Metivier, F. Brochard, M. Darbouli, and A. Magnin, "Oscillatory Rayleigh-Bénard convection in elasto-viscoplastic gels," Journal of Non-Newtonian Fluid Mechanics 286, 104428 (2020).</footnote>

<footnote n="17">¹⁷Z. Kebiche, C. Castelain, and T. Burghelea, "Experimental investigation of the Rayleigh-Bénard convection in a yield stress fluid," Journal of Non-Newtonian Fluid Mechanics 203, 9-23 (2014).</footnote>

<footnote n="18">¹⁸P. Wei, R. Ni, K.-Q. Xia, et al., "Enhanced and reduced heat transport in turbulent thermal convection with polymer additives," Physical Review E 86, 016325 (2012).</footnote>

<footnote n="19">¹⁹W. Cai, Q. Ye, J. Chen, H. Zhang, T. Wei, and F.-c. Li, "Study on the characteristics of Rayleigh-Bénard convection with viscoelastic fluids," in *Fluids Engineering Division Summer Meeting*, Vol. 50299 (American Society of Mechanical Engineers, 2016) p. V01BT14A006.</footnote>

<footnote n="20">²⁰J.-P. Cheng, H.-N. Zhang, W.-H. Cai, S.-N. Li, and F.-C. Li, "Effect of polymer additives on heat transport and large-scale circulation in turbulent Rayleigh-Bénard convection," Physical Review E 96, 013111 (2017).</footnote>

<footnote n="21">²¹X. Zheng, M. Boutaous, S. Xin, D. A. Siginer, F. Hagani, and R. Knikker, "A new approach to the numerical modeling of the viscoelastic Rayleigh-Bénard convection," in *ASME International Mechanical Engineering Congress and Exposition*, Vol. 59445 (American Society of Mechanical Engineers, 2019) p. V007T08A027.</footnote>

<footnote n="22">²²S.-X. Guo, S.-Q. Zhou, X.-R. Cen, L. Qu, Y.-Z. Lu, L. Sun, and X.-D. Shang, "The effect of cell tilting on turbulent thermal convection in a rectangular cell," Journal of Fluid Mechanics 762, 273-287 (2015).</footnote>

<footnote n="23">²³O. Shishkina and S. Horn, "Thermal convection in inclined cylindrical containers," Journal of Fluid Mechanics 790 (2016).</footnote>

<footnote n="24">²⁴L. Zwirner and O. Shishkina, "Confined inclined thermal convection in low-Prandtl-number fluids," Journal of Fluid Mechanics 850, 984-1008 (2018).</footnote>

<footnote n="25">²⁵J. F. Torres, D. Henry, A. Komiya, and S. Maruyama, "Transition from multiplicity to singularity of steady natural convection in a tilted cubical enclosure," Physical Review E 92, 023031 (2015).</footnote>

<footnote n="26">²⁶Q. Wang, Z.-H. Wan, R. Yan, and D.-J. Sun, "Multiple states and heat transfer in two-dimensional tilted convection with large aspect ratios," Physical Review Fluids 3, 113503 (2018).</footnote>

<footnote n="27">²⁷P. Wei and K.-Q. Xia, "Viscous boundary layer properties in turbulent thermal convection in a cylindrical cell: the effect of cell tilting," Journal of Fluid Mechanics 720, 140-168 (2013).</footnote>

<footnote n="28">²⁸A. Y. Vasil'ev, I. Kolesnichenko, A. Mamykin, P. Frick, R. Khalilov, S. Rogozhkin, and V. Pakholkov, "Turbulent convective heat transfer in an inclined tube filled with sodium," Technical physics 60, 1305-1309 (2015).</footnote>

<footnote n="29">²⁹Q. Wang, S.-N. Xia, W. Bo-Fu, S. De-Jun, Q. Zhou, and W. Zhen-Hua, "Flow reversals in two-dimensional thermal convection in tilted cells," Journal of Fluid Mechanics 849, 355-372 (2018).</footnote>

<footnote n="30">³⁰I. Vinogradov, L. Khezzar, and D. A. Siginer, "Heat transfer of non-Newtonian dilatant power law fluids in square and rectangular cavities," Journal of Applied Fluid Mechanics 4, 37-42 (2012).</footnote>

<footnote n="31">³¹L. Khezzar, D. A. Siginer, and I. Vinogradov, "Natural convection of power law fluids in inclined cavities," International Journal of Thermal Sciences 53, 8-17 (2012).</footnote>

<footnote n="32">³²J. Neustupa and D. A. Siginer, "Existence and structure of steady solutions of the bénard problem in a two dimensional quadrangular cavity," Nonlinear</footnote>
Analysis: Theory, Methods & Applications **123**, 68–88 (2015).

<footnote n="33"><sup>33</sup>J. Neustupa and D. A. Siginer, “Structure of the set of stationary solutions to the equations of motion of a class of generalized newtonian fluids,” Non-linear Analysis: Real World Applications **45**, 704–720 (2019).</footnote>

<footnote n="34"><sup>34</sup>X. Zheng, M. Boutaous, S. Xin, W. Cai, and D. A. Siginer, “Time-dependent oscillating viscoelastic Rayleigh-Bénard convection: viscoelastic kinetic energy budget analysis,” Physical Review Fluids , In press (2022).</footnote>

<footnote n="35"><sup>35</sup>X. Zheng, F. Hagani, M. Boutaous, R. Knikker, S. Xin, and D. A. Siginer, “Pattern selection in Rayleigh-Bénard convection with nonlinear viscoelastic fluids,” Physical Review Fluids 7, 023301 (2022).</footnote>

<footnote n="36"><sup>36</sup>H. Park and D. Ryu, “Rayleigh-Bénard convection of viscoelastic fluids in finite domains,” Journal of non-Newtonian fluid mechanics **98**, 169–184 (2001).</footnote>

<footnote n="37"><sup>37</sup>H. Park, K. Shin, and H. Sohn, “Numerical simulation of thermal convection of viscoelastic fluids using the grid-by-grid inversion method,” International journal of heat and mass transfer **52**, 4851–4861 (2009).</footnote>

<footnote n="38"><sup>38</sup>N. P. Thien and R. I. Tanner, “A new constitutive equation derived from network theory,” Journal of Non-Newtonian Fluid Mechanics **2**, 353–365 (1977).</footnote>

<footnote n="39"><sup>39</sup>D. A. Siginer, *Developments in the Flow of Complex Fluids in Tubes* (Springer, 2015).</footnote>

<footnote n="40"><sup>40</sup>J. F. Schoonen, F. H. Swartjes, G. W. Peters, F. Baaijens, and H. E. Meijer, “A 3d numerical/experimental study on a stagnation flow of a polyisobutylene solution,” Journal of non-Newtonian fluid mechanics **79**, 529–561 (1998).</footnote>

<footnote n="41"><sup>41</sup>L. M. Quinzani, R. C. Armstrong, and R. A. Brown, “Birefringence and laser-doppler velocimetry (ldv) studies of viscoelastic flow through a planar contraction,” Journal of Non-Newtonian Fluid Mechanics **52**, 1–36 (1994).</footnote>

<footnote n="42"><sup>42</sup>J. Azaiez, R. Guenette, and A. Ait-Kadi, “Entry flow calculations using multi-mode models,” Journal of Non-Newtonian Fluid Mechanics **66**, 271–281 (1996).</footnote>

<footnote n="43"><sup>43</sup>C. Carrot, J. Guillet, P. Revenu, and A. Arsac, “Experimental validation of non linear network models,” in *Rheology Series*, Vol. 5 (Elsevier, 1996) pp. 141–198.</footnote>

<footnote n="44"><sup>44</sup>R. G. Larson, *Constitutive equations for polymer melts and solutions: But-terworths series in chemical engineering* (Butterworth-Heinemann, 2013).</footnote>

<footnote n="45"><sup>45</sup>F. Hagani, M. Boutaous, R. Knikker, S. Xin, and D. A. Siginer, “Numerical modeling of non-affine viscoelastic fluid flow including viscous dissipation through a square cross-section duct: Heat transfer enhancement due to the inertia and the elastic effects,” in *ASME International Mechanical Engineering Congress and Exposition*, Vol. 84584 (American Society of Mechanical Engineers, 2020) p. V010T10A009.</footnote>

<footnote n="46"><sup>46</sup>C.-W. Shu, “Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws,” in *Advanced numerical approximation of nonlinear hyperbolic equations* (Springer, 1998) pp. 325–432.</footnote>

<footnote n="47"><sup>47</sup>R. R. Nourgaliev and T. G. Theofanous, “High-fidelity interface tracking in compressible flows: unlimited anchored adaptive level set,” Journal of Computational Physics **224**, 836–866 (2007).</footnote>

<footnote n="48"><sup>48</sup>X. Zheng, *Numerical study on viscoelastic Rayleigh-Bénard convection*, Ph.D. thesis, Université de Lyon (2021).</footnote>

<footnote n="49"><sup>49</sup>H. Park and H. Lee, “Nonlinear hydrodynamic stability of viscoelastic fluids heated from below,” Journal of non-Newtonian fluid mechanics **60**, 1–26 (1995).</footnote>

<footnote n="50"><sup>50</sup>C.-K. Li, Z.-Y. Zheng, H.-N. Zhang, F.-C. Li, S. Qian, S. W. Joo, and L. V. Kulagina, “Numerical study on secondary flows of viscoelastic fluids in straight ducts: Origin analysis and parametric effects,” Computers & Fluids **152**, 57–73 (2017).</footnote>