Stability of a Buoyant Oldroyd-B Flow Saturating a Vertical Porous Layer with Open Boundaries

Abstract

The performance of several engineering applications are strictly connected to the rheology of the working fluids and the Oldroyd-B model is widely employed to describe a linear viscoelastic behaviour. In the present paper, a buoyant Oldroyd-B flow in a vertical porous layer with permeable and isothermal boundaries is investigated. Seepage flow is modelled through an extended version of Darcy’s law which accounts for the Oldroyd-B rheology. The basic stationary flow is parallel to the vertical axis and describes a single-cell pattern where the cell has an infinite height. A linear stability analysis of such a basic flow is carried out to determine the onset conditions for a multicellular pattern. This analysis is performed numerically by employing the shooting method. The neutral stability curves and the values of the critical Rayleigh number are evaluated for different retardation time and relaxation time characteristics of the fluid. The study highlights the extent to which the viscoelasticity has a destabilising effect on the buoyant flow. For the limiting case of a Newtonian fluid, the known results available in the literature are recovered, namely a critical value of the Darcy–Rayleigh number equal to 197.081 and a corresponding critical wavenumber of 1.05950.

1. Introduction

The convective instability in a plane vertical porous layer saturated by a fluid is a topic of great interest for its applications spanning from geophysical systems to building insulation. Although the analysis of geophysical systems typically involves groundwater or hydrocarbons, the study of thermal insulation panels features air or another gas as a working fluid. Hydrocarbons saturating a porous medium may occasionally display a non-Newtonian viscoelastic behaviour. This circumstance will be that envisaged in the present study.

The classical paper by Gill [1] offered a rigorous mathematical proof that a vertical porous layer, saturated by a Newtonian fluid endowed with impermeable boundaries, kept at different uniform temperatures, displays a stationary conduction regime with a stable parallel vertical buoyant flow, regardless of the temperature gap between the boundaries. In other words, such basic buoyant flow is always stable for every value of the Rayleigh number. The linear stability analysis carried out by Gill [1] is relative to a Newtonian fluid subject to Darcy’s law. In the context of Newtonian fluids, other authors further developed Gill’s important result by including other features such as the non-linearity of the perturbation dynamics or other momentum transfer models of the seepage flow in the porous medium [2,3,4,5,6,7]. Further investigations provided evidence that the stability theorem proved by Gill does not hold when the hypothesis of impermeable boundaries for the porous layer is released [8,9,10,11,12]. As proved in these studies, with perfectly or partly permeable boundaries, the onset of convective instability is possible for sufficiently large values of the Rayleigh number.

Stability analyses of the buoyant flow in a vertical porous layer envisage a non-Newtonian saturating fluid only in a few instances [13,14,15]. Barletta and Alves [13] considered a power-law fluid, while Shankar and Shivakumara [14,15] considered a viscoelastic fluid described according to the Oldroyd-B model as developed, for saturated porous media, by Khuzhayorov et al. [16] and employed by other authors [17,18]. However, the non-Newtonian versions of the stability analysis of Gill’s system developed by Barletta and Alves [13] and by Shankar and Shivakumara [14,15] are relative to a layer with impermeable boundaries.

Following the path devised by Barletta [8], the aim of this paper is to relax the assumption of impermeable boundaries and reconsider the linear stability analysis of the vertical buoyant flow in a vertical layer saturated by a viscoelastic Oldroyd-B fluid. The special case analysed by Barletta [8] is retrieved when the relaxation time of the viscoelastic fluid coincides with the retardation time. On the other hand, the viscoelastic behaviour shows up when such times are different. In this case, the threshold to convective instability is affected by the viscoelastic rheology of the fluid. The neutral stability condition and the critical Rayleigh number will be obtained for different values of the characteristic governing parameters of viscoelasticity.

2. Mathematical Model

A vertical porous layer saturated by a viscoelastic fluid is studied. The layer is infinitely wide in the y and z directions and has thickness L in the horizontal direction x. A sketch of the layer is presented in Figure 1. The Oldroyd-B model is employed to describe the viscoelasticity. Darcy’s law for Oldroyd-B type of fluids [16] is written

