Onset of Thermal Instabilities in the Plane Poiseuille Flow of Weakly Elastic Fluids: Viscous Dissipation Effects

Abstract

The onset of thermal instabilities in the plane Poiseuille flow of weakly elastic fluids is examined through a linear stability analysis by taking into account the effects of viscous dissipation. The destabilizing thermal gradients may come from the different temperatures imposed on the external boundaries and/or from the volumetric heating induced by viscous dissipation. The rheological properties of the viscoelastic fluid are modeled using the Oldroyd-B constitutive equation. As in the Newtonian fluid case, the most unstable structures are found to be stationary longitudinal rolls (modes with axes aligned along the streamwise direction). For such structures, it is shown that the viscoelastic contribution to viscous dissipation may be reduced to one unique parameter: $$\gamma ={\lambda }_1(1-\mathsf{\Gamma })$$, where $${\lambda }_1$$ and $$\mathsf{\Gamma }$$ represent the relaxation time and the viscosity ratio of the viscoelastic fluid, respectively. It is found that the influence of the elasticity parameter $$\gamma$$ on the linear stability characteristics is non-monotonic. The fluid elasticity stabilizes (destabilizes) the basic Poiseuille flow if $$\gamma <{\gamma }^{*}$$ ($$\gamma >{\gamma }^{*}$$) where $${\gamma }^{*}$$ is a particular value of $$\gamma$$ that we have determined. It is also shown that when the temperature gradient imposed on the external boundaries is zero, the critical Reynolds number for the onset of such viscous dissipation/viscoelastic-induced instability may be well below the one needed to trigger the pure hydrodynamic instability in weakly elastic solutions.

1. Introduction

Hydrodynamic instabilities in the plane Poiseuille flow of Newtonian fluids have been extensively investigated in the literature. In the absence of thermal gradients, their linear and nonlinear behaviors are driven by inertia, the impact of which can be measured by a single dimensionless parameter, the Reynolds number $Re$. Orszag [1] showed that the plane Poiseuille flow becomes linearly unstable at the critical value $R{e}_c=5772$ ($R{e}_c$ based on the layer half-depth and the centerline velocity). It is well-known that the dissolution of polymer molecules in shear flows of Newtonian liquids significantly modify their stability characteristics. Sadanandan and Sureshkumar [2] theoretically investigated the viscoelastic plane Poiseuille flow of an Oldroyd-B fluid. The authors showed that departing from its Newtonian value, the critical Reynolds number $R{e}_c$ first decreases with increasing fluid elasticity until a critical value is reached, and then it increases back again. In the absence of inertia, the channel Couette and Poiseuille flows of Maxwell and Oldroyd-B fluids are linearly stable at any level of elasticity [3,4].

On the other hand, non-isothermal shear flows, such as the Poiseuille-Rayleigh-Bénard (PRB) system, have also been extensively investigated. In this case, two types of instabilities can set in (i) buoyancy-driven (thermal) instabilities, which give rise to structures in the form of thermoconvective rolls, and (ii) shear-flow (hydrodynamic) instabilities, which can lead to Tollmien-Schlichting waves at high flow rates. Thermal instabilities were studied in shear flows of Newtonian fluids [5,6], binary fluid mixtures [7,8], as well as fluid-saturated porous media [9,10]. Transient growth mechanisms have also been investigated [11]. Recently, the linear stability analyses of non-isothermal shear flows have been extended to non-Newtonian fluids [12,13,14,15]. In particular, Hirata et al. [14] investigated the linear stability of the PRB flow of an Oldroyd-B fluid. It has been shown that longitudinal rolls (LRs) are the preferred mode at the convection onset. For weakly elastic solutions, the emerging LRs are stationary, while for strongly elastic solutions, they are oscillatory.

The effects of viscous dissipation were ignored in all the aforementioned investigations. Recent studies concerning fluid and porous media [16,17,18,19,20,21,22,23,24,25,26,27] have shown that viscous dissipation, often considered a negligible effect, can induce instabilities even without external forcing. In other words, differently from the classical Rayleigh-Bénard (RB) system where the vertical temperature gradient is a consequence of the prescribed boundary conditions, viscous dissipation is an internal heat generation mechanism that can also lead to unstable temperature gradients. The buoyancy-induced instabilities caused by viscous dissipation in plane Couette flow and plane Poiseuille flow were analyzed respectively by Barletta and Nield [24] and Barletta, Celli and Nield [25]. The authors considered an adiabatic lower horizontal wall and an upper isothermal wall. With such boundary conditions, no external thermal forcing is prescribed, and the temperature gradient of the basic flow is only a consequence of the flow rate through the volumetric heating induced by viscous dissipation. Differently, Requilé et al. [27] considered fixed temperature walls to investigate how viscous dissipation changes the linear stability characteristics of PRB mixed convection flows of Newtonian fluids. Results showed that the most unstable perturbations are three-dimensional stationary LRs. The authors showed that two distinct modes of instability could occur. For moderate values of the dissipation parameter, the linear instability properties are almost similar to RB convection. When the dissipation parameter is large enough, the viscous dissipation mechanism induces a strong destabilization and triggers instability even if heating is from above. A thermal energy budget analysis is proposed to understand the physical mechanism underlying the smooth transition between the two regimes of convection.

The present paper can be viewed either as an extension of the work by Requilé et al. [27] to viscoelastic fluids, or as an extension of the study by Hirata et al. [14] to account for viscous dissipation effects. Our analysis is oriented to cover small flow rates (or small Reynolds numbers) in order to avoid the appearance of Tollmien-Schlichting modes. Furthermore, we consider the flow of diluted polymeric solutions (represented by the Oldroyd-B rheological model) presenting low levels of elasticity.

2. Governing Equations and Linear Stability Analysis

Consider the plane Poiseuille flow of an incompressible Oldroyd-B fluid in a layer of depth h in the vertical z-direction. The layer has an infinite extent in the x and y directions. The upper wall is kept at temperature $T_1$, while the lower wall is kept at temperature $T_0>T_1$, and a constant pressure gradient is imposed in the x-direction. We take $\lambda$ for the thermal conductivity, ${\rho }_0$ for the reference mass density, $C_v$ for the specific heat of the fluid at constant volume, $\mu$ for the dynamic viscosity, $\beta$ for the thermal expansion coefficient of the fluid, and g for the modulus of the gravitational acceleration. ${\lambda }_1^{*}$ and ${\lambda }_2^{*}$ are the relaxation and retardation times of the viscoelastic fluid. Assuming that the Oberbeck-Boussinesq approximation holds and taking into account viscous dissipation, the mass, momentum and energy conservation can be written as [27]:

(1)${\displaystyle u_{i,i}^{*}=0}$

