Onset of Primary and Secondary Instabilities of Viscoelastic Fluids Saturating a Porous Layer Heated from below by a Constant Flux

Abstract

We analyze the thermal convection thresholds and linear characteristics of the primary and secondary instabilities for viscoelastic fluids saturating a porous horizontal layer heated from below by a constant flux. The Galerkin method is used to solve the eigenvalue problem by taking into account the elasticity of the fluid, the ratio between the viscosity of the solvent and the total viscosity of the fluid and the lateral confinement of the medium. For the primary instability, we found out that depending on the rheological parameters, two types of convective structures may appear when the basic conductive solution loses its stability: stationary long wavelength instability as for Newtonian fluids and oscillatory convection. The effect of the lateral confinement of the porous medium by adiabatic walls is to stabilize the oblique and longitudinal rolls and therefore selects transverse rolls at the onset of convection. In the range of the rheological parameters where stationary long wave instability develops first, we use a parallel flow approximation to determine analytically the velocity and temperature fields associated with the monocellular convective flow. The linear stability analysis of the monocellular flow is performed, and the critical conditions above which the flow becomes unstable are determined. The combined influence of the viscoelastic parameters and the lateral confinement on the characteristics of the secondary instability is quantified. The major new findings concerning the secondary instabilities may be summarized as follows: (i) For concentrated viscoelastic fluids, computations showed that the most amplified mode of convection corresponds to oscillatory transverse rolls, which appears via a Hopf bifurcation. This pattern selection is independent of both the fluid elasticity and the lateral confinement of the porous medium; (ii) For diluted viscoelastic fluids, the preferred mode of convection is found to be oscillatory transverse rolls for a very laterally-confined medium. Otherwise, stationary or oscillatory longitudinal rolls may develop depending on the fluid elasticity. Results also showed the destabilizing effect of the relaxation fluid elasticity and the stabilizing effect of the viscosity ratio for the onset of both primary and secondary instabilities.

1. Introduction

The study of viscoelastic fluids has applications in a number of processes that occur in industry, such as the extrusion of polymer fluids, solidification of liquid crystals, suspension solutions and petroleum activities. In contrast to the case of Newtonian fluids, the study of the thermal convection of viscoelastic fluids in porous media is limited. In rheology, one crucial problem is the formulation of the constitutive equations regarding viscoelastic fluid flows in porous media. Recently, a modified Darcy’s law was employed to study the stability of a viscoelastic fluid in a horizontal porous layer using linear and nonlinear stability theory [1,2,3,4,5,6,7,8,9]. Kim et al. [1] and Yoon et al. [2] performed a linear stability analysis and showed that in viscoelastic fluids, such as polymeric liquids, a Hopf bifurcation, as well as a stationary bifurcation may occur depending on the magnitude of the viscoelastic parameter. From the nonlinear point of view, Kim et al. [1] carried out a nonlinear stability analysis by assuming a densely-packed porous layer and found that both stationary and Hopf bifurcations are supercritical relative to the critical heating rate. The question of whether standing or traveling waves are preferred at onset has been fully addressed by Hirata et al. [4]. The three-dimensional convective and absolute instabilities of a viscoelastic fluid in the presence of a horizontal pressure gradient have been analyzed by Hirata and Ouarzazi [5]. Alves et al. [6] studied the effect of viscous dissipation of viscoelastic fluids at the onset of convection. In addition to its theoretical interest, Delenda et al. [7] have shown that viscoelastic convection in porous media may be useful for industrial applications for the separation of species of viscoelastic solutions. The introduction of a porous packing allows one to control the average vertical convective velocity and to generate a homogeneous convection current, improving the separation of species. Fu et al. [8] performed direct numerical simulations on two-dimensional thermal convection of a viscoelastic fluid saturating a porous square cavity. Their numerical experiments revealed the existence of a second transition from oscillatory convection to a stationary one followed by a third transition to oscillatory convection for some combinations of rheological parameters, while these successive transitions never occur for other combinations of viscoelastic parameters. Taleb et al. [9] used both theoretical and numerical approaches and obtained a global picture and bifurcations diagrams on possible successive bifurcations of convection patterns in a square porous cavity saturated by a viscoelastic fluid.

All of the above investigations considered conventional boundary conditions, namely impermeable isothermal horizontal plates, commonly known as Horton–Rogers–Lapwood convection. With a configuration of an isothermal upper wall and impermeable adiabatic side walls, Niu et al. [10] were the first to introduce at the bottom plate a general variable heat flux boundary condition expressed in terms of the Biot number ($Bi$). They proved by a linear stability analysis that when $Bi\to 0$, the thermal convection properties of a viscoelastic fluid of an Oldroyd fluid type approach that subject to uniform heat flux heating. On the other hand, when $Bi\to \infty$, linear properties of convection coincide with those corresponding to isothermal heating. Therefore, convection with different thermal boundary conditions may be analyzed by adjusting the Biot number. However, to the best of our knowledge, no results have been published for thermal convection of viscoelastic fluids in situations more commonly met in engineering applications, namely when the porous medium is subjected, both in the upper and lower walls, to a general variable heat flux boundary condition. Therefore, the objective of this work is to fill a part of this gap by investigating the onset of three-dimensional primary and secondary instabilities of a viscoelastic fluids under the assumption that the upper and lower horizontal walls are kept at a constant flux as a limiting case of a more general variable heat flux boundary condition.

For Newtonian fluids, the stability of an infinite porous layer with different boundary conditions was studied by Nield [11] and is well documented in Section $6.2$ of the book by Nield and Bejan [12]. For the case of a porous medium heated from the bottom and cooled from the top by a constant heat flux, Nield [11] found that the critical Rayleigh number at the onset of convection is approximately 12 with a vanishing critical wavenumber. Mamou et al. [13] extended the work of Nield [11] by taking into account the effect of the anisotropy of the porous medium. Mojtabi and Rees [14] studied the case where the impermeable boundary walls have a finite thickness. They analyzed the combined influence on the onset of convection of the ratio between the thermal conductivity of the horizontal walls and the thermal conductivity of the porous medium, as well as the ratio between the thickness of the horizontal walls and the thickness of the porous layer. The problem of thermal convection in a rotating horizontal layer of a porous medium was analyzed by Falsaperla et al. [15] for a single fluid and by Falsaperla et al. [16] in the case of a binary fluid mixture. In both papers, the authors examined the effects on the linear properties of the convection of boundary conditions for the temperature of the Newton–Robin type with heat flux prescribed as a limiting case.