(1)$\begin{array}{c}\hfill \begin{array}{c}\hfill \frac{\mu }{K}\left(1+{\tau }_2{\displaystyle \frac{\partial }{\partial t}}\right)\mathbf{u}=-\left(1+{\tau }_1{\displaystyle \frac{\partial }{\partial t}}\right)\left[\mathbf{\nabla }p+{\rho }_f\mathbf{g}\beta \left(T-T_0\right)\right],\end{array}\end{array}$

where ${\tau }_1$ and ${\tau }_2$ are two characteristic relaxation/retardation times, respectively, with ${\tau }_2⩽{\tau }_1$. The buoyancy force is modelled by the Oberbeck–Boussinesq approximation, $\mathbf{x}=\left(x,y,z\right)$ is the Cartesian position vector, $\mathbf{u}=\left(u,v,w\right)$ is the filtration velocity vector, T is the temperature and $\mathbf{g}$ is the gravitational acceleration vector (opposite to the z direction).

The energy balance equation employed to model the heat transfer is the convection/conduction equation where no source or sink term is considered.

The vertical boundaries of the layer are permeable and the external environments in the regions $x<-1/2$ and $x>1/2$ are considered as isothermal fluid reservoirs in a motionless state, so that the pressure distribution along the boundaries is purely hydrostatic. Since the external fluid reservoirs are kept at different uniform temperatures $T_1$ and $T_2$, where $T_2>T_1$, the layer is subject to side heating, as depicted in Figure 1.

Thus, the governing equations are given by

(2)$\begin{array}{c}\hfill \begin{array}{c}\mathbf{\nabla }·\mathbf{u}=0,\hfill \\ \frac{\mu }{K}\left(1+{\tau }_2{\displaystyle \frac{\partial }{\partial t}}\right)\mathbf{u}=-\left(1+{\tau }_1{\displaystyle \frac{\partial }{\partial t}}\right)\left[\mathbf{\nabla }p-{\rho }_{f}g\beta \left(T-T_0\right){\mathbf{e}}_z\right],\hfill \\ \sigma {\displaystyle \frac{\partial T}{\partial t}}+\mathbf{u}·\mathbf{\nabla }T=\chi {\nabla }^{2}T,\hfill \\ x=-L/2:p=p_0,T=T_1,\hfill \\ x=+L/2:p=p_0,T=T_2,\hfill \end{array}\end{array}$

where g and ${\mathbf{e}}_z$ are the modulus of $\mathbf{g}$ and the unit vector along the z direction, respectively, while $p_0$ is the reference hydrostatic pressure of the external environments. In Equation (2) $\mu$ is the dynamic viscosity of the fluid, K is the permeability of the porous medium, ${\rho }_f$ is the fluid density evaluated at the reference temperature $T_0$, $\beta$ is the thermal expansion coefficient of the fluid, $\sigma$ is the heat capacity ratio, $\chi$ is the ratio between the effective thermal conductivity $k_{eff}$ of the fluid saturated porous medium and the volumetric heat capacity of the fluid, $T_0=(T_1+T_2)/2$. The parameters $\sigma$, $k_{eff}$ and $\chi$ are defined as [19]

(3)$\begin{array}{c}\hfill \begin{array}{c}\hfill \sigma ={\displaystyle \frac{\phi {\rho }_f{c}_f+(1-\phi ){\rho }_s{c}_s}{{\rho }_f{c}_f}},k_{eff}=\phi k_f+(1-\phi )k_s,\chi =\frac{k_{eff}}{{\rho }_f{c}_f},\end{array}\end{array}$

where $\phi$ is the porosity, $c_f$ is the specific heat capacity of the fluid, $c_s$ is the specific heat capacity of the porous medium, ${\rho }_s$ is the density of the solid phase, $k_f$ is the thermal conductivity of the fluid and $k_s$ is the thermal conductivity of the solid phase.

The following scaling, where $\Delta T=T_2-T_1$, is employed to obtain the dimensionless formulation of the problem:

(4)$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle \frac{\mathbf{x}}{L}}\to \mathbf{x},{\displaystyle \frac{\chi }{\sigma L^2}}t\to t,{\displaystyle \frac{K}{\chi \mu }}\left(p-p_0\right)\to p,{\displaystyle \frac{L}{\chi }}\mathbf{u}\to \mathbf{u},{\displaystyle \frac{T-T_0}{\Delta T}}\to T,\hfill \\ L\mathbf{\nabla }\to \mathbf{\nabla },L^2{\nabla }^2\to {\nabla }^2.\hfill \end{array}\end{array}$

By substituting Equation (4) in Equation (2), the following dimensionless governing equations are obtained

(5)$\begin{array}{c}\hfill \begin{array}{c}\mathbf{\nabla }·\mathbf{u}=0,\hfill \\ \left(1+{\lambda }_2{\displaystyle \frac{\partial }{\partial t}}\right)\mathbf{u}=-\left(1+{\lambda }_1{\displaystyle \frac{\partial }{\partial t}}\right)\left(\mathbf{\nabla }p-RT{\mathbf{e}}_z\right),\hfill \\ {\displaystyle \frac{\partial T}{\partial t}}+\mathbf{u}·\mathbf{\nabla }T={\nabla }^{2}T,\hfill \\ x=\pm 1/2:p=0,T=\pm 1/2,\hfill \end{array}\end{array}$

where

(6)$\begin{array}{c}\hfill \begin{array}{c}\hfill R={\displaystyle \frac{{\rho }_{f}g\beta \Delta TKL}{\mu \chi }},{\lambda }_{1,2}={\displaystyle \frac{\chi {\tau }_{1,2}}{\sigma L^2}}.\end{array}\end{array}$

Basic State

Let us consider as basic state the steady and fully developed buoyant flow in the vertical direction z and adopt the subscript b to denote the quantities relative to the basic state. Thus, the only non-vanishing component of the filtration velocity vector ${\mathbf{u}}_b$ is $w_b$, along the z-direction. According to Equation (5), ${\mathbf{u}}_b$ results in a solenoidal field. Thus, the basic stationary solution is characterised by a zero net mass flow rate throughout the layer, namely

(7)$\begin{array}{c}\hfill \begin{array}{c}\hfill {\mathbf{u}}_b=(0,0,R{T}_b),p_b=0,T_b=x.\end{array}\end{array}$

The solution for the basic state describes the behaviour of a stationary flow which is parallel to the vertical axis, thus yielding a single-cell vertical pattern where the cell has an infinite vertical height.

3. Linear Stability Analysis

To perform a linear stability analysis, we proceed by perturbing the basic state just defined by applying small-amplitude disturbances, namely

(8)$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathbf{u}={\mathbf{u}}_b+\epsilon \mathbf{U},p=p_b+\epsilon P,T=T_b+\epsilon \mathrm{\Theta },\end{array}\end{array}$

where $\epsilon \ll 1$. By substituting Equation (8) into Equation (5) and by neglecting $O\left({\epsilon }^2\right)$ terms, one obtains

(9)$\begin{array}{c}\hfill \begin{array}{c}\mathbf{\nabla }·\mathbf{U}=0,\hfill \\ \left(1+{\lambda }_2{\displaystyle \frac{\partial }{\partial t}}\right)\mathbf{U}=-\left(1+{\lambda }_1{\displaystyle \frac{\partial }{\partial t}}\right)\left(\mathbf{\nabla }P-R\mathrm{\Theta }{\mathbf{e}}_z\right),\hfill \\ {\displaystyle \frac{\partial \mathrm{\Theta }}{\partial t}}+U+Rx{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial z}}={\nabla }^2\mathrm{\Theta },\hfill \\ x=\pm 1/2:P=0,\mathrm{\Theta }=0.\hfill \end{array}\end{array}$

By applying the divergence operator to the perturbed momentum balance equation and by applying the linear differential operator $(1+{\lambda }_2\partial /\partial t)$ to the perturbed energy balance equation, one obtains the following pressure–temperature formulation, which is more convenient than Equation (9) since the boundary conditions are relative to P and $\mathrm{\Theta }$