(2)${\displaystyle {\rho }_0\left(\frac{\partial u_i^{*}}{\partial t^{*}}+u_j^{*}{u}_{i,j}^{*}\right)=-P_{,i}^{*}+{\tau }_{ij,j}^{*}-{\rho }_0\beta (T^{*}-T_0)g{\overrightarrow{e}}_z}$

(3)${\rho }_0{C}_v\left(\frac{\partial T^{*}}{\partial t^{*}}+u_j^{*}{T}_{,j}^{*}\right)=\lambda T_{,ii}^{*}+{\tau }_{ij}^{*}{u}_{i,j}^{*}$

The Oldroyd-B constitutive relation is:

(4)${\tau }_{ij}^{*}+{\lambda }_1^{*}\frac{\delta {\tau }_{ij}^{*}}{\delta t^{*}}=\mu \left({\gamma }_{ij}^{*}+{\lambda }_2^{*}\frac{\delta {\gamma }_{ij}^{*}}{\delta t^{*}}\right)$

where $\frac{\delta }{\delta t^{*}}$ is the Oldroyd time-derivative, ${\overrightarrow{e}}_z$ is the unit vector in the vertical direction, $u^{*}$, $T^{*}$ and $P^{*}$ are the velocity vector, temperature and hydrostatic pressure, respectively, ${\tau }_{ij}^{*}$ is the stress tensor and ${\gamma }_{ij}^{*}$ is the strain rate tensor. The last term in the energy equation ${\mathsf{\Phi }}^{*}={\tau }_{ij}^{*}{u}_{i,j}^{*}$ denotes the power per unit volume generated by viscous dissipation. The boundary conditions at the upper and lower walls are given by:

(5)$T^{*}=T_1\mathrm{at}{z}^{*}=\frac{h}{2},T^{*}=T_0\mathrm{at}{z}^{*}=-\frac{h}{2}\mathrm{and}{u}_i^{*}=0\mathrm{at}{z}^{*}=\pm \frac{h}{2}.$

We consider the following dimensionless variables: $z=z^{*}/h$, $t=\kappa t/{h^2}^{*}$, $u=h{u}^{*}/\kappa$, $P=h^2{P}^{*}/\mu {\kappa }^2$, $T=(T^{*}-T_1)\beta g{h}^3/\left(\nu \kappa \right)$, where $\kappa =\lambda /\left({\rho }_0{C}_v\right)$ is the fluid’s thermal diffusivity and $\nu =\mu /{\rho }_0$ is the kinematic viscosity. By doing so, the dimensionless governing equations take the form:

(6)${\displaystyle u_{i,i}=0}$

(7)$\left(\frac{\partial u_i}{\partial t}+u_j{u}_{i,j}\right)=-Pr{P}_{,i}+Pr{\tau }_{ij,j}+PrT{\overrightarrow{e}}_z$

(8)$\left(\frac{\partial T}{\partial t}+u_j{T}_{,j}\right)=T_{,ii}+Ge{\tau }_{ij}{u}_{i,j}$

(9)${\tau }_{ij}+{\lambda }_1\frac{\delta {\tau }_{ij}}{\delta t}={\gamma }_{ij}+{\lambda }_2\frac{\delta {\gamma }_{ij}}{\delta t}$

These equations are subjected to the following dimensionless boundary conditions:

(10)$T=0\mathrm{at}z=\frac{1}{2},T=Ra\mathrm{at}z=-\frac{1}{2}\mathrm{and}{u}_i=0\mathrm{at}z=\pm \frac{1}{2}.$

The dimensionless parameters appearing in the above equations are the Prandtl number $Pr$, the Gebhart number $Ge$, the Rayleigh number $Ra$ as well as the dimensionless relaxation ${\lambda }_1$ and retardation ${\lambda }_2$ times, defined as follows.

(11)$Pr=\frac{\nu }{\kappa }$

(12)$Ge=\frac{\beta gh}{C_v}$

(13)$Ra=\frac{\beta g{h}^3(T_0-T_1)}{\left(\nu \kappa \right)}$

(14)${\lambda }_1=\frac{{\lambda }_1^{*}\kappa }{h^2}$

(15)${\lambda }_2=\frac{{\lambda }_2^{*}\kappa }{h^2}$

System (6)–(10) presents a steady state solution that corresponds to a combination of a plane Poiseuille flow and a thermal stratification:

(16)$u_b\left(z\right)=Pef\left(z\right)$

(17)$T_b\left(z\right)=\left(\frac{1}{2}-z\right)Ra+P{e}^{2}Geg\left(z\right)$

(18)$P_b(x,z)=\frac{1}{2}Raz(1-z)-8Pex+\frac{1}{3}z\left(1-\frac{16}{5}{z}^4\right)P{e}^{2}Ge$

(19)${\tau }_{b11}={\lambda }_1(1-\mathsf{\Gamma }){\left(Pe{f}^{\prime }\left(z\right)\right)}^2$

(20)${\tau }_{b13}={\tau }_{b31}=Pe{f}^{\prime }\left(z\right)$

with:

$f\left(z\right)=(1-4z^2),g\left(z\right)=\frac{1}{3}(1-16z^4)$

Here, the Péclet number is defined as the product of the Reynolds and Prandtl numbers $Pe=RePr$, with $Re=(U_{0}h/\nu )$ and $U_0$ as the centerline velocity of the Poiseuille flow.

In the present study, the relaxation time ${\lambda }_1^{*}$ is scaled by the thermal diffusion time $h^2/\kappa$. The dimensionless parameter ${\lambda }_1$ (also known as the elasticity number) is a measure of fluid elasticity or, in other words, a measure of the fluid delayed reaction to its own advection. One may equivalently scale the relaxation time by the characteristic hydrodynamic time to yield the Deborah number $De={\lambda }_1^{*}/(h/U_0^{*})$, which is taken as a measure of normal stress effects. The elasticity and Deborah numbers are related by the formula: $De={\lambda }_{1}RePr$. Parameter $\mathsf{\Gamma }={\lambda }_2/{\lambda }_1$ is the ratio between the retardation and relaxation times of the viscoelastic fluid. Typical fluids satisfying constitutive relation (9) are solutions composed of a Newtonian solvent and a polymeric solute with viscosities ${\mu }_s$ and ${\mu }_p$, respectively. The total viscosity of the solution is given by $\mu ={\mu }_s+{\mu }_p$ and, hence, the ratio between retardation and relaxation times is $\mathsf{\Gamma }={\mu }_s/\mu$ [28]. Finally, it should be noted that ${\lambda }_1⩾{\lambda }_2$, i.e., $\mathsf{\Gamma }\in [0,1]$ and that the Maxwellian model is recovered when $\mathsf{\Gamma }=0$, whereas the Newtonian model is recovered when $\mathsf{\Gamma }=1$. Note that the base thermal gradient is the result of two contributions. The linear part comes from the different temperatures imposed on the external boundaries. The non-linear part instead comes from the viscous dissipation contribution, and vanishes in the absence of throughflow (i.e., $Pe=0$) or when viscous dissipation is neglected (i.e., $Ge=0$).