Kimura et al. [17] investigated secondary instabilities for a Newtonian fluid saturating a porous medium heated from below by a constant flux. For a Rayleigh number larger than its critical value 12 above which the conduction state looses its stability against long wave instability, these authors used the parallel flow approximation and obtained a nonlinear solution, which corresponds to a monocellular flow. Two-dimensional numerical results were also presented to test the validity of the approximated nonlinear solution. In addition, they analyzed its stability against three-dimensional disturbances and showed that the monocellular flow is linearly stable to transverse disturbances for Rayleigh number as high as 506, at which point, a Hopf bifurcation sets in. However, further analysis indicated that an exchange of stability due to longitudinal disturbances will occur much sooner at a Rayleigh number equal to $311.53$.

This contribution aims to understand how the viscoelastic character of the fluid influences the properties of convection at the onset of primary and secondary instabilities, when the porous layer is heated from below by a constant flux. Therefore, this work may be viewed as an extension to viscoelastic fluids of the work done by Kimura et al. [17].

The paper is organized as follows. After presenting the governing equations in Section 2, the stability of the conductive state is studied in Section 3 by considering steady, as well as oscillatory three-dimensional perturbations. Section 4 is devoted to the discussion of the combined effects of the viscoelastic parameters and the lateral aspect ratio of the porous medium on the pattern selection at the onset of secondary instabilities. Finally, in Section 5, the main conclusions of the present study are presented.

2. Mathematical Formulation

Let us consider an isotropic and homogeneous porous cavity of thickness e, height H and width W (see Figure 1). The porous medium is saturated by an Oldroyd-B fluid, and we assume that the solid matrix is in local thermal equilibrium with the fluid. The upper and lower horizontal walls are kept at constant flux, while the lateral vertical walls are considered adiabatic. The solid walls of the domain $\mathsf{\Omega }=[0,W]\times [0,e]\times [0,H]$ are considered impermeable. We assume that the Oberbeck–Boussinesq approximation holds.

There are several ways to obtain macroscopic laws for polymeric flows in a porous medium: by direct numerical simulations of viscoelastic flows in a specific pore geometry model (a good review of these studies can be found in [18]) or by analytical means. In general, the former is the most commonly-used method for the derivation of macroscopic laws. It can be divided in two techniques: the REV method (representative elementary volume method) and the homogenization theory. The starting point for the two techniques is a local description in a pore scale. The pore space is assumed to be saturated by an incompressible viscoelastic fluid. For slow flows, the momentum balance equation can be linearized:

$$\rho \frac{\partial {\mathbf{U}}^{*}}{\partial t^{*}}=-\nabla p^{*}+\rho \mathbf{g}+\nabla ·\tilde{\tau }$$

(1)

where ${\mathbf{U}}^{*}$ is the fluid velocity field, $p^{*}$ is the pressure, $\tilde{\tau }$ is the stress tensor and $\mathbf{g}$ is the gravity field.

In Newtonian incompressible fluids, the constitutive relation between stress tensor $\tilde{\tau }$ and strain tensor $\tilde{D}(D_{i,j}=[u_{i,j}^{*}+u_{j,i}^{*}]/2)$ is the Newtonian law $\tilde{\tau }=2{\mu }_N\tilde{D}$, where ${\mu }_N$ is the dynamic viscosity, and, in this case, the relation $\nabla ·\tilde{\tau }={\mu }_N{\nabla }^2{\mathbf{U}}^{*}$ is obtained.

The simplest rheological model useful for describing dilute viscoelastic solutions at high rates of deformation is the Oldroyd-B model. In this model, $\tilde{\tau }$ and $\tilde{D}$ for viscoelastic fluids composed of a Newtonian solvent and a polymeric solute of “Newtonian” viscosity ${\mu }_s$ and “elastic viscosity” ${\mu }_p$ [19], respectively, are given by:

$$\tilde{\tau }={\tilde{\tau }}_s+{\tilde{\tau }}_p$$

(2)

with:

$${\tilde{\tau }}_s=2{\mu }_s\tilde{D}$$

(3)

and:

$$(1+{\lambda }_1^{*}\frac{\partial }{\partial t^{*}}){\tilde{\tau }}_p=2{\mu }_p\tilde{D}$$

(4)

where ${\lambda }_1^{*}$ represents the relaxation time. Then, by combining (2)–(4), we obtain the constitutive equation:

$$(1+{\lambda }_1^{*}\frac{\partial }{\partial t^{*}}){\tilde{\tau }}_p=2\mu (1+{\lambda }_2^{*}\frac{\partial }{\partial t^{*}})\tilde{D}$$

(5)

where the dynamic viscosity $\mu$ and the retardation time ${\lambda }_2^{*}$ are related to ${\mu }_s$ and ${\mu }_p$ by:

$$\mu ={\mu }_s+{\mu }_{p}and{\lambda }_2^{*}/{\lambda }_1^{*}={\mu }_s/({\mu }_s+{\mu }_p).$$

An Oldroyd-B fluid may thus be characterized by three parameters: the dynamic viscosity $\mu$, the relaxation ${\lambda }_1^{*}$ and the retardation ${\lambda }_2^{*}$ times. The relation $\mathsf{\Gamma }={\lambda }_2^{*}/{\lambda }_1^{*}$ may also be used instead of ${\lambda }_2^{*}$.

In order to derive a macroscopic filtration law based on the Oldroyd constitutive equation, we have to introduce, in a fashion similar to that done in the derivation of the classical Darcy law for Newtonian fluids, the superficial velocity ${\mathbf{V}}^{*}$ defined by:

$$<{\mathbf{V}}^{*}>=\frac{1}{\mathit{v}}\underset{\mathit{v}}{\int }{\mathbf{U}}^{*}d\mathit{v}$$

(6)

where $\mathit{v}$ is a representative elementary volume (REV) of the medium, incorporating both solid and fluid. Using the REV method by averaging Equation (1) and taking into account Equation (6) leads to:

$$\frac{\mu }{K}(1+{\lambda }_2^{*}\frac{\partial }{\partial t^{*}}){\mathbf{V}}^{*}+(1+{\lambda }_1^{*}\frac{\partial }{\partial t^{*}})(\nabla P^{*}-\rho \mathbf{g})=0,$$

(7)

where K is the permeability.

Under the assumption of low Reynolds number based on the pore dimension, the generalized Darcy’s law (7) is also derived by [20] using a homogenization theory.

The fluid density $\rho$ obeys the state law:

$$\rho ={\rho }_0(1-{\beta }_T(T^{*}-T_0^{*}))$$

(8)

where ${\rho }_0$ is the fluid density at temperature $T_0^{*}$, which is chosen here as the temperature at the geometric center of the cavity, and ${\beta }_T$ is the thermal expansion coefficient.

Assuming thermal equilibrium between the solid and the fluid phases and neglecting the viscous dissipation effects, energy and continuity equations can then be written as [1,2,3,4,5,6,7,8,9,10]:

$$\frac{{\left(\rho c\right)}_{sf}}{{\left(\rho c\right)}_f}\frac{\partial T^{*}}{\partial t^{*}}+{\mathbf{V}}^{*}·\nabla T^{*}=\nabla ·(\alpha \nabla T^{*})$$

(9)

$$\nabla ·{\mathbf{V}}^{*}=0$$

(10)

The boundary conditions at the impermeable horizontal walls kept at a constant flux q and the impermeable insulated vertical walls are:

$$\begin{array}{c}-k_T\frac{\partial T^{*}}{\partial z}=q\mathrm{at}z=0,H,\hfill \\ \frac{\partial T^{*}}{\partial x}=0\mathrm{at}x=0,W,\hfill \\ \frac{\partial T^{*}}{\partial y}=0\mathrm{at}y=0,e,\hfill \\ \mathbf{V}·\mathbf{n}=0\mathrm{at}\partial \mathsf{\Omega }.\hfill \end{array}$$

(11)

Here, ($\rho c)$, $\mu$, $\nu$, $k_T$, $\alpha =k_T/{\left(\rho c\right)}_f$ are respectively the heat capacity per unit volume, the dynamic and kinematic viscosity of the fluid, the effective thermal conductivity and the effective thermal diffusivity. Subscript sf refers to an effective quantity, while f refers to the fluid alone.

We choose H, $k_T/(H{\left(\rho c\right)}_f)$, $H^2{\left(\rho c\right)}_{sf}/k_T$, $k_T\mu /(K{\left(\rho c\right)}_f)$ and $qH/k_T$ as reference quantities for the length, velocity, time, pressure and temperature difference ($T^{*}-T_0^{*}$). With this scaling, the following set of dimensionless equations is obtained:

$$\nabla ·\mathbf{V}=0$$

(12)

$$(1+\mathsf{\Gamma }{\lambda }_1\frac{\partial }{\partial t})\mathbf{V}+(1+{\lambda }_1\frac{\partial }{\partial t})(\nabla P-RaT{\mathbf{e}}_z)=0,$$

(13)

$$\frac{\partial T}{\partial t}+\mathbf{V}·\nabla T={\nabla }^{2}T$$

(14)

The dimensionless boundary conditions are:

$$\begin{array}{c}\frac{\partial T}{\partial z}=-1\mathrm{at}z=0,1,\hfill \\ \frac{\partial T}{\partial x}=0\mathrm{at}x=\pm \frac{A}{2},\hfill \\ \frac{\partial T}{\partial y}=0\mathrm{at}y=0,a,\hfill \\ \mathbf{V}·\mathbf{n}=0\mathrm{at}\partial \mathsf{\Omega }.\hfill \end{array}$$

(15)

The dimensionless parameters are: the filtration Rayleigh number:

$$Ra=\frac{{\beta }_{T}gK{H}^{2}q}{\alpha \nu k_T}$$

(16)

the horizontal and lateral aspect ratios:

$$A=W/H,a=e/H$$

(17)

the relaxation time:

$${\lambda }_1={\lambda }_1^{*}{k}_T/(H^2{\left(\rho c\right)}_{sf})$$

(18)

and the ratio $\mathsf{\Gamma }$ that varies in the interval $[0,1]$:

$$\mathsf{\Gamma }={\lambda }_2^{*}/{\lambda }_1^{*}$$

(19)

This model reduces to the Maxwell model in the limit $\mathsf{\Gamma }\to 0$ and to the Newtonian model in the limit $\mathsf{\Gamma }\to 1$.

In the following, we will examine the stability of the conductive state (the primary instability), as well as the stability of the monocellular flow (the secondary instability).

3. Primary Stationary and Oscillatory Instabilities

In the conductive regime, the basic solution is a motionless state $\mathbf{V}=\mathbf{0}$ with a vertical thermal stratification $T_0=-z+\frac{1}{2}$.

The aim of this section is to perform a temporal stability analysis of the conductive state with respect to both stationary and oscillatory disturbances.

3.1. Infinite Aspect Ratios

To investigate the stability of the basic solution, infinitesimal three-dimensional perturbations are super-imposed onto the basic solution:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\hfill \mathbf{V}={\mathbf{V}}_0+\mathbf{v}(x,y,z,t)\\ \hfill T=T_0+\theta (x,y,z,t)\\ \hfill P=P_0+p(x,y,z,t)\end{array}\right.\end{array}$$

(20)

We first assume very large aspect ratios $A(A\to \infty )$ and $a(a\to \infty )$. The three-dimensional disturbance quantities are expressed as:

$$(u,v,w,\theta ,p)=\left[\tilde{u}\left(z\right),\tilde{v}\left(z\right),\tilde{w}\left(z\right),\tilde{\theta }\left(z\right),\tilde{p}\left(z\right)\right]exp(ikx+ily-i\omega t)$$

(21)

where k and l are the wave numbers in the x and y directions, respectively, and the temporal growth rate of unstable perturbations is given by the imaginary part of the complex frequency $\omega ={\omega }_r+i{\omega }_i$. Therefore, the neutral temporal stability curve is obtained for ${\omega }_i=0$, which selects dominant modes at the onset of convection.

Substituting Equations (20), (21) into (12)–(15), linearizing the equations and applying the curl twice to the momentum balance equation, one can obtain:

$$(1-i\omega \mathsf{\Gamma }{\lambda }_1)(D^2-{\tilde{k}}^2)\tilde{w}+Ra(1-i\omega {\lambda }_1){\tilde{k}}^2\tilde{\theta }=0$$

(22)

$$-i\omega \tilde{\theta }-\tilde{w}-(D^2-{\tilde{k}}^2)\tilde{\theta }=0$$

(23)

where $D=\frac{d}{dz}$ and ${\tilde{k}}^2=k^2+l^2$. The corresponding boundary conditions take the form:

$$\tilde{w}=0,\frac{d\tilde{\theta }}{dz}=0\mathrm{at}z=0,1.$$

(24)

The system (40), (41) is solved by means of the Galerkin method using the following expansions:

$$\tilde{w}\left(z\right)=\sum _{n=1}^M{w}_n\sin \left(n\pi z\right)$$

(25)

$$\tilde{\theta }\left(z\right)=\sum _{n=1}^M{\theta }_n\cos \left[(n-1)\pi z\right]$$

(26)

The number M of modes is chosen so that the quantitative convergence is secured.

We first examine the single-term Galerkin expansion, or single-term truncation, for the present problem. It allows the generation of analytical expressions that can guide future laboratory experiments, owing to its simplicity. Even when their error in quantitative predictions was high, their qualitative behavior prediction was accurate enough to guide additional numerical efforts. Setting $M=1$ in (25), (26), an algebraic dispersion relation is obtained:

$$\begin{array}{c}\hfill (1-i\omega \mathsf{\Gamma }{\lambda }_1)({\pi }^2+{\tilde{k}}^2)({\tilde{k}}^2-i\omega )-\frac{8}{{\pi }^2}Ra({\tilde{k}}^2(1-i\omega {\lambda }_1)=0,\end{array}$$

(27)

As the viscoelastic parameters appear only in front of a time derivative in the momentum Equation (13) and in front of $-i\omega$ in the approximated dispersion Equation (27), the elasticity of the fluid cannot influence the properties of a stationary instability. Consequently, the characteristics of the stationary instability are the same as for Newtonian fluids. For such fluids, linear instability analysis has been considered by Nield [11] and has provided quantitative information on the stability condition when the porous layer is supposed infinite in the x and y directions.

On evaluating the minimum of the Rayleigh number defined by (27), we obtain the critical values for the onset of stationarity.

$$\begin{array}{c}\hfill R{a}_c^s=\frac{{\pi }^4}{8}\simeq 12.176and{\tilde{k}}_c^s=0\end{array}$$

(28)

The critical value $R{a}_c^s=\frac{{\pi }^4}{8}$ obtained with a single-term truncation is only $1.5$ percent higher than the exact critical $R{a}_c^s=12$ found in [11].

By setting $\omega ={\omega }_r\ne 0$ in (27) and minimizing $Ra$ with respect to $\tilde{k}$, we obtain a Hopf bifurcation giving rise to an oscillatory convection if the Rayleigh number exceeds the critical value:

$$\begin{array}{c}\hfill R{a}_c^{osc}=\frac{{\pi }^2}{8{\lambda }_1}{(1+\pi \sqrt{{\lambda }_1\mathsf{\Gamma }})}^2.\end{array}$$

(29)

The corresponding critical wave number and critical frequency for oscillatory mode are:

$$\begin{array}{c}\hfill {\tilde{k}}_c^2=\frac{\pi }{\sqrt{{\lambda }_1\mathsf{\Gamma }}}\end{array}$$

(30)

and

$$\begin{array}{c}\hfill {\omega }_c^2=\frac{1}{{\lambda }_1^2\mathsf{\Gamma }}[{\tilde{k}}_c^2({\lambda }_1-{\lambda }_1\mathsf{\Gamma })-1]\end{array}$$

(31)

As ${\omega }_c$ is real, the system exhibits oscillatory convection only if ${\omega }_c^2>0$. This condition is satisfied if:

$$\begin{array}{c}\hfill {\lambda }_1>\frac{\mathsf{\Gamma }}{{\pi }^2{(1-\mathsf{\Gamma })}^2}={\lambda }_1^{*}\end{array}$$

(32)

On the other hand, the bifurcation to oscillatory convection precedes the bifurcation to steady convection, i.e., the critical Rayleigh number for oscillatory instability is less than its value for stationary instability if:

$$\begin{array}{c}\hfill {\lambda }_1>\frac{1}{{\pi }^2{(1-\sqrt{\mathsf{\Gamma }})}^2}={\lambda }_1^{**}\end{array}$$

(33)

As computations showed that ${\lambda }_1^{*}<{\lambda }_1^{**}$, we concluded that the system selects oscillatory convection as a first instability if ${\lambda }_1>{\lambda }_1^{**}$. The border separating stationary and oscillatory instabilities is given by the curve ${\lambda }_1={\lambda }_1^{**}$, as is shown in the plane ($\mathsf{\Gamma }$, ${\lambda }_1$) of Figure 2.

The single-term Galerkin truncation does not present accurate results over the entire range of parameters under consideration, and higher-order truncation becomes necessary. Having used the Galerkin expansion (25) and (26) with $M=5$, we obtain results with a very good agreement with those obtained in [11] in the case of stationary instability. Figure 3a represents the marginal stability curve in the ($\tilde{k}$, Ra) plane and shows that a long wave instability (i.e., the critical wave number ${\tilde{k}}_c=0$) may develop if the Darcy–Rayleigh number exceeds the critical value $R{a}^s=12.009$ in accordance with the critical value $R{a}^s=12$ found in [11].

In Figure 3b, we plot the curve of neutral stability for oscillatory mode of convection in the $(\tilde{k},Ra)$ plane for $\mathsf{\Gamma }=0.02$ and different values of the elasticity number ${\lambda }_1=0.4;0.5;0.6$. It can be seen from this figure that the minimum value of the Rayleigh number is lower than the critical Rayleigh number $Ra=12$ needed to trigger steady long wave instability. Therefore, oscillatory instability may set up as a first convective pattern instead of steady long wave instability.

The dependence of the critical Rayleigh number and the critical frequency at the onset of oscillatory convection on the elasticity number ${\lambda }_1$ for fixed values of $\mathsf{\Gamma }$ is numerically determined, and the results are plotted in Figure 4a,b respectively.

It is clear from Figure 4a that an increase in ${\lambda }_1$ leads to flow destabilization, i.e., to a reduction in the respective critical Rayleigh number. Figure 4a also shows the stabilizing effect of the ratio $\mathsf{\Gamma }$. Moreover, as is seen in Figure 4a, for a fixed value of $\mathsf{\Gamma }$, there exists a particular value of ${\lambda }_1={\lambda }_1^f$ where the critical Rayleigh numbers for the onsets of both oscillatory and stationary convection coincide, and therefore, a codimension two bifurcation occurs. For ${\lambda }_1>{\lambda }_1^f$, Figure 4b shows that the critical frequency decreases with the decrease of the fluid elasticity or the increase of the viscosity ratio.

3.2. Effect of Lateral Confinement on Pattern Selection

This section is devoted to investigating the effect of the lateral confinement of the porous cavity by assuming a very large aspect ratio $A(A\to \infty )$ and finite lateral aspect ratio a. The three-dimensional disturbance quantities respecting the boundary conditions (15) are expressed as:

$$(u,w,\theta ,p)=\left[\tilde{u}\left(z\right),\tilde{w}\left(z\right),\tilde{\theta }\left(z\right),\tilde{p}\left(z\right)\right]exp(ikx-i\omega t)cos(L\pi y/a)$$

(34)

$$v=\tilde{v}\left(z\right)exp(ikx-i\omega t)sin(L\pi y/a)$$

(35)

The governing equations are still the system (40) and (41), except with l now replaced by $L\pi /a$ where the integer L is the number of rolls in the y direction.

We begin the study by considering the stability of the conductive state against stationary rolls with axes parallel to the x direction, called longitudinal rolls (LRs). Steady longitudinal rolls are characterized by $k=0$, $L\ne 0$ and ${\omega }_r=0$. The dependence of the critical Rayleigh number at the onset of LRs on the lateral aspect ratio a for different numbers L of rolls is displayed in Figure 5a. For comparison, we also represent in the same figure the threshold of the steady long wave instability. The threshold of steady three-dimensional patterns in the form of oblique rolls (i.e., $k\ne 0$, $L\ne 0$ and ${\omega }_r=0$) is bounded by the thresholds of the two limiting cases: the steady long wave instability and steady LRs.

We remark that the mode $L=1$ is the most unstable mode for LRs. As is expected, we note that the critical Rayleigh number increases as a decreases, meaning that the lateral confinement stabilizes the conductive state against longitudinal rolls. We also note that as $a\to \infty$, the limiting value of $Ra=12$ is reached monotonically, and an infinity of modes may be simultaneously unstable in this limit. Consequently, a relatively moderate lateral confinement is necessary to select the long wave instability, which corresponds in real experiments to a monocellular flow in the x direction.

Now, we consider the effect of the lateral confinement on the stability of the conductive state against oscillatory LRs defined by $k=0$, $L\ne 0$ and ${\omega }_r\ne 0$. Numerical results for neutral stability curves of oscillatory LRs with $L=1$, $L=2$, $L=3$ and $L=4$ are shown in Figure 5b as functions of the lateral aspect ratio a. These curves have a parabolic shape and intersect at some particular values of a, indicating that the true critical Rayleigh number strongly depends on both a and L for fixed rheological parameters. The destabilizing oscillatory LR mode changes in the intersection points of neural curves from a mono-cellular flow to a two-cellular flow, and so on, as the lateral aspect ratio a increases. In addition, the behavior of the critical Rayleigh number is non-monotonic as a increases. We also note that the maximum of critical Rayleigh number decreases as a increases and tends asymptotically to the critical threshold found in the unbounded case ($a\to \infty$). The results are therefore in contrast to the case of stationary LRs where the dominant mode corresponds to $L=1$ independently of the lateral confinement.

In Figure 5b, the critical Rayleigh number at the onset of oscillatory transverse rolls (TRs) characterized by $k\ne 0$, $L=0$ and ${\omega }_r\ne 0$ is indicated by the horizontal line. As can be seen from this figure, finite values of a stabilize oscillatory LRs and may select oscillatory TRs as a dominant mode of convection.

4. Secondary Instabilities

4.1. Nonlinear Solution and Formulation of Its Linear Stability

According to the above linear stability analysis, we found that a stationary bifurcation occurs giving rise to a convective pattern in the form of a long wave instability in the x direction provided that the elasticity number ${\lambda }_1$ does not exceed a particular value ${\lambda }_1^f$, which depends on the viscosity ratio $\mathsf{\Gamma }$. In that case, the viscoelastic fluid behaves like a Newtonian fluid. Consequently, the nonlinear solution in the regime of steady long wave convection is the same regardless of weather or not the fluid is viscoelastic.

By assuming a shallow cavity $A\gg 1$ and by using the parallel flow approximation [21,22,23], Kimura et al. [17] found that the analytical solution for the monocellular flow consists of: a horizontal asymmetric velocity with a zero mean along any vertical section,

$$U\left(z\right)=\frac{1}{2}RaC(1-2z)$$

(36)

and a vertical, as well as a horizontal stratification of the temperature,

$$T_0(x,y,z)=Cx+\Theta \left(z\right)$$

(37)

with:

$$\Theta \left(z\right)=\frac{1}{2}Ra{C}^2(\frac{z^2}{2}-\frac{z^3}{3}-\frac{1}{12})-z+\frac{1}{2}$$

(38)

and:

$$C=\pm \sqrt{\frac{10}{Ra}(1-\frac{12}{Ra})}$$

(39)

where C is negative or positive according to whether the flow is clockwise or counter-clockwise, and both solutions are possible depending on the initial conditions.

From Equation (39), it is seen that no motion may be induced inside the cavity for $Ra<12$. For the case of a porous medium heated from the bottom and cooled from the top by a constant heat flux, a critical Rayleigh number of $Ra=12$ for the onset of convection was predicted by Nield [11]. This result is in agreement with the prediction of Equation (39).

For finite aspect ratio, Kimura et al. [17] performed two-dimensional numerical simulations of the full problem. Their numerical results show that the conductive state is stable when the Rayleigh number is smaller than 12. Computations carried out for Ra in excess of 12 were found to agree with analytical solutions (36)–(38).

The equations governing the linear stability of the monocellular flow are obtained by the same previous approach used for the stability of the conductive basic solution. By assuming very large aspect ratios $A(A\to \infty )$ and $a(a\to \infty )$, the following system is obtained:

$$(1-i\omega \mathsf{\Gamma }{\lambda }_1)(D^2-{\tilde{k}}^2)\tilde{w}+Ra(1-i\omega {\lambda }_1){\tilde{k}}^2\tilde{\theta }=0$$

(40)

$$-i\omega \tilde{\theta }+\tilde{w}D{T}_0+ik\tilde{\theta }{U}_0-(D^2-{\tilde{k}}^2)\tilde{\theta }=0$$

(41)

where we substitute $U_0$ and $T_0$ by their explicit expressions (36)–(38).

The corresponding boundary conditions take the form:

$$\tilde{w}=0,\frac{d\tilde{\theta }}{dz}=0\mathrm{at}z=0,1.$$

(42)

On the other hand, if we assume a very large aspect ratio A and a finite value of the lateral aspect ratio a, the governing equations are still the system (40) and (41) where ${\tilde{k}}^2$ is replaced by $k^2+L^2{\pi }^2/a^2$.

The resulting linear stability problem is solved by means of the Galerkin method, using the expansion (25) and (26) at the order $M=30$.

4.2. Results for Newtonian Fluids

To verify the accuracy of our numerical results based on the Galerkin expansion to the order $M=30$, we perform a test for the limiting case of a Newtonian fluid and compare the results with those obtained by Kimura et al. [17]. In the first instance, two-dimensional disturbances, corresponding to $l=0$, were considered. We found out that for the Newtonian fluid, the base velocity and temperature profiles (36)–(38) are stable for values of $Ra$ less than $R{a}_{c2}^T=506.27$ as shown by the neutral stability curve represented in Figure 6a. At this critical Rayleigh number occurs an instability via a Hopf bifurcation to oscillatory TRs with a critical frequency ${\omega }_{c2}^T=138.24$ and a critical wave number $k_{c2}^T=4.80$. These results are in a good agreement with those obtained in [17] by using a shooting method, namely $R{a}_{c2}^T=506.07$, ${\omega }_{c2}^T=138.92$ and $k_{c2}^T=4.82$.

On the other hand, Kimura et al. [17] considered three-dimensional disturbances with the value of the y-wave number l being gradually increased from zero. For $l>0$, the stability analysis indicates that the monocellular flow will be destabilized not by a Hopf bifurcation, but by an exchange of stability for which the x-wave number k vanishes. In that case the threshold of the appearance of steady longitudinal rolls as a secondary instability is found to be $R_{c2,s}^L\approx 311.53$. Since this critical Rayleigh number is much lower than any of those for the Hopf bifurcations obtained when $k\ne 0$, Kimura et al. [17] concluded that the monocellular flow will in fact be destabilized by longitudinal, rather than transverse, disturbances.

In the second instance, three-dimensional disturbances, corresponding to $k\ne 0$ and $l\ne 0$, were considered in this study. Numerical results performed by assuming infinite aspects ratios A and a indicated that the most unstable mode corresponds to $k=0$ and $l\ne 0$. The corresponding critical Rayleigh number above which this most unstable mode is in the form of steady LRs is $R_{c2,s}^L=313.107$. Once again, this critical value agrees very well with $R_{c2,s}^L\approx 311.53$ obtained in [17].

In the third instance, the effect of the confinement of the porous medium in the y direction is explored. We plot in Figure 6b the critical Rayleigh number against the aspect ratio for several of the leading modes, from which it is clear that $(L=1)$ remains the destabilizing mode, ahead of the other modes $(L>1)$, and that the order of these modes, in the sense that $R{a}_c\left(L\right)<R{a}_c(L+1)$, is preserved as a increases. In particular, we also note that as $a\to \infty$, the limiting value of $R_{c2,s}^L=313.107$ is approached monotonically. Figure 6b also shows that the curve corresponding to steady longitudinal mode with $L=1$ intersects the line $R{a}_{c2}^T=506.07$ at a particular value of the lateral aspect ratio $a=a^{*}$. This means that perturbations promote the appearance of oscillatory TRs provided that $a<a^{*}$ or stationary LRs otherwise.

4.3. Results for Viscoelastic Fluids

4.3.1. Hopf Bifurcation to Transverse Rolls

In order to study the influence of viscoelastic parameters on the secondary instability, we first computed the bifurcation line from a stationary monocellular convective pattern to oscillatory TRs $(l=0)$ for either a fixed value of the elasticity number ${\lambda }_1$ with varying values of the viscosity ratio $\mathsf{\Gamma }$ or a fixed value of $\mathsf{\Gamma }$ with varying values of ${\lambda }_1$. With regard to the question of the influence of the viscosity ratio $\mathsf{\Gamma }$ for a viscoelastic fluid with a relaxation time ${\lambda }_1=0.1$ on the onset of a secondary instability in the form of oscillatory TRs, Figure 7a illustrates the behavior of neutral stability curves in the $(k,Ra)$ plane for $\mathsf{\Gamma }=0.75$, $\mathsf{\Gamma }=0.5$ and $\mathsf{\Gamma }=0.3$. For a comparison, the Newtonian case $(\mathsf{\Gamma }=1)$ is also represented on Figure 7a.

We note in this figure that the minimum of neutral stability curves increases when $\mathsf{\Gamma }$ is augmented to reach the critical value for Newtonian fluids in the limit of $\mathsf{\Gamma }=1$. Physically, this result means that concentrated polymeric solutions with a small viscosity ratio $\mathsf{\Gamma }$ favor the appearance of oscillatory multicellular flow convection as a secondary instability. On the other hand, for diluted viscoelastic solutions, more heating is needed to trigger the secondary instability.

We report in

Table 1. Critical Rayleigh number $R{a}_{c2}^T$, frequency ${\omega }_{c2}^T$ and wave number $k_{c2}^T$ at the onset of moving transverse rolls as a secondary instability for ${\lambda }_1=0.1$ and different values of $\mathsf{\Gamma }$.

$\mathsf{\Gamma }$	${\mathit{Ra}}_{\mathit{c}\mathbf{2}}^{\mathit{T}}$	${\mathit{\omega }}_{\mathit{c}\mathbf{2}}^{\mathit{T}}$	${\mathit{k}}_{\mathit{c}\mathbf{2}}^{\mathit{T}}$
Newtonian	506.27	138.24	4.8
0.75	358.62	115.209	4.660
0.70	329.48	110.448	4.630
0.65	300.89	105.918	4.610
0.60	272.90	101.395	4.590
0.55	245.58	96.897	4.570
0.50	219.04	92.825	4.570
0.45	193.39	88.902	4.575
0.40	168.81	85.603	4.610
0.35	145.47	83.112	4.685
0.30	123.55	81.578	4.805

the computed results of critical Rayleigh number $R{a}_{c2}^T$, critical frequency ${\omega }_{c2}^T$ and critical wave number $k_{c2}^T$ at the onset of the secondary instability organized as oscillatory TRs for ${\lambda }_1=0.1$ and different values of $\mathsf{\Gamma }$. Table 1 shows a strong stabilizing effect of the viscosity ratio. The values of the critical oscillatory frequency decrease with decreasing $\mathsf{\Gamma }$. This implies that emerging transversal convection rolls have a larger time-period and move with a larger phase velocity when the polymer concentration is high.

We now present results corresponding to the influence of the fluid elasticity by considering the properties of the emerging oscillatory TRs at different values of ${\lambda }_1$ for a fixed value of $\mathsf{\Gamma }$. Figure 7b presents neutral stability curves for $\mathsf{\Gamma }=0.75$, a typical viscosity ratio value for Boger fluids and different values of the relaxation time ${\lambda }_1=0.1$, ${\lambda }_1=0.35$ and ${\lambda }_1=0.5$. We note from this figure that the neural stability curves are nearly superposed when ${\lambda }_1$ is increased, meaning that beyond ${\lambda }_1=0.1$, the increase in the fluid elasticity has little influence on the critical Rayleigh number at the onset of oscillatory TRs.

Table 2. Critical Rayleigh number $R{a}_{c2}^T$, frequency ${\omega }_{c2}^T$ and wave number $k_{c2}^T$ at the onset of moving transverse rolls as a secondary instability for $\mathsf{\Gamma }=0.75$ and different values of ${\lambda }_1$.

${\mathit{\lambda }}_{\mathbf{1}}$	${\mathit{Ra}}_{\mathit{c}\mathbf{2}}^{\mathit{T}}$	${\mathit{\omega }}_{\mathit{c}\mathbf{2}}^{\mathit{T}}$	${\mathit{k}}_{\mathit{c}\mathbf{2}}^{\mathit{T}}$
0.7	354.21	110.819	4.545
0.6	354.31	110.979	4.550
0.5	354.45	111.042	4.550
0.4	354.66	111.251	4.555
0.3	355.05	111.642	4.565
0.2	355.83	112.584	4.590
0.1	358.62	115.209	4.660

gathers the results for seven values of ${\lambda }_1$. It can be observed from Table 2 that critical Rayleigh number $R{a}_{c2}^T$, critical frequency ${\omega }_{c2}^T$ and critical wave number $k_{c2}^T$ at the onset of the secondary instability are weakly dependent on the elasticity number number ${\lambda }_1$. We conclude that the preponderant effect on the properties of the emerging oscillatory TRs is mainly linked to the variations in the viscosity ratio, while the effect of the elasticity remains very weak.

4.3.2. Bifurcation to Steady or Oscillatory Longitudinal Rolls

Finally, we present in the second part of this section the secondary instability results in the case where disturbances are assumed in the form of longitudinal rolls LRs. We mention that as for the primary instability, the onset of stationary LRs convection is not affected by the two viscoelastic parameters. Consequently, the critical Rayleigh number above which stationary LRs convection develops as a secondary instability is the same as that found for Newtonian fluids, namely $R_{c2,s}^L=313.107$. However, the computations indicate Hopf bifurcation from steady unicellular flow to oscillatory LRs convection. We emphasize that the Hopf bifurcation to oscillatory LRs is not observed for Newtonian fluids and is due solely to the viscoelastic character of the fluids. The effects of the two viscoelastic parameters on the linear properties of the oscillatory LRs convection are examined in the remainder of this subsection. In order to evaluate the effect of elasticity alone, the ${\lambda }_1=0.1$, ${\lambda }_1=0.3$ and ${\lambda }_1=0.5$ cases are investigated for a fixed $\mathsf{\Gamma }=0.75$. On the other hand, the effect of viscosity ratio alone is studied by fixing ${\lambda }_1=0.1$ and investigating the $\mathsf{\Gamma }=0.75$, $\mathsf{\Gamma }=0.6$, $\mathsf{\Gamma }=0.5$ and $\mathsf{\Gamma }=0.3$ cases. The computed results for the six different cases are reported in

Table 3. Critical Rayleigh number $R{a}_{c2}^L$, frequency ${\omega }_{c2}^L$ and wave number $k_{c2}^L$ at the onset of oscillatory longitudinal rolls as the secondary instability for different values of $\mathsf{\Gamma }$ and ${\lambda }_1$.

${\mathit{\lambda }}_{\mathbf{1}}$	$\mathsf{\Gamma }$	${\mathit{Ra}}_{\mathit{c}\mathbf{2},\mathbf{osc}}^{\mathit{L}}$	${\mathit{\omega }}_{\mathit{c}\mathbf{2}}^{\mathit{L}}$	${\mathit{k}}_{\mathit{c}\mathbf{2}}^{\mathit{L}}$	${\mathit{R}}_{\mathit{c}\mathbf{2}}^{\mathit{T}}$	${\mathit{R}}_{\mathit{c}\mathbf{2},\mathit{s}}^{\mathit{L}}$
Newtonian		-	-	-	506.27	313.107
0.1	0.75	426.27	1.53	5.8	358.62	313.107
0.3	0.75	317.55	3.58	4.5	355.03	313.107
0.5	0.75	291.34	2.65	3.9	354.45	313.107
0.1	0.6	333.47	12.35	6.3	272.90	313.107
0.1	0.5	288.08	17.53	6.5	219.04	313.107
0.1	0.3	194.20	33.62	7.0	123.55	313.107

, which indicates the critical Rayleigh number, wave number and oscillatory frequency at the onset of oscillatory LRs’ secondary instability. As has already been highlighted in the previous sections considering the primary instability and the TRs’ secondary instability, we recognize the destabilizing effect of the elasticity number ${\lambda }_1$ and the stabilizing effect of the viscosity ratio $\mathsf{\Gamma }$. Moreover, a comparison of Table 1, Table 2 and Table 3 attests that the frequencies of oscillatory LRs are much smaller than those corresponding to oscillatory TRs.

An additional remark about Table 3 is necessary. For comparison purposes, we also indicate in this table the threshold of both stationary and oscillatory LRs. It is clear that the true critical Rayleigh number depends on the combination of the rheological parameters. The least stable mode of convection is the one with smallest critical Rayleigh number and is identified in Table 3 with a bold character. For instance, we consider diluted viscoelastic solutions with $\mathsf{\Gamma }=0.75$ with different elasticity number ${\lambda }_1$. For the combination of the rheological parameters (${\lambda }_1=0.1$, $\mathsf{\Gamma }=0.75$), the true critical Rayleigh number is $R_{c2,s}^L$, indicating that the secondary instability pattern is in the form of steady LRs. In that case, polymeric solutions are almost inelastic and evolve as a Newtonian fluid. In contrast, for the combination (${\lambda }_1=0.5$, $\mathsf{\Gamma }=0.75$), the least stable mode of convection changes from steady LRs to oscillatory LRs, meaning that elastic effects become the most important ones in this range. In the same way, the preferred pattern as a secondary instability depends on the viscosity ratio $\mathsf{\Gamma }$. Table 3 shows that by keeping ${\lambda }_1=0.1$ and decreasing gradually $\mathsf{\Gamma }$ from $\mathsf{\Gamma }=0.75$ (diluted solutions) to $\mathsf{\Gamma }=0.3$ (concentrated solutions), the most amplified mode of convection evolves from steady LRs to oscillatory LRs and eventually to oscillatory TRs.

All of the results stated in the Section 4.3 are obtained by assuming infinite aspects ratios in the x and y directions. For the sake of brevity, we exemplify the effect of the lateral aspect ratio a on the pattern selection for two combinations of rheological parameters ($\mathsf{\Gamma }=0.75$, ${\lambda }_1=0.3$) and ($\mathsf{\Gamma }=0.5$, ${\lambda }_1=0.1$). We plot in Figure 8a,b the variation of the critical Rayleigh number for both stationary LRs and oscillatory LRs as a function of the lateral aspect ratio a in the cases ($\mathsf{\Gamma }=0.75$, ${\lambda }_1=0.3$) and ($\mathsf{\Gamma }=0.5$, ${\lambda }_1=0.1$), respectively. Computations showed that there is a competition between the two patterns in the sense that depending on the magnitude of lateral confinement, the system may select either stationary LRs or oscillatory LRs. For a fixed value of L and by increasing a, the following behavior is observed for the curves representing the critical Rayleigh number for LRs and all values of rheological parameters (see Figure 8a,b): (i) the curve associated with the critical Rayleigh number of oscillatory LRs decreases to reach a minimum equal to its value for infinite a. This minimum point moves to the right in the (a, $Ra$) plane when the number of rolls L is increased; (ii) then, the same curve increases to intersect another branch corresponding to the critical Rayleigh number of steady LRs at a particular value of a; (iii) finally, when a exceeds this particular value, the curve associated with the critical Rayleigh number of steady LRs becomes the lower curve, decreases monotonically and tends asymptotically to the critical Rayleigh number $R_{c2,s}^L=313.107$ of steady LRs found in the case of infinite a.

For the particular combination ($\mathsf{\Gamma }=0.75$, ${\lambda }_1=0.3$) and in the limit of infinite a, the critical Rayleigh number $R_{c2,s}^L=313.107$ of steady LRs is less than the critical Rayleigh number $R_{c2,osc}^L=317.55$ of oscillatory LRs. Consequently, the decreasing curve of the critical Rayleigh number of steady LRs (ultra thick curve) with $L=1$ crosses the absolute minimum $R_{c2,osc}^L=317.55$ of oscillatory LRs at a particular value $a=a^{**}$ ($a^{**}\approx 2$ as shown in Figure 8a). This means that for all values of a larger than $a^{**}$, the dominant mode of convection is a steady monocellular LRs. Otherwise, the system may select oscillatory LRs or steady monocellular LRs depending on a. It is important to note that Figure 8a also shows that the curve corresponding to oscillatory longitudinal mode (thick curve) with $L=1$ intersects the horizontal line representing the critical Rayleigh number of oscillatory TRs $R{a}_{c2}^T=355.03$ at a particular value of the lateral aspect ratio $a=a^{*}\approx 0.4$. This means that perturbations promote the appearance of oscillatory TRs if $a<a^{*}$, oscillatory LRs or a steady monocellular LRs if $a^{*}<a<a^{**}$ and stationary LRs if $a>a^{**}$.

In the case of the combination ($\mathsf{\Gamma }=0.5$, ${\lambda }_1=0.1$), this behavior is not observed, since as can be seen from Figure 8b, the critical Rayleigh number of oscillatory TRs is much smaller than the critical Rayleigh number for both stationary and oscillatory LRs. For this particular combination, the system selects oscillatory TRs independently of the lateral confinement.

5. Conclusions

In the present paper, the Galerkin method was used to investigate the primary and secondary instabilities of viscoelastic fluids saturating a porous layer heated from below by a constant flux. The modified Darcy’s law based on the Oldroyd-B model was used for modeling the momentum equation. In addition to the Darcy–Rayleigh number $Ra$, two viscoelastic parameters play a key role when characterizing the temporal behavior of the instability, namely the relaxation time ${\lambda }_1$, which measures the elasticity of the fluid and the ratio $\mathsf{\Gamma }$ between the viscosity of the solvent and the total viscosity of the fluid. In the first part of the paper, three-dimensional disturbances were considered in order to study the stability of the basic motionless solution. For sufficiently elastic fluids, we found that the primary instability is oscillatory. Otherwise, the primary bifurcation gives rise to stationary long wave instability. Results indicated that the lateral confinement of the porous layer by isolated side walls eliminates oblique or longitudinal rolls in favor of two-dimensional transverse rolls. Based on a fully-developed parallel flow assumption, a nonlinear analytical solution for the velocity and temperature fields was developed in the range of the rheological parameters where stationary long wave instability develops first. In the second part of the paper, we reported findings on the linear stability analysis of the monocellular flow, which is performed with special attention given to the interplay between the viscoelastic parameters and the lateral aspect ratio a of the porous layer. For weakly-elastic fluids, we determined a second critical value of Rayleigh number above which the system exhibits a Hopf bifurcation from steady monocellular flow to oscillatory transverse rolls convection. The well-known limit of $R{a}_{c2}^T\approx 506$ for Newtonian fluids is recovered, and the fluid elasticity effect is found to delay the onset of the Hopf bifurcation.

Three-dimensional analysis showed that for the diluted solutions, as the Boger fluid type (i.e., $\mathsf{\Gamma }=0.75$), the monocellular flow is more unstable for either stationary longitudinal disturbances for weakly-elastic fluids (${\lambda }_1=0.1$) or for oscillatory longitudinal rolls for strongly-elastic fluids (${\lambda }_1=0.5$). This pattern selection holds if the lateral walls are pushed to infinity. When a finite lateral confinement is taken into account, there exist particular values $a^{*}$ and $a^{**}$ of the lateral aspect ratio a such that perturbations promote the appearance of oscillatory transverse rolls if $a<a^{*}$, stationary or oscillatory longitudinal rolls if $a^{*}<a<a^{**}$ and stationary longitudinal rolls if $a>a^{**}$. Computations proved that the interval $[a^{*},a^{**}]$ is enlarged by increasing the fluid elasticity.

For concentrated viscoelastic fluids ($\mathsf{\Gamma }=0.6$, $\mathsf{\Gamma }=0.5$ and $\mathsf{\Gamma }=0.3$), it is found that oscillatory transverse rolls are the preferred mode of convection even for weakly-elastic fluids and independently of the lateral confinement of the porous medium.