(10)$\begin{array}{c}\hfill \begin{array}{c}\left(1+{\lambda }_1{\displaystyle \frac{\partial }{\partial t}}\right)\left({\nabla }^{2}P-R{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial z}}\right)=0,\hfill \\ \left(1+{\lambda }_2{\displaystyle \frac{\partial }{\partial t}}\right)\left({\displaystyle \frac{\partial \mathrm{\Theta }}{\partial t}}+Rx{\displaystyle \frac{\partial \mathrm{\Theta }}{\partial z}}-{\nabla }^2\mathrm{\Theta }\right)-\left(1+{\lambda }_1{\displaystyle \frac{\partial }{\partial t}}\right){\displaystyle \frac{\partial P}{\partial x}}=0,\hfill \\ x=\pm 1/2:P=0,\mathrm{\Theta }=0.\hfill \end{array}\end{array}$

We assume the pressure and temperature perturbations to have the form of normal modes,

(11)$\begin{array}{c}\hfill \begin{array}{c}\{P(x,y,z,t),\mathrm{\Theta }(x,y,z,t)\}=\{f\left(x\right),h\left(x\right)\}{e}^{(\eta -i\omega )t}{e}^{i(k_{y}y+k_{z}z)},\hfill \end{array}\end{array}$

where $\eta$ is the growth or decay rate, $\omega$ is the angular frequency and $\mathbf{k}=(0,k_y,k_z)$ is the wave vector. The parameters $(k_y,k_z,\eta ,\omega )$ are real while $(f,h)$ are, in general, complex functions. The growth rate $\eta$ marks the difference between stability $(\eta <0)$ and instability $(\eta >0)$. The  neutrally stable configuration is identified by $\eta =0$. The condition of minimum R among the neutrally stable modes defines the critical values $k_c$ and $R_c$.

By substituting definitions (11) into Equation (10), one obtains the following eigenvalue problem for neutrally stable modes,

(12)$\begin{array}{c}\hfill \begin{array}{c}{f}^{\prime \prime }-k^{2}f-ikSh=0,\hfill \\ h^{\prime \prime }-\left(k^2+ikSx-i\omega \right)h+\frac{{\lambda }_1\omega +i}{\gamma {\lambda }_1\omega +i}{f}^{\prime }=0,\hfill \\ x=\pm 1/2:f=0,h=0,\hfill \end{array}\end{array}$

where

(13)$\begin{array}{c}\hfill \begin{array}{c}\hfill \gamma =\frac{{\lambda }_2}{{\lambda }_1},k^2=k_y^2+k_z^2,S=\frac{k_z}{k}R.\end{array}\end{array}$

We mention that $0<\gamma ⩽1$ where, according to Equation (1), the limiting case $\gamma =1$ represents the case of a Newtonian fluid.

Moreover, the rescaled Darcy–Rayleigh number S accounts for the inclination of the wave vector to the vertical z axis. When $k_y=0$, one has a transverse mode propagating along the z axis and $S=R$. On the other hand, when $k_z=0$, one has a longitudinal mode propagating along the y axis and $S=0$. Since the latter result can be obtained also when $R=0$, horizontally propagating modes are equivalent to transverse modes having a vanishing Darcy–Rayleigh number. Thus, modes having $k_z<k$ and a given value of S correspond to transverse modes $(k_z=k)$ with the same S and a larger value of R. As a consequence, the transverse modes necessarily are the most unstable.

4. Numerical Solution and Discussion of the Results

The eigenvalue problem formulated through the dimensionless governing Equation (12) can be solved numerically by means of the shooting method. In the considered approach, the associated initial value problem is solved with the Runge–Kutta method, while Brent’s method is adopted for the root searching procedure. The numerical solution of the problem is here found through the built-in functions NDSolve and FindRoot available within the Mathematica 12.0 environment.

In detail, the numerical approach used to solve the eigenvalue problem given by Equation (12) requires the input values $({\lambda }_1,\gamma )$ and determines both $S\left(k\right)$ and $\omega \left(k\right)$ as eigenvalue data. The minimum of $S\left(k\right)$ yields the critical values $(k_c,S_c,{\omega }_c)$ that define the onset of the convective instability.

The main requirement for an efficient solution of Equation (12) through the shooting method is an accurate initial guess for the eigenvalue quantities S and $\omega$, for a given prescribed input dataset $({\lambda }_1,\gamma )$.