Following the classical linear stability procedure, we impose small perturbations to the base stationary solution and linearize the equations. Then, we seek solutions in the form of normal modes:

(21)$\{u,v,w,\theta ,p,{\tau }_{(i,j)}\}=\{u_n,v_n,w_n,{\theta }_n,p_n,{\tau }_{n(i,j)}\}\left(z\right)exp\left(i(k_{x}x+k_{y}y-\omega t)\right)$

where $k_x$ and $k_y$ are the wavenumbers in the x and y-directions, and $\omega$ is the oscillation frequency. The perturbation equations are then:

(22)$i{k}_x{u}_n+i{k}_y{v}_n+\mathrm{D}{w}_n=0$

(23)$\frac{1}{Pr}\left(-i\omega u_n+i{k}_{x}Pef\left(z\right)u_n+w_{n}Pe{f}^{\prime }\left(z\right)\right)=-i{k}_x{p}_n+i{k}_x{\tau }_{n11}+i{k}_y{\tau }_{n12}+\mathrm{D}{\tau }_{n13}$

(24)$\frac{1}{Pr}\left(-i\omega v_n+i{k}_{x}Pef\left(z\right)v_n\right)=-i{k}_y{p}_n+i{k}_x{\tau }_{n21}+i{k}_y{\tau }_{n22}+\mathrm{D}{\tau }_{n23}$

(25)$\frac{1}{Pr}\left(-i\omega w_n+i{k}_{x}Pef\left(z\right)w_n\right)=-\mathrm{D}{p}_n+{\theta }_n+i{k}_x{\tau }_n^{31}+i{k}_y{\tau }_{n32}+\mathrm{D}{\tau }_{n33}$

(26)$\begin{array}{c}\hfill -i\omega {\theta }_n+i{k}_{x}Pef\left(z\right){\theta }_n-Ra{w}_n+P{e}^{2}Ge{g}^{\prime }\left(z\right)w_n-({\mathrm{D}}^2-k^2){\theta }_n\\ \hfill -Ge\left[Pe{f}^{\prime }\left(z\right){\tau }_{n13}+i{k}_x{\tau }_{b11}{u}_n+{\tau }_{b13}\left(i{k}_x{w}_n+\mathrm{D}{u}_n\right)\right]=0\end{array}$

where $k^2=k_x^2+k_y^2$, $\mathrm{D}=\partial /\partial z$ and ${\mathrm{D}}^2={\partial }^2/\partial z^2$. The boundary conditions become:

(27)${\theta }_n=0\mathrm{at}z=\pm \frac{1}{2}\mathrm{and}{u}_n=v_n=w_n=\mathrm{D}{w}_n=0\mathrm{at}z=\pm \frac{1}{2}.$

Equations (22)–(26) with boundary conditions (27) form a generalized eigenvalue problem involving, in general, the six dimensionless parameters: ${\lambda }_1$, $\mathsf{\Gamma }$, $Ra$, $Pe$, $Pr$ and $Ge$. This eigenvalue problem was solved numerically by considering disturbances in the form of longitudinal rolls ($k_x=0$, $k_y=k$), transverse rolls ($k_x=k$, $k_y=0$) and oblique rolls ($k_x\ne 0$, $k_y\ne 0$) by means of a shooting method (please refer to [14,29] for details on the numerical method).

2.1. Transverse Rolls

For TRs, whose axes are perpendicular to the main flow direction, we pose $k_y=0$ and $k_x=k$. In this case, we use the continuity Equation (22) to eliminate $u_n$ and reduce the problem to the following two-dimensional system of equations:

(28)$\begin{array}{c}\hfill \frac{i}{Pr}\left(\right(k_{x}RePrf\left(z\right)-i\omega )\left({\mathrm{D}}^2-k^2\right)-k_{x}RePrf{}^{″}\left(z\right))w_n+k^2\mathrm{D}{\tau }_{n33}\\ \hfill +\mathrm{D}\mathrm{L}\left(w_n\right)+\left({\mathrm{D}}^2+k^2\right)\mathrm{J}\left(w_n\right)+k^2{\theta }_n=0\end{array}$

(29)$\begin{array}{c}\hfill i\omega {\theta }_n+i{k}_{x}Pef\left(z\right){\theta }_n-Ra{w}_n+\mathsf{\Lambda }{g}^{\prime }\left(z\right)w_n-({\mathrm{D}}^2-k^2){\theta }_n\\ \hfill -GePe\left[f^{\prime }\left(z\right){\tau }_{n13}-{\lambda }_1(1-\mathsf{\Gamma })Pe{\left(f^{\prime }\left(z\right)\right)}^{2}D{w}_n+f^{\prime }\left(z\right)\left(i{k}_x{w}_n-\frac{1}{i{k}_x}{D}^2{w}_n\right)\right]=0\end{array}$

with:

(30)$\begin{array}{ccc}\hfill \mathrm{L}\left(w_n\left(z\right)\right)& =& -k^2{\tau }_{n11}\left(z\right)\hfill \\ & =& \left[2ikPe{\lambda }_{1}f{}^{\prime }\left(z\right)\mathrm{J}\right(w_n\left(z\right))+(2ikPe{\lambda }_1(2\mathsf{\Gamma }-1)f{}^{\prime }\left(z\right){\mathrm{D}}^2\hfill \\ & +& 4k^{2}P{e}^2{\lambda }_1^2(1-\mathsf{\Gamma })(f{}^{\prime }\left(z\right){)}^2+2k^2+2i{k}^{2}kPe{\lambda }_1\mathsf{\Gamma }f\left(z\right)\hfill \\ & -& 2i{k}^2\omega {\lambda }_1\mathsf{\Gamma }\mathrm{D}+4k^{2}P{e}^2{\lambda }_1^2(1-\mathsf{\Gamma })f{}^{\prime }\left(z\right)f{}^{″}\left(z\right)\hfill \\ & +& 2i{k}^{3}Pe{\lambda }_1\mathsf{\Gamma }f{}^{\prime }\left(z\right)\left)w_n\right]/[1+ikPe{\lambda }_{1}f\left(z\right)-i\omega {\lambda }_1]\hfill \end{array}$