The investigated cases refer to wide ranges of values of the dimensionless parameters ${\lambda }_1$ and $\gamma$. In detail, the considered cases are those obtained from the possible pairs $({\lambda }_1,\gamma )$ among the values ${\lambda }_1=(0.015,0.025,0.05,$$0.075,0.1,0.25,0.5)$ and $\gamma =(0.1,0.25,0.5,0.75,1)$. It is worth noticing that these values include the limiting case of Newtonian fluid $(\gamma =1)$ and the special case of Boger’s fluid $({\lambda }_1=0.1,\gamma =0.75)$ [20].

4.1. Limiting Case of a Newtonian Fluid

As already emphasised, the limiting case of a Newtonian fluid occurs when $\gamma =1$. In this case, regardless of the value of ${\lambda }_1$, the neutral stability curve, displayed in Figure 2, shows up a minimum in the plane $(k,S)$ for the critical values

(14)$\begin{array}{c}\hfill \begin{array}{c}{k}_c=1.05950,S_c=197.081,{\omega }_c=0,\hfill \end{array}\end{array}$

which hold for every choice of ${\lambda }_1={\lambda }_2$. It is worth noticing that, as it can be inferred from Equation (12), such a neutral stability solution with $\omega =0$ is always present. In other words, among all the possible values assigned to $({\lambda }_1,\gamma )$, Equation (14) yields the maximum of $S_c$.

4.2. Results for Viscoelastic Fluids

The neutral stability curves for the considered data $({\lambda }_1,\gamma )$ are reported in Figure 3. Each frame refers to a different value of $\gamma$, while the solid lines are drawn for given values of ${\lambda }_1$. For decreasing values of ${\lambda }_1$, the curves move upward, as well as for increasing values of $\gamma$. It is worth noticing that the vertical axis has the range $0⩽S⩽197.081$, since this value corresponds to the onset of instability for a Newtonian fluid, namely the maximum possible critical value of S.

The threshold values of $k,S,\omega$ that have been found numerically are reported in Table 1. The obtained results $(k_c,S_c)$ are insensitive to the sign of ${\omega }_c$. Thus, in Table 1, the values of ${\omega }_c$ are to be intended as $|{\omega }_c|$. As can be inferred from Table 1, for a given value of $\gamma$, an increase in the dimensionless relaxation time ${\lambda }_1$ leads to a decrease in the threshold values $k_c,S_c,{\omega }_c$ (at least, whenever ${\omega }_c\ne 0$). The minimum $(k_c,S_c)$ of the neutral stability curve moves to the left and becomes smaller.

5. Conclusions

The effects of viscoelasticity on the onset of convection in a fluid saturated vertical porous layer have been analysed. An extended version of Darcy’s law has been utilised in order to implement the Oldroyd-B rheological model for the seepage flow in the porous medium. The boundaries of the vertical porous layer have been considered as permeable to external fluid reservoirs kept at different temperatures. The boundary temperature difference causes a side-heating mechanism and, hence, a stationary buoyant flow driven by the uniform horizontal temperature gradient. Such conditions define the basic conduction state. The linear stability of this basic state has been studied by employing a three-dimensional normal mode scheme, with a suitable Squire transformation mapping all possible modes to equivalent transverse rolls. Finally, the stability eigenvalue problem has been formulated and solved numerically by employing the shooting method.

The main highlights of our study are the following:

• The dimensionless governing parameters identifying the viscoelastic behaviour are the relaxation parameter, ${\lambda }_1$, and the ratio between the retardation time and the relaxation time, $\gamma$. The physically significant domain is one where $0<\gamma ⩽1$. When $\gamma \to 1$, the Newtonian fluid behaviour is recovered: the critical value of the modified Darcy–Rayleigh number is 197.081 and the corresponding critical value of the wavenumber is 1.05950;

• The neutral stability threshold to the convective instability obtained for the Newtonian case always exists, whatever is the choice of ${\lambda }_1$ and $\gamma$. However, in most cases, the Newtonian branch of neutral stability is not the lowest one. In these cases, the neutral stability condition is characterised by travelling modes, i.e., modes with a non-zero angular frequency. Thus, the effect of viscoelasticity is generally destabilising with respect to the Newtonian fluid case;

• There exist input values of $({\lambda }_1,\gamma )$, with $\gamma <1$, such that the critical conditions for the onset of the convective instability coincide with those for a Newtonian fluid or, equivalently, the Newtonian branch of neutral stability is the lowest one.

There are several possible directions where the study presented in this paper can be extended. In Authors’ opinion, one of the most interesting developments of this study is the introduction of a model for partially permeable boundary conditions. With such a model one can investigate the gradual transition from permeable to perfectly impermeable boundaries, and its effects on the onset conditions for the instability.

Author Contributions

Conceptualization: S.L., M.C., and A.B.; Formal analysis: S.L., M.C., and A.B.; Funding acquisition: A.B.; Investigation: S.L., M.C., and A.B.; Methodology: S.L., M.C., and A.B.; Project administration: A.B.; Software: S.L.; Supervision: A.B.; Validation: S.L. and M.C.; Writing—original draft: S.L., M.C., and A.B.; Writing—review and editing: S.L., M.C., and A.B. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Italian Ministry of Education and Scientific Research grant number PRIN 2017F7KZWS.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

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

Figures and Table

Figure 1. A sketch of the vertical porous layer.

Figure 2. Neutral stability curve [Forumla omitted. See PDF.] for a Newtonian fluid.

Figure 3. Neutral stability curves [Forumla omitted. See PDF.] for different values of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Table 1

Critical values of $$(k,S,\omega )$$ versus $${\lambda }_1,\gamma$$.

	$$\gamma =0.1$$			$$\gamma =0.25$$
$${\lambda }_1$$	$$k_c$$	$$S_c$$	$${\omega }_c$$	$$k_c$$	$$S_c$$	$${\omega }_c$$
0.015	6.23219	125.630	301.680	4.35669	149.556	166.030
0.025	5.48899	92.0028	186.630	3.77358	114.751	95.0600
0.05	4.56254	59.4450	96.4569	3.24560	78.3898	47.4367
0.075	4.25166	45.9027	71.5633	3.11482	63.1441	35.6140
0.1	4.05828	38.3924	59.2173	3.03338	54.7411	29.7598
0.25	3.50236	23.0465	34.1704	2.77928	37.6406	17.7300
0.5	3.16277	17.0453	23.0225	2.61721	31.0081	12.2442
	$$\gamma =0.5$$			$$\gamma =0.75$$
$${\lambda }_1$$	$$k_c$$	$$S_c$$	$${\omega }_c$$	$$k_c$$	$$S_c$$	$${\omega }_c$$
0.015	1.05950	197.081	0	1.05950	197.081	0
0.025	2.43604	183.662	53.9920	1.05950	197.081	0
0.05	2.21927	125.192	25.2047	1.05950	197.081	0
0.075	2.20008	102.173	19.2658	1.50760	183.721	8.72516
0.1	2.18714	90.0757	16.3564	1.53835	159.824	8.11327
0.25	2.12657	66.8434	10.1776	1.58390	118.328	5.94810
0.5	2.08008	58.4180	7.19023	1.59069	104.712	4.42120

References