(31)$\begin{array}{ccc}\hfill \mathrm{J}\left(w_n\left(z\right)\right)& =& ik{\tau }_{n13}\left(z\right)\hfill \\ & =& [ikPe{\lambda }_{1}f{}^{\prime }\left(z\right){\tau }_{n33)}\left(z\right)-(2ikPe{\lambda }_1\mathsf{\Gamma }f{}^{\prime }\left(z\right)\mathrm{D}-i\omega {\lambda }_1\mathsf{\Gamma }{\mathrm{D}}^2\hfill \\ & +& (1+k^2)(1+ikPe{\lambda }_1\mathsf{\Gamma }f\left(z\right))+ikPe{\lambda }_1(1-\mathsf{\Gamma })f{}^{″}\left(z\right)-i\omega k^2{\lambda }_1\mathsf{\Gamma }\hfill \\ & +& 2k^{2}P{e}^2{\lambda }_1^2(1-\mathsf{\Gamma })\left(f{}^{\prime }\left(z\right){)}^2\right)w_n]/[1+ikPe{\lambda }_{1}f\left(z\right)-i\omega {\lambda }_1]\hfill \end{array}$

(32)$\begin{array}{ccc}\hfill {\tau }_{n33}\left(z\right)& =& \frac{2(\mathrm{D}{w}_n(i+\mathsf{\Gamma }{\lambda }_1\omega -k\mathsf{\Gamma }{\lambda }_1{u}_b)+k{\lambda }_1{w}_n(\mathsf{\Gamma }\mathrm{D}{u}_b-{\tau }_{b31}))}{i+{\lambda }_1\omega -k{\lambda }_1{u}_b}\hfill \end{array}$