1. Gill, A.E. A proof that convection in a porous vertical slab is stable. J. Fluid Mech.; 1969; 35, pp. 545-547. [DOI: https://dx.doi.org/10.1017/S0022112069001273]

2. Rees, D.A.S. The stability of Prandtl–Darcy convection in a vertical porous layer. Int. J. Heat Mass Transf.; 1988; 31, pp. 1529-1534. [DOI: https://dx.doi.org/10.1016/0017-9310(88)90260-8]

3. Straughan, B. A nonlinear analysis of convection in a porous vertical slab. Geophys. Astrophys. Fluid Dyn.; 1988; 42, pp. 269-275. [DOI: https://dx.doi.org/10.1080/03091928808213611]

4. Lewis, S.; Bassom, A.P.; Rees, D.A.S. The stability of vertical thermal boundary–layer flow in a porous medium. Eur. J. Mech. B Fluids; 1995; 14, pp. 395-407.

5. Rees, D.A.S.; Lage, J.L. The effect of thermal stratification on natural convection in a vertical porous insulation layer. Int. J. Heat Mass Transf.; 1996; 40, pp. 111-121. [DOI: https://dx.doi.org/10.1016/S0017-9310(96)00060-9]

6. Rees, D.A.S. The effect of local thermal nonequilibrium on the stability of convection in a vertical porous channel. Transp. Porous Media; 2011; 87, pp. 459-464. [DOI: https://dx.doi.org/10.1007/s11242-010-9694-5]

7. Scott, N.L.; Straughan, B. A Nonlinear stability analysis of convection in a porous vertical channel including local thermal nonequilibrium. J. Math. Fluid Mech.; 2013; 15, pp. 171-178. [DOI: https://dx.doi.org/10.1007/s00021-012-0109-y]

8. Barletta, A. A proof that convection in a porous vertical slab may be unstable. J. Fluid Mech.; 2015; 770, pp. 273-288. [DOI: https://dx.doi.org/10.1017/jfm.2015.154]

9. Celli, M.; Barletta, A.; Rees, D.A.S. Local thermal non–equilibrium analysis of the instability in a vertical porous slab with permeable sidewalls. Transp. Porous Media; 2017; 119, pp. 539-553. [DOI: https://dx.doi.org/10.1007/s11242-017-0897-x]

10. Barletta, A.; Celli, M.; Ouarzazi, M.N. Unstable buoyant flow in a vertical porous layer with convective boundary conditions. Int. J. Therm. Sci.; 2017; 120, pp. 427-436. [DOI: https://dx.doi.org/10.1016/j.ijthermalsci.2017.05.028]

11. Barletta, A.; Rees, D.A.S. On the onset of convection in a highly permeable vertical porous layer with open boundaries. Phys. Fluids; 2019; 31, 074106. [DOI: https://dx.doi.org/10.1063/1.5110484]

12. Barletta, A.; Celli, M. Anisotropy and the onset of the thermoconvective instability in a vertical porous layer. ASME J. Heat Transf.; 2021; 143, 102601. [DOI: https://dx.doi.org/10.1115/1.4051322]

13. Barletta, A.; de Alves, L.S.B. On Gill’s stability problem for non–Newtonian Darcy’s flow. Int. J. Heat Mass Transf.; 2014; 79, pp. 759-768. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2014.08.051]

14. Shankar, B.M.; Shivakumara, I.S. On the stability of natural convection in a porous vertical slab saturated with an Oldroyd–B fluid. Theor. Comput. Fluid Dyn.; 2017; 31, pp. 221-231. [DOI: https://dx.doi.org/10.1007/s00162-016-0415-8]

15. Shankar, B.M.; Shivakumara, I.S. Effect of local thermal nonequilibrium on the stability of natural convection in an Oldroyd–B fluid saturated vertical porous layer. ASME J. Heat Transf.; 2017; 139, 041001. [DOI: https://dx.doi.org/10.1115/1.4035199]

16. Khuzhayorov, B.; Auriault, J.L.; Royer, P. Derivation of macroscopic filtration law for transient linear viscoelastic fluid flow in porous media. Int. J. Eng. Sci.; 2000; 38, pp. 487-504. [DOI: https://dx.doi.org/10.1016/S0020-7225(99)00048-8]

17. Hirata, S.C.; Ouarzazi, M.N. Three–dimensional absolute and convective instabilities in mixed convection of a viscoelastic fluid through a porous medium. Phys. Lett. A; 2010; 374, pp. 2661-2666. [DOI: https://dx.doi.org/10.1016/j.physleta.2010.04.039]

18. de Alves, L.S.B.; Barletta, A.; Hirata, S.; Ouarzazi, M.N. Effects of viscous dissipation on the convective instability of viscoelastic mixed convection flows in porous media. Int. J. Heat Mass Transf.; 2014; 70, pp. 586-598. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.11.041]

19. Nield, D.A.; Bejan, A. Convection in Porous Media; 5th ed. Springer: New York, NY, USA, 2017.

20. Boger, D.V. Viscoelastic flows through contractions. Annu. Rev. Fluid Mech.; 1987; 19, pp. 157-182. [DOI: https://dx.doi.org/10.1146/annurev.fl.19.010187.001105]