(33)$\begin{array}{ccc}\hfill {\tau }_{n13}\left(z\right)& =& \frac{1}{i+{\lambda }_1\omega -k{\lambda }_1{u}_b}[\mathrm{D}{u}_n(i+\mathsf{\Gamma }{\lambda }_1\omega -k\mathsf{\Gamma }{\lambda }_1{u}_b)\hfill \\ & -& i(-{\lambda }_1(\mathrm{D}{w}_n+ik{u}_n){\tau }_{b31}+{\lambda }_1\mathrm{D}{u}_b(3\mathsf{\Gamma }\mathrm{D}{w}_n+ik\mathsf{\Gamma }{u}_n-{\tau }_{n33})\hfill \\ & +& w_n(-k(i+\mathsf{\Gamma }{\lambda }_1\omega )+{\lambda }_1(k^2\mathsf{\Gamma }{u}_b-ik{\tau }_{b11}+\mathrm{D}{\tau }_{b31}-\mathsf{\Gamma }{\mathrm{D}}^2{u}_b)\left)\right]\hfill \end{array}$

(34)$\begin{array}{ccc}\hfill {\tau }_{n11}\left(z\right)& =& \frac{-i}{I+{\lambda }_1\omega -k{\lambda }_1{u}_b}[2k{u}_n(-i-\mathsf{\Gamma }{\lambda }_1\omega +k\mathsf{\Gamma }{\lambda }_1{u}_b-i{\lambda }_1{\tau }_{b11})\hfill \\ & -& 2{\lambda }_1\mathrm{D}{u}_n{\tau }_{b31}+2{\lambda }_1\mathrm{D}{u}_b(2\mathsf{\Gamma }\mathrm{D}{u}_n+ik\mathsf{\Gamma }{w}_n-{\tau }_{n31})+{\lambda }_1{w}_n\mathrm{D}{\tau }_{b11}]\hfill \end{array}$

2.2. Longitudinal Rolls

For LRs, which are helicoidal rolls with axes parallel to the main direction, one has $k_y=k$ and $k_x=0$. Equations (22)–(26) can then be reduced to:

(35)$\begin{array}{c}\hfill \frac{L_1}{Pr}\left(-i\omega u_n+w_{n}Pe{f}^{\prime }\left(z\right)\right)-L_2\left({\mathrm{D}}^2-k^2\right)u_n\\ \hfill -{\lambda }_{1}Pe\left(\frac{L_2}{L_1}-\mathsf{\Gamma }\right)\left(f{}^{\prime }\left(z\right)\left({\mathrm{D}}^2-k^2\right)w_n+2f{}^{″}\left(z\right)\mathrm{D}{w}_n\right)=0\end{array}$

(36)$\begin{array}{c}\hfill L_1\left(\frac{-i\omega }{Pr}\left({\mathrm{D}}^2-k^2\right)w_n+k^2{\theta }_n\right)-L_2{\left({\mathrm{D}}^2-k^2\right)}^2{w}_n=0\end{array}$

(37)$\begin{array}{c}\hfill L_1\left(-i\omega {\theta }_n-Ra{w}_n+\mathsf{\Lambda }{g}^{\prime }\left(z\right)w_n-({\mathrm{D}}^2-k^2){\theta }_n\right)=\\ \hfill \mathsf{\Lambda }{f}^{\prime }\left(z\right)\left[{\lambda }_{1}f{}^{\prime }\left(z\right)(1+2\frac{L_2}{L_1}-3\mathsf{\Gamma })\mathrm{D}{w}_n-{\lambda }_1(1-\mathsf{\Gamma })f{}^{″}\left(z\right)w_n+\frac{(1+L_2)}{Pe}\mathrm{D}{u}_n\right]\end{array}$

with:

(38)$\begin{array}{c}\hfill L_1=1-i\omega {\lambda }_1{L}_2=1-i\omega {\lambda }_1\mathsf{\Gamma }\end{array}$

Note that we have introduced the parameter $\mathsf{\Lambda }=GeP{e}^2$ for convenience. When viscous dissipation is not accounted for ($\mathsf{\Lambda }=0$), the above problem reduces to two Equations (36) and (37) for $w_n$ and ${\theta }_n$, and the velocity component $u_n$ does not interfere in the calculation of the critical values. When $\mathsf{\Lambda }\ne 0$, it is possible to reduce the viscoelastic contribution to viscous dissipation by one unique parameter $\gamma ={\lambda }_1(1-\mathsf{\Gamma })$. As in [27], we introduce the wall-normal voticity ${\eta }_n$, defined by ${\eta }_n=ik{u}_n=Re{\tilde{\eta }}_n$. For stationary LRs ($\omega =0$), the system (35)–(37) then takes the form:

(39)$\begin{array}{c}\hfill \left({\mathrm{D}}^2-k^2\right){\tilde{\eta }}_n=ik{w}_n{f}^{\prime }\left(z\right)-ik\gamma Pr\left(f{}^{\prime }\left(z\right)\left({\mathrm{D}}^2-k^2\right)w_n+2f{}^{″}\left(z\right)\mathrm{D}{w}_n\right)\end{array}$

(40)$\begin{array}{c}\hfill k^2{\theta }_n-{\left({\mathrm{D}}^2-k^2\right)}^2{w}_n=0\end{array}$

(41)$\begin{array}{c}\hfill -Ra{w}_n+\mathsf{\Lambda }{g}^{\prime }\left(z\right)w_n-({\mathrm{D}}^2-k^2){\theta }_n-\mathsf{\Lambda }{f}^{\prime }\left(z\right)\left[3\gamma f^{\prime }\left(z\right)\mathrm{D}{w}_n-w_n\gamma f^{″}\left(z\right)-\frac{2i}{kPr}D{\tilde{\eta }}_n\right]=0\end{array}$

We note that the last term in the right hand side of Equation (39) represents a fluid elasticity contribution to the wall-normal vorticity and may induce qualitative changes in the linear characteristics of the instability.

3. Discussion

Before discussing the effects of viscous heating, it is important to recall the results by Li and Khayat [30] on the linear stability of the Rayleigh-Bénard (RB) convection of viscoelastic fluids. In the absence of throughflow (i.e., $Re=0$), two different convection regimes can be distinguished: the weakly elastic regime corresponds to a parameter set where the conductive base state bifurcates to a stationary convective state; in the strongly elastic regime, however, the onset of convection is oscillatory. Hirata et al. [14] extended the previous analysis by adding a Poiseuille flow to the problem (PRB configuration). It has been shown that stationary/oscillatory LRs are the selected patterns at the convection onset for weakly and strongly elastic solutions, respectively. Their critical thresholds are not affected by the parallel flow, and therefore, their stability characteristics are the same as those obtained without throughflow. On the other hand, TRs are stabilized by increasing the main flow intensity.

3.1. Effects of Viscous Dissipation and Viscoelasticity on Critical Parameters

We shall now study the effects of viscous dissipation on the thresholds of convective instability of viscoelastic fluids in the PRB configuration. As a first step, it is convenient to determine what type of structure will prevail when viscous dissipation is considered. Figure 1 illustrates the evolution of the critical Rayleigh number $R{a}_c$ with the Reynolds number $Re$ for LRs and TRs in weakly and strongly elastic solutions. At $Re=0$, the thresholds for LRs and TRs are the same. As expected, the critical Rayleigh number for strongly elastic solutions is lower than the one for weakly elastic solutions (which is the same as the one of the RB problem: $R{a}_c=1707.76$). The critical value $R{a}_c$ for TRs significantly increases with $Re$, either in the weakly elastic regime or in the strongly elastic one. On the other hand, the critical value $R{a}_c$ for LRs shows a very mild decrease with $Re$. It has been verified that the thresholds for ORs always lie in between the ones for LRs and TRs. Hence, we conclude that LRs are the selected pattern at the convection onset, either in the weakly elastic or in the strongly elastic regime. In the weakly elastic regime, the critical frequency ${\omega }_c$ of LRs is found to be zero for all values of $\mathsf{\Lambda }$ considered.

From this point on, we focus our results on the stability characteristics of stationary LRs in the weakly elastic regime. Figure 2a presents the evolution of the critical Rayleigh number $R{a}_c$ with the dissipation parameter $\mathsf{\Lambda }$ for stationary LRs and different values of the fluid elasticity $\gamma$. The Prandtl number is kept fixed to $Pr=450$ (please refer to [27] for a discussion on the effect of Prandtl number for realistic values of $Ge$ and $\mathsf{\Lambda }$). Parameters ${\lambda }_1$ and $\mathsf{\Gamma }$ have been chosen so that the starting regime (without throughflow) is always weakly elastic. Figure 2a shows that for all $\gamma$, the critical Rayleigh number decreases regularly with $\mathsf{\Lambda }$, demonstrating the destabilizing effects of viscous dissipation. It is interesting to note in Figure 2a that the critical Rayleigh number decreases with increasing dissipation until it crosses the x-axis and becomes negative. This means that the combined effect of viscous dissipation and the fluid elasticity may induce instability for a zero external thermal gradient ($R{a}_c=0$) or even for a stabilizing thermal gradient imposed on the external walls ($R{a}_c<0$). In these cases, the instability is triggered by the internal heating due to viscous dissipation.

On the other hand, the elasticity effect on the critical Rayleigh number for a fixed value of $\mathsf{\Lambda }$ may be stabilizing or destabilizing depending on the viscoelestic parameter $\gamma$. Indeed, Figure 2a shows that for any assigned value of $\mathsf{\Lambda }$ and for $\gamma =0.00625$, the critical Rayleigh number is always higher than the one corresponding to a Newtonian fluid, while the opposite is observed for $\gamma =0.0125$ and $\gamma =0.02125$. The nonmonotonic dependence on the critical Rayleigh number on $\gamma$ is also observed for the critical wave number. As can be seen from Figure 2b, the critical wave number curve is located below the Newtonian one for $\gamma =0.00625$, while this curve is above the Newtonian one for $\gamma =0.0125$ and $\gamma =0.02125$.

The nonmonotonic influence of the elasticity parameter $\gamma$ on the linear stability characteristics reported above is now studied in detail for a fixed value of $\mathsf{\Lambda }$. Figure 3a,b present respectively the dependence of the critical Rayleigh number $R{a}_c$ and the critical wave number $k_c$ on the viscoelastic parameter $\gamma$ for the dissipation parameter $\mathsf{\Lambda }=7500$. It can be seen from Figure 3a that as $\gamma$ increases, the critical Rayleigh number increases monotonically up to its maximum value $R{a}_c^{Max}$ at $\gamma ={\gamma }^{*}\approx 0.005$. It is found that $R{a}_c^{Max}$ is similar to the critical Rayleigh number for the onset of three-dimensional longitudinal rolls in the Poiseuille-Rayleigh-Bénard (PRB) problem in the absence of viscous dissipation $R{a}_c^{Max}\simeq R{a}_c^{RB}=1707.76$. This is accompanied by a monotonic decrease in the critical wave number up to $k_c\simeq k_c^{RB}=3.12$, as shown in Figure 3b. When $\gamma$ exceeds ${\gamma }^{*}$, the critical Rayleigh number decreases monotonically while the critical wave number increases. We then conclude that the viscoelastic parameter $\gamma$ induces a stabilizing effect if $\gamma <{\gamma }^{*}$ and a destabilizing effect otherwise. It is instructive to note that for $\gamma ={\gamma }^{*}$, the stabilizing effect of the elasticity cancels the destabilizing effect of the viscous dissipation and one recovers the linear characteristics of the PRB problem without viscous dissipation. Using the above results, it is possible to estimate the parameter set that locates the boundary separating stabilizing and destabilizing elastic regimes. For fixed values of the viscosity ratio $\mathsf{\Gamma }$, there exists a value of the relaxation time ${\lambda }_1^{*}$, coresponding to ${\gamma }^{*}$ where the critical Rayleigh number reaches its maximum value, namely, ${\lambda }_1^{*}\approx 0.005/(1-\mathsf{\Gamma })$. The elasticity is therefore, stabilizing (destabilizing) if ${\lambda }_1<{\lambda }_1^{*}$ (${\lambda }_1>{\lambda }_1^{*}$).

3.2. Spatial Distribution of the Eigen-Fields

The dependence of the eigenfunctions on the vertical coordinate is shown in Figure 4 with the dissipation parameter $\mathsf{\Lambda }=7500$ both for Newtonian fluids and viscoelastic fluids with different values of the parameter $\gamma$. In order to have a qualitative understanding of the elasticity effect, we fix the parameter $\gamma$ to three values: $\gamma =0.0025$ in the region of stabilizing effect of the elasticity, $\gamma =0.005$ for which the linear stability properties are found to be the same as for the PRB problem and $\gamma =0.02125$ in the region of the destabilizing elasticity. Computations indicate that the wall-normal vorticity and the stream function are imaginary valued while the temperature deviation is a pure real value. As shown in Figure 4 by continuous black lines for Newtonian fluids (i.e., $\gamma =0$), the wall-normal vorticity, the temperature deviation and the stream function have asymmetric profiles with respect to the midplane. We remark that the maximum of the three eigenfunctions is located in the upper part of the cell. The observed asymmetric profile is a consequence of the presence of viscous dissipation. When elasticity is taken into account, drastic changes occur in the vertical distribution of the eigenfunctions. Figure 4 shows that for the special value $\gamma =0.005$, where the stabilizing effect of the elasticity reaches its maximum, the system tends to recover a midplane reflection symmetry. In this case, the asymmetric tendency due to the destabilizing effect of viscous dissipation is counter-balanced by the stabilizing contribution of the fluid elasticity and the emerging convection profile is qualitatively similar to the one observed in the classical PRB problem. Increasing $\gamma$ to $\gamma =0.02125$, both the elasticity and the viscous dissipation destabilizes the basic Poiseuille flow. The mirror symmetry observed for $\gamma =0.005$ is lost and the absolute value of the wall-normal vorticity increases considerably as shown by the red-dashed-dotted line of Figure 4a. The maximum of the absolute value of the wall-normal vorticity is reached in the lower part of the cell. As a consequence, the maximum value of the stream function and the temperature deviation are also located in the lower part of the medium as can be seen by the red-dashed-dotted lines of Figure 4b,c, respectively. Moreover, we observe that the temperature deviation changes sign in the upper part of the cell.

All the above results may be well visualized by representing the streamlines and the isotherms of the stationary longitudinal rolls at the onset of convection. We focus on the qualitative change of such structures as a result of the coupling between the dissipation and the viscoelasticity of the fluid. Computed streamlines and isotherms are presented for $\mathsf{\Lambda }=7500$ and $Pr=450$ in Figure 5, Figure 6 and Figure 7. Figure 5 represents the Newtonian case, while Figure 6 and Figure 7 correspond to the viscoelastic fluids with $\gamma =0.005$ and $\gamma =0.02125$, respectively. For comparison purposes, the horizontal y length in all of the figures is kept equal to the wavelength of the longitudinal rolls for Newtonian fluids. As the critical wave number for viscoelastic fluids with $\gamma =0.005$ ($\gamma =0.02125$) is smaller (larger) than the one associated to the Newtonian fluids, the emerging longitudinal rolls have a larger (smaller) cross section as seen from Figure 6 and Figure 7. Further qualitative assessments are referred from Figure 5, Figure 6 and Figure 7. One can note that while convection develops almost symmetrically with respect to the horizontal mid-plane for stabilizing fluid elasticity with $\gamma =0.005$, the core of the convective cells moves upwards for Newtonian fluids and downwards for destabilizing fluid elasticity with $\gamma =0.02125$. In addition, Figure 7b shows a change of sign of the isotherms, implying a vertical alternation between hot and cold currents.

3.3. Relation to Hydrodynamic Instability

Sureshkumar and Beris [31] as well as Sadanan and Sureshkumar [2] theoretically studied the stability of the same problem in the absence of viscous heating (i.e., the isothermal plane Poiseuille flow of Oldroyd-B fluids). As in the present work, these authors took small values of the elasticity parameter ${\lambda }_1⩽0.0035$ and $0⩽\mathsf{\Gamma }⩽1$, characterizing fluids in the weakly elastic regime. The thresholds found for the appearance of the pure hydrodynamic instability were at least $R{e}_c\sim \mathcal{O}\left({10}^3\right)$ (Reynolds number based on the total layer depth and centerline velocity), which is much larger than the values displayed in Figure 1. This means that the dissipation-induced thermal instability may precede the pure hydrodynamic instability for weakly elastic fluids. In the case $Ra=0$ (where the upper and lower walls are kept at the same temperature), the parameters $\mathsf{\Lambda }=GeP{e}^2$ and $Re=Pe/Pr$ are related by the following expression:

(42)$Re=\frac{1}{Pr}\sqrt{\frac{\mathsf{\Lambda }}{Ge}}$

We take $Re(\gamma ,Pr,Ge)$ as the control parameter in order to be able to study the case $Ra=0$ and compare our results to previous works. Figure 8a displays the influence of parameter $\gamma$ on the critical Reynolds number for a weakly elastic fluid, with ${\lambda }_1=0.025$, $Pr=450$ and a fixed reasonable value of the Gebhart number $Ge={10}^{-4}$ to account for viscous dissipation effects. The critical Newtonian value of the Reynolds number is $R{e}_c^N=36.68$. It can be seen that parameter $\gamma$ has a stabilizing effect ($R{e}_c>R{e}_c^N$) for $\gamma <0.01$, otherwise it plays a destabilizing role on the system ($R{e}_c<R{e}_c^N$). The critical Reynolds number shows a peak near $\gamma \sim 0.005$, followed by an abrupt decrease and a gentler slope after $\gamma \sim 0.02$. This same behavior is observed in the wave number plot of Figure 8b: it grows quickly untill $\gamma \sim 0.005$, and after an abrupt decrease it shows a mild increase with $\gamma$.

4. Conclusions

In the present paper, we investigate the effects of viscous dissipation on linear instability of viscoelastic fluids subjected simultaneously to a basic Poiseuille flow and a vertical temperature gradient. In the framework of the Oldroyd-B model, we find that the base state is composed of a non-vanishing normal stress component induced by viscoelasticity, in addition to the shear stress component induced by the Poiseuille flow. The linear stability of the base state is investigated numerically using a shooting method. The results indicate that stationary longitudinal rolls aligned with the direction of the Poiseuille flow represent the preferred mode for the onset of convection for weakly viscoelastic fluids. In addition to Prandtl number, the critical Rayleigh number and wave number of stationary longitudinal rolls depend on a dimensionless parameter $\mathsf{\Lambda }$, which measures the strength of viscous dissipation, and a unique parameter $\gamma ={\lambda }_1(1-\mathsf{\Gamma })$, where ${\lambda }_1$ and $\mathsf{\Gamma }$ represent the relaxation time and the viscosity ratio of the viscoelastic fluid, respectively. It is found that the viscous dissipation is always destabilizing, while the fluid elasticity induces a stabilization provided that $\gamma$ remains below a particular value ${\gamma }^{*}$. However, when $\gamma$ exceeds ${\gamma }^{*}$, a strong destabilization is observed. For $\gamma ={\gamma }^{*}$, we showed that the competition between the fluid elasticity and the viscous dissipation leads to the appearance of longitudinal rolls with linear characteristics similar to those associated with the Poiseuille-Rayleigh-Bénard problem without viscous dissipation. The influence of viscous dissipation and fluid elasticity on the structural changes in the perturbation eigen-fields is also analyzed. In particular, we found that the elasticity increases significantly at the intensity of the wall-normal vorticity. Finally, when the temperature gradient imposed on the external boundaries is zero, the combined effect of viscous dissipation and fluid elasticity may also trigger the instability. We found that the critical Reynolds number needed to develop this thermal instability may be significantly lower than the critical Reynolds number associated the pure hydrodynamic instability.

Author Contributions

Conceptualization, methodology, software, validation, formal analysis, investigation, resources, data curation, writing, visualization, supervision, funding acquisition, S.C.H. and M.N.O. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Unité de Mécanique de Lille, URL 7512, Université de Lille.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

$$C_v$$	specific heat of the fluid at constant volume
g	modulus of the gravitational acceleration
$$Ge$$	Gebhart number, Equation (12)
h	fluid layer depth
k	dimensionless wave number, $$k={(k_x^2+k_y^2)}^{1/2}$$
$$(k_x,k_y)$$	dimensionless wave vector
$$Pe$$	Peclet number, $$Pe=RePr$$
$$Pr$$	Prandtl number, Equation (11)
$$Ra$$	Rayleigh number, Equation (13)
$$Re$$	Reynolds number, $$Re=U_{0}h/\nu$$
$$U_0$$	centerline velocity of the Poiseuille flow
Greek symbols
$$\beta$$	thermal expansion coefficient
$$\mathsf{\Gamma }$$	viscosity ratio, $$\mathsf{\Gamma }={\lambda }_2/{\lambda }_1={\mu }_s/\mu$$
$$\gamma$$	viscoelastic parameter, $$\gamma ={\lambda }_1(1-\mathsf{\Gamma })$$
$$\kappa$$	thermal diffusivity
$$\lambda$$	thermal conductivity
$${\lambda }_1$$	dimensionless relaxation time, Equation (14)
$${\lambda }_2$$	dimensionless retardation time, Equation (15)
$$\mathsf{\Lambda }$$	dissipation parameter, $$\mathsf{\Lambda }=GeP{e}^2$$
$$\mu$$	total dynamic viscosity, $$\mu ={\mu }_s+{\mu }_p$$
$${\mu }_s$$	solvent dynamic viscosity
$${\mu }_p$$	polymeric solute dynamic viscosity
$$\nu$$	total kinematic viscosity, $$\nu =\mu /{\rho }_0$$
$$\omega$$	dimensionless oscillation frequency
$${\rho }_0$$	reference mass density

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. Influence of the Reynolds number on the stability of transverse (TRs) and longitudinal (LRs) rolls for a weakly ([Forumla omitted. See PDF.]) and a strongly ([Forumla omitted. See PDF.]) elastic fluid with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 2. (a) Critical Rayleigh number and (b) corresponding wave number as functions of the dissipation parameter [Forumla omitted. See PDF.] for stationary LRs in the weakly elastic regime with [Forumla omitted. See PDF.], and various [Forumla omitted. See PDF.]. The dotted lines correspond to the Newtonian fluid case.

Figure 3. (a) Critical Rayleigh number and (b) corresponding wave number as functions of the viscoelastic parameter [Forumla omitted. See PDF.], with [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]. The dashed lines correspond to the Newtonian fluid case.

Figure 4. Dependence of the eigenfunctions on the vertical coordinate: (a) the wall-normal vorticity, (b) the stream function and (c) the temperature deviation with [Forumla omitted. See PDF.], and [Forumla omitted. See PDF.]. The Newtonian fluid is represented in the inset of (a) and by continuous black lines in (b,c). The viscoelastic fluids are depicted with [Forumla omitted. See PDF.] (yellow-dashed), [Forumla omitted. See PDF.] (green-dotted) and [Forumla omitted. See PDF.] (red-dashed-dotted).

Figure 5. (a) Streamfunction and (b) isotherms for a Newtonian fluid ([Forumla omitted. See PDF.]) with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 6. (a) Streamfunction and (b) isotherms for a weakly elastic fluid ([Forumla omitted. See PDF.]) with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 7. (a) Streamfunction and (b) isotherms for a weakly elastic fluid ([Forumla omitted. See PDF.]) with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Figure 8. (a) Critical Reynolds number and (b) corresponding wave number as functions of parameter [Forumla omitted. See PDF.] for LRs in the weakly elastic regime, with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

References

1. Orszag, S.A. Accurate solution of the Orr-Sommerfeld stability equation. J. Fluid Mech.; 1971; 50, pp. 689-703. [DOI: https://dx.doi.org/10.1017/S0022112071002842]

2. Sadanandan, B.; Sureshkumar, R. Viscoelastic effects on the stability of wall-bounded shear flows. Phys. Fluids; 2002; 14, pp. 41-48. [DOI: https://dx.doi.org/10.1063/1.1425847]

3. Denn, M.M. Issues in Viscoelastic Fluid Mechanics. Annu. Rev. Fluid Mech.; 1990; 22, pp. 13-32. [DOI: https://dx.doi.org/10.1146/annurev.fl.22.010190.000305]

4. Larson, R.G. Instabilities in viscoelastic flows. Rheol. Acta; 1992; 31, pp. 213-263. [DOI: https://dx.doi.org/10.1007/BF00366504]

5. Muller, H.W.; Lucke, M.; Kamps, M. Transversal convection patterns in horizontal shear flow. Phys. Rev. A; 1992; 45, pp. 3714-3726. [DOI: https://dx.doi.org/10.1103/PhysRevA.45.3714] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/9907419]

6. Carrière, P.; Monkewitz, P.A. Convective versus absolute instability in mixed Rayleigh-Bénard-Poiseuille convection. J. Fluid Mech.; 1999; 384, pp. 243-262. [DOI: https://dx.doi.org/10.1017/S0022112098004145]

7. Büchel, P.; Lücke, M. Influence of through-flow on binary fluid convection. Phys. Rev. E; 2000; 61, 3793. [DOI: https://dx.doi.org/10.1103/PhysRevE.61.3793]

8. Hu, J.; Ben Hadid, H.; Henry, D. Linear stability analysis of Poiseuille-Rayleigh-Bénard flows in binary fluids with Soret effect. Phys. Fluids; 2007; 19, pp. 034101(1)-034101(17). [DOI: https://dx.doi.org/10.1063/1.2709931]

9. Delache, A.; Ouarzazi, M.N.; Combarnous, M. Spatio-temporal stability analysis of mixed convection flows in porous media heated from below: Comparison with experiments. Int. J. Heat Mass Transf.; 2007; 50, pp. 1485-1499. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2006.08.040]

10. Ouarzazi, M.N.; Mejni, F.; Delache, A.; Labrosse, G. Nonlinear global modes in inhomogeneous mixed convection flows in porous media. J. Fluid Mech.; 2008; 595, pp. 367-377. [DOI: https://dx.doi.org/10.1017/S0022112007009159]

11. Jerome, J.J.S.; Chomaz, J.M.; Huerre, P. Transient growth in Rayleigh-Benard-Poiseuille/Couette convection. Phys. Fluids; 2012; 24, 044103. [DOI: https://dx.doi.org/10.1063/1.4704642]

12. Nouar, C.; Benaouda-Zouaoui, B.; Desaubry, C. Laminar mixed convection in a horizontal annular duct. case of thermodependent non-Newtonian fluid. Eur. J. Mech. B Fluids; 2000; 19, pp. 423-452. [DOI: https://dx.doi.org/10.1016/S0997-7546(00)00120-5]

13. Métivier, C.; Nouar, C. Linear stability of Rayleigh-Bénard Poiseuille flow for thermodependent viscoplastic fluids. J. Non-Newton. Fluid Mech.; 2009; 163, pp. 1-8. [DOI: https://dx.doi.org/10.1016/j.jnnfm.2009.06.001]

14. Hirata, S.; Alves, L.S.d.B.; Delenda, N.; Ouarzazi, M. Convective and absolute instabilities in Rayleigh-Bénard-Poiseuille mixed convection for viscoelastic fluids. J. Fluid Mech.; 2015; 765, pp. 167-210. [DOI: https://dx.doi.org/10.1017/jfm.2014.721]

15. Alves, L.S.d.B.; Hirata, S.; Ouarzazi, M. Linear onset of convective instability for Rayleigh-Bénard-Couette flows of viscoelastic fluids. J. Non-Newton. Fluid Mech.; 2016; 231, pp. 79-90. [DOI: https://dx.doi.org/10.1016/j.jnnfm.2016.03.007]

16. Roy, K.; Ponalagusamy, R.; Murthy, P. The effects of double-diffusion and viscous dissipation on the oscillatory convection in a viscoelastic fluid saturated porous layer. Phys. Fluids; 2020; 32, 094108. [DOI: https://dx.doi.org/10.1063/5.0020076]

17. Straughan, B. Nonlinear Stability for Thermal Convection in a Brinkman Porous Material with Viscous Dissipation. Transp. Porous Media; 2020; 134, pp. 303-314. [DOI: https://dx.doi.org/10.1007/s11242-020-01446-5]

18. Norouzi, M.; Dorrani, S.; Shokri, H.; Anwar Beg, O. Effects of viscous dissipation on miscible thermo-viscous fingering instability in porous media. Int. J. Heat Mass Transf.; 2019; 129, pp. 212-223. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2018.09.048]

19. Miklavcic, M. Bistable fully developed mixed convection flow with viscous dissipation in a vertical channel. R. Soc. Open Sci.; 2018; 5, 171880. [DOI: https://dx.doi.org/10.1098/rsos.171880] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30109046]

20. Roy, K.; Murthy, P. Linear Stability of the Double-Diffusive Convection in a Horizontal Porous Layer with Open Top: Soret and Viscous Dissipation Effects. Transp. Porous Media; 2018; 122, pp. 693-712. [DOI: https://dx.doi.org/10.1007/s11242-018-1022-5]

21. Celli, M.; Kuznetsov, A. A new mechanism for buoyancy driven convection in pulsating viscous flows: A theoretical study. Int. J. Heat Mass Transf.; 2018; 118, pp. 340-348. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.10.112]

22. Celli, M.; Barletta, A. Onset of convection in a non-Newtonian viscous flow through a horizontal porous channel. Int. J. Heat Mass Transf.; 2018; 117, pp. 1322-1330. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2017.10.089]

23. Barletta, A. On the thermal instability induced by viscous dissipation. Int. J. Therm. Sci.; 2015; 88, pp. 238-247. [DOI: https://dx.doi.org/10.1016/j.ijthermalsci.2014.02.009]

24. Barletta, A.; Nield, D.A. Convection-dissipation instability in the horizontal plane Couette flow of a highly viscous fluid. J. Fluid Mech.; 2010; 662, pp. 475-492. [DOI: https://dx.doi.org/10.1017/S0022112010003320]

25. Barletta, A.; Celli, M.; Nield, D.A. On the onset of dissipation thermal instability for the Poiseuille flow of a highly viscous fluid in a horizontal channel. J. Fluid Mech.; 2011; 681, pp. 499-514. [DOI: https://dx.doi.org/10.1017/jfm.2011.213]

26. Requilé, Y.; Hirata, S.C.; Ouarzazi, M.N.; Barletta, A. Weakly nonlinear analysis of viscous dissipation thermal instability in plane Poiseuille and plane Couette flows. J. Fluid Mech.; 2020; 886, A26. [DOI: https://dx.doi.org/10.1017/jfm.2019.1062]

27. Requilé, Y.; Hirata, S.C.; Ouarzazi, M.N. Viscous dissipation effects on the linear stability of Rayleigh-Bénard-Poiseuille/Couette convection. Int. J. Heat Mass Transf.; 2020; 146, 118834. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2019.118834]

28. Bird, R.B.; Armstrong, R.C.; Hassager, O. Dynamics of Polymeric Liquids. Vol. 1: Fluid Mechanics; John Wiley and Sons Inc.: New York, NY, USA, 1987.

29. Alves, L.D.B.; Hirata, S.C.; Schuabb, M.; Barletta, A. Identifying linear absolute instabilities from differential eigenvalue problems using sensitivity analysis. J. Fluid Mech.; 2019; 870, pp. 941-969. [DOI: https://dx.doi.org/10.1017/jfm.2019.275]

30. Li, Z.; Khayat, R. Finite-amplitude Rayleigh-Bénard convection and pattern selection for viscoelastic fluids. J. Fluid Mech.; 2005; 529, pp. 221-251. [DOI: https://dx.doi.org/10.1017/S0022112005003563]

31. Sureshkumar, R.; Beris, A. Linear stability analysis of viscoelastic Poiseuille flow using an Arnoldi-based orthogonalization algorithm. J. Non-Newton. Fluid Mech.; 1994; 56, pp. 151-182. [DOI: https://dx.doi.org/10.1016/0377-0257(94)01279-Q]