Rayleigh–Bénard Instability of an Ellis Fluid Saturated Porous Channel with an Isoflux Boundary

Abstract

The onset of the thermal instability is investigated in a porous channel with plane parallel boundaries saturated by a non–Newtonian shear–thinning fluid and subject to a horizontal throughflow. The Ellis model is adopted to describe the fluid rheology. Both horizontal boundaries are assumed to be impermeable. A uniform heat flux is supplied through the lower boundary, while the upper boundary is kept at a uniform temperature. Such an asymmetric setup of the thermal boundary conditions is analysed via a numerical solution of the linear stability eigenvalue problem. The linear stability analysis is developed for three–dimensional normal modes of perturbation showing that the transverse modes are the most unstable. The destabilising effect of the non–Newtonian shear–thinning character of the fluid is also demonstrated as compared to the behaviour displayed, for the same flow configuration, by a Newtonian fluid.

1. Introduction

The emergence of convection Rayleigh–Bénard cells in a fluid saturated porous medium heated from below is a cornerstone topic in the research on convection heat transfer in porous solids. There is an impressively wide literature on this topic, often called Darcy–Bénard instability, which is surveyed in Chapter 6 of Nield and Bejan [1]. Other surveys relative to this area of research are provided by Straughan [2], Barletta [3]. Except for a minor part of the papers available in the literature regarding the Darcy–Bénard instability, most investigators have focussed their attention on the Newtonian fluids saturating a porous medium.

The Darcy–Bénard instability for saturating fluids with a non–Newtonian rheology has been investigated with reference to viscoelastic fluids [4,5,6], to fluids with yield–stress [7] and to purely viscous fluids [8,9,10,11,12,13,14]. In particular, Barletta and Nield [8], Alves and Barletta [9], Celli et al. [10] and Petrolo et al. [11] discuss the onset of the Darcy–Bénard instability in a porous medium by considering the power–law model. A more sophisticated model, i.e. Carreau–Yasuda, is employed by Brandão and Ouarzazi [12] and by Brandão et al. [13] to model purely viscous, non–Newtonian flow in porous media.

The usual setup envisaged in the studies on the onset of Darcy–Bénard instability in a saturated porous layer is one where the plane parallel boundaries are considered as isothermal with different temperatures. In the present paper, the objective is to extend the analysis carried out by Celli et al. [14]. In fact, Celli et al. [14] employ Ellis’ model to describe the behaviour of a shear–thinning non–Newtonian fluid saturating a porous medium. The temperature boundary conditions are modified, with respect to the study developed by Celli et al. [14], by assuming a uniform heat flux on the lower boundary of the layer. Such a modification is accompanied by a numerical solution of the stability eigenvalue problem, whereas an analytical dispersion relation is available with isothermal boundary conditions [14]. The shear–thinning rheology of the fluid is formulated by adopting the Ellis model extension of Darcy’s law. Results will be discussed for the instability onset in terms of the neutral stability curves and of the critical values of the perturbation wavenumber and of the Darcy–Rayleigh number. The main result is the demonstration of the destabilising effect of the Ellis fluid rheology as compared to the Newtonian fluid behaviour, with the transverse normal modes leading the transition to the instability.

2. Mathematical Formulation

A porous layer of height H saturated by a non–Newtonian fluid is considered. The boundary walls are considered to be impermeable. The layer is heated from below by a constant heat–flux $q_0$, while the upper boundary is kept at a constant temperature $T_0$. A sketch of the porous layer with a scheme of the boundary conditions and of the coordinate axes is displayed in Figure 1. The Oberbeck–Boussinesq approximation is employed and local thermal equilibrium between the fluid and the solid matrix is assumed. A basic uniform throughflow is imposed along the horizontal x direction.

2.1. Rheological Model

The rheological behaviour of the non–Newtonian fluid is described by employing the Ellis model. More precisely, such a model is based on three–parameters that describe the rheology of a time–independent, shear–thinning and non–yield–stress fluid [15]. According to Sochi [15], the Ellis model is more reliable than the power–law model in matching the experimental data. The apparent viscosity $\eta$, according to the Ellis model, is given by [15,16].

(1)$\begin{array}{c}\begin{array}{c}\eta ={\displaystyle \frac{{\eta }_0}{1+{\left({\displaystyle \frac{\tau }{{\tau }_{1/2}}}\right)}^{\zeta -1}}},\end{array}\end{array}$

where $\zeta$ is a positive parameter such that $\zeta =1/n$, where n is the power–law index. The apparent viscosity at zero shear–stress is represented by ${\eta }_0$, while ${\tau }_{1/2}$ represents the value of $\tau$ at which $\eta ={\eta }_0/2$. By employing the power–law index, Equation (2) can be rewritten as

(2)$\begin{array}{c}\begin{array}{c}\eta ={\displaystyle \frac{{\eta }_0}{1+{\left({\displaystyle \frac{\tau }{{\tau }_{1/2}}}\right)}^{\frac{1-n}{n}}}}.\end{array}\end{array}$

The apparent viscosity predicted by the Ellis model, for some limiting cases, is given by the following list:

(3)$\begin{array}{c}\left\{\begin{array}{cc}{\tau }_{1/2}\to 0\hfill & ⟹\eta =0,\hfill \\ {\tau }_{1/2}\to \infty \hfill & ⟹\eta ={\eta }_0,\hfill \\ \tau <{\tau }_{1/2},n\to 0\hfill & ⟹\eta ={\eta }_0,\hfill \\ \tau >{\tau }_{1/2},n\to 0\hfill & ⟹\eta =0,\hfill \\ n\to 1\hfill & ⟹\eta ={\eta }_0/2.\hfill \end{array}\right.\end{array}$

On the other hand, in the limiting case of strong shear stresses, $\tau \gg {\tau }_{1/2}$, Equation (2) reduces to

(4)$\begin{array}{c}\begin{array}{c}\eta ={\eta }_0{\left({\displaystyle \frac{\tau }{{\tau }_{1/2}}}\right)}^{\frac{n-1}{n}}.\end{array}\end{array}$

2.2. Generalization of Darcy’s Law

For a Newtonian fluid saturating a porous medium, Darcy’s law can be formulated as [1]

(5)$\begin{array}{c}\begin{array}{c}\overline{\mathbf{u}}={\displaystyle \frac{K}{\eta }}{\overline{\mathbf{f}}}_d,\end{array}\end{array}$

where the bars over the symbols represent the dimensional fields, time and coordinates. Here, $\overline{\mathbf{u}}$ is the seepage velocity vector with components $(\overline{u},\overline{v},\overline{w})$, K is the permeability of the porous medium and ${\overline{\mathbf{f}}}_d$ is the drag force which contains the pressure gradient and the buoyancy term based on the Oberbeck–Boussinesq approximation, namely [1]

(6)$\begin{array}{c}\begin{array}{c}{\overline{\mathbf{f}}}_d=-\mathbf{\nabla }\overline{p}-{\rho }_0\mathbf{g}\beta (\overline{T}-T_0).\end{array}\end{array}$

In Equation (6), $\overline{p}$ is the local difference between the pressure and the hydrostatic pressure, ${\rho }_0$ is the fluid density evaluated at the reference temperature $T_0$, $\mathbf{g}$ is the gravity acceleration vector and $\beta$ is the thermal expansion coefficient of the fluid. Darcy’s law can be generalised for non–Newtonian fluids by replacing the fluid viscosity with an effective viscosity [16,17]. Hence, Equation (5) becomes

(7)$\begin{array}{c}\begin{array}{c}\overline{\mathbf{u}}={\displaystyle \frac{K}{{\eta }_{eff}}}{\overline{\mathbf{f}}}_d,\end{array}\end{array}$

where ${\eta }_{eff}$ is the effective viscosity given by

(8)$\frac{1}{{\eta }_{eff}}=\frac{1}{{\eta }_0}\left[1+\frac{4n}{3n+1}{\left(\frac{|{\overline{\mathbf{f}}}_d|r_h}{{\tau }_{1/2}}\right)}^{\frac{1-n}{n}}\right].$

Here, $r_h$ is the mean hydraulic radius, which is directly proportional to the square root of the permeability of the porous medium divided by its porosity $\phi$ [16]

(9)$r_h=\sqrt{\frac{5K}{\phi }}.$

Thus, the modified Darcy’s law for a porous medium saturated by an Ellis fluid can be rewritten in a more convenient way as

(10)$\begin{array}{c}\begin{array}{c}\overline{\mathbf{u}}={\displaystyle \frac{K}{{\eta }_0}}\left(1+A|{\overline{\mathbf{f}}}_d{|}^{\frac{1-n}{n}}\right){\overline{\mathbf{f}}}_d,\end{array}\end{array}$

where A is a coefficient determined by the properties of the porous medium and of the fluid,

(11)$\begin{array}{c}\begin{array}{c}A=\frac{4n}{3n+1}{\left(\frac{r_h}{{\tau }_{1/2}}\right)}^{\frac{1-n}{n}}.\end{array}\end{array}$

It is worth noting that, near the instability threshold, the drag forces have, usually, a small magnitude. In fact, in the limit of vanishing drag forces, ${\overline{\mathbf{f}}}_d\to 0$, Equation (10) reduces to Equation (5).

2.3. Governing Equations

The mass balance, momentum balance and energy balance equations for the present problem are given by

(12)$\begin{array}{c}\begin{array}{c}\mathbf{\nabla }·\overline{\mathbf{u}}=0,\\ \frac{{\eta }_0}{K}\overline{\mathbf{u}}=\left(1+A|{\overline{\mathbf{f}}}_d{|}^{\frac{1-n}{n}}\right){\overline{\mathbf{f}}}_d,\\ {\overline{\mathbf{f}}}_d=-\mathbf{\nabla }\overline{p}-{\rho }_0\mathbf{g}\beta (\overline{T}-T_0),\\ \sigma \frac{\partial \overline{T}}{\partial \overline{t}}+\overline{\mathbf{u}}·\mathbf{\nabla }\overline{T}=\alpha {\nabla }^2\overline{T},\\ \overline{y}=0:\overline{v}=0,-\chi \frac{\partial \overline{T}}{\partial \overline{y}}=q_0,\\ \overline{y}=H:\overline{v}=0,\overline{T}=T_0.\end{array}\end{array}$

In Equation (12), $\sigma$ is the ratio between the average volumetric heat capacity of the porous medium and the volumetric heat capacity of the fluid; $\alpha$ and $\chi$ are the average thermal diffusivity and the average thermal conductivity of the saturated porous medium, respectively.

The following scaling allows one to express Equation (12) in a dimensionless form

(13)$\begin{array}{c}\begin{array}{c}\mathbf{x}=\frac{\overline{\mathbf{x}}}{H},\mathbf{u}={\displaystyle \frac{H}{\alpha }}\overline{\mathbf{u}},p={\displaystyle \frac{K}{{\eta }_0\alpha }}\overline{p},t={\displaystyle \frac{\alpha }{\sigma H^2}}\overline{t},T={\displaystyle \frac{\overline{T}-T_0}{\Delta T}},\end{array}\end{array}$

where $\Delta T=q_{0}H/\chi$ and $\mathbf{x}$ is the position vector with Cartesian components $(x,y,z)$. By substituting Equation (13) into Equation (12), one may write

(14)$\begin{array}{c}\begin{array}{c}\mathbf{\nabla }·\mathbf{u}=0,\\ \mathbf{u}=\left(1+\mathrm{El}|{\mathbf{f}}_d{|}^{\frac{1-n}{n}}\right){\mathbf{f}}_d,\\ \frac{\partial T}{\partial t}+\mathbf{u}·\mathbf{\nabla }T={\nabla }^{2}T,\\ y=0:v=0,-\frac{\partial T}{\partial y}=1,\\ y=1:v=0,T=0,\end{array}\end{array}$

where

(15)$\begin{array}{c}\begin{array}{c}{\mathbf{f}}_d=-\mathbf{\nabla }p+\mathrm{R}T{\mathbf{e}}_y.\end{array}\end{array}$

The parameter $\mathrm{El}$ is the Darcy–Ellis number and the parameter $\mathrm{R}$ is the Darcy–Rayleigh number. They are defined as follows:

(16)$\begin{array}{c}\mathrm{El}=A{\left({\displaystyle \frac{\alpha {\eta }_0}{HK}}\right)}^{\frac{1-n}{n}},\mathrm{R}={\displaystyle \frac{\rho g\beta HK\Delta T}{\alpha {\eta }_0}}.\end{array}$

2.4. Basic State

The stationary solution of the problem is given by a uniform throughflow in the horizontal direction and by a negative vertical temperature gradient. The latter is due to the supplied heat flux on the bottom and the prescribed temperature on the top. The uniform throughflow is generated by a prescribed constant pressure drop along the x axis, $\partial p_b/\partial x$. Thus, the basic solution is expressed as

(17)$\begin{array}{c}\begin{array}{c}{u}_b=-\frac{\partial p_b}{\partial x}\left(1+\mathrm{El}{\left|\frac{\partial p_b}{\partial x}\right|}^{\frac{1-n}{n}}\right),v_b=0,\\ w_b=0,\frac{\partial p_b}{\partial y}=\mathrm{R}{T}_b,\frac{\partial p_b}{\partial z}=0,T_b=1-y,\end{array}\end{array}$

where the subscript b denotes the basic state. Without any loss of generality, we can assume that $\partial p_b/\partial x$ is negative, so that $u_b$ assumes only positive values.

2.5. Linear Stability Analysis

In order to simplify the numerical treatment, the governing equations can be rewritten by employing a pressure–temperature formulation. By considering the first Equation (14) and by applying the divergence operator to the second Equation (14), we can get the rid of the velocity field in the momentum balance equation. For the energy balance equation, we employ the second Equation (14) to express the velocity as function of pressure. Thus, the pressure–temperature formulation of the governing equations is given by

(18)$\begin{array}{c}\begin{array}{c}\mathbf{\nabla }·\left[\left(1+\mathrm{El}|{\mathbf{f}}_d{|}^{\frac{1-n}{n}}\right){\mathbf{f}}_d\right]=0,\\ \frac{\partial T}{\partial t}+\left[\left(1+\mathrm{El}|{\mathbf{f}}_d{|}^{\frac{1-n}{n}}\right){\mathbf{f}}_d\right]·\mathbf{\nabla }T={\nabla }^{2}T,\\ {\mathbf{f}}_d=-\mathbf{\nabla }p+\mathrm{R}T{\mathbf{e}}_y,\\ y=0:\frac{\partial p}{\partial y}=\mathrm{R}T,-\frac{\partial T}{\partial y}=1,\\ y=1:\frac{\partial p}{\partial y}=0,T=0,\end{array}\end{array}$

where the boundary conditions based on the pressure–temperature formulation are obtained by employing the governing Equation (14).

2.6. Normal Modes

The stability of the system is investigated by decomposing pressure and temperature into two parts: one relative to the basic equilibrium solution, and the other relative to the infinitesimal disturbances. The stability analysis consists in observing the linear evolution in time of such disturbances in order to determine the threshold values of the governing parameters for the onset of the instability. Such disturbances are investigated through a normal mode analysis. The decomposed pressure and temperature are given by

(19)$\begin{array}{c}\begin{array}{c}p(x,y,z,t)=p_b(x,y)+\epsilon f\left(y\right)e^{\lambda t}{e}^{i\left(k_{x}x+k_{z}z-\omega t\right)},\\ T(x,y,z,t)=T_b\left(y\right)+\epsilon h\left(y\right)e^{\lambda t}{e}^{i\left(k_{x}x+k_{z}z-\omega t\right)},\end{array}\end{array}$

where $\lambda$ is the growth rate, $k_x$ and $k_z$ are the wavenumbers in the x and z direction, $\omega$ is the angular frequency, f and h are the eigenfunctions of the problem, and $\epsilon$ is a positive parameter that defines the amplitude of the disturbances. Since $\epsilon \ll 1$, all terms $\mathcal{O}\left({\epsilon }^2\right)$ are neglected in the present analysis.

The instability threshold is determined by seeking the neutral stability, which is the transitional condition given by the normal modes with a zero growth rate. The present analysis is focussed on the dynamics of the most unstable normal modes arising in the system. By imposing $\lambda =0$, we ensure the neutral stability condition. Hence, the differential dispersion relation for the disturbances is given by

(20)$\begin{array}{c}{\tilde{f}}^{″}-\tilde{n}{k}^2\tilde{f}-\tilde{\mathrm{R}}{h}^{\prime }=0,\\ h^{″}-(k^2-\tilde{\mathrm{R}}-i\tilde{\omega })h-{\tilde{f}}^{\prime }=0,\end{array}$

subject to the boundary conditions

(21)$\begin{array}{c}y=0:{\tilde{f}}^{\prime }=\tilde{\mathrm{R}}h,h^{\prime }=0,\\ y=1:{\tilde{f}}^{\prime }=0,h=0,\end{array}$

where the following quantities have been defined:

(22)$\begin{array}{c}\begin{array}{c}\tilde{f}=(1+\widetilde{\mathrm{El}})f,\tilde{\mathrm{R}}=(1+\widetilde{\mathrm{El}})\mathrm{R},\widetilde{\mathrm{El}}=\mathrm{El}{\left|\frac{\partial p_b}{\partial x}\right|}^{\frac{1-n}{n}},\\ k_x=kcos\varphi ,k_z=ksin\varphi ,\\ \tilde{n}={\displaystyle \frac{\widetilde{\mathrm{El}}+n(\widetilde{\mathrm{El}}+2)+\widetilde{\mathrm{El}}(1-n)cos\left(2\varphi \right)}{2n(\widetilde{\mathrm{El}}+1)}}.\end{array}\end{array}$

Here, $\varphi$ is the angle between the wave vector of the normal mode and the x axis. In Appendix A, a proof that the eigenvalue problem is to be solved with $\tilde{\omega }=0$ is provided. In Appendix B, it is shown that the most unstable modes are the transverse rolls, whose wave vector is parallel to the basic flow direction $(\varphi =0)$. These findings are the basis for the discussion of the results presented in the forthcoming section.

3. Results and Discussion

Equations (20) together with the boundary conditions (21) are solved numerically by using the procedure described in Appendix C. The results of the stability analysis are presented by displaying the neutral stability curves and the plots of the critical values of $\mathrm{R}$ and k versus the governing parameters ($\widetilde{\mathrm{El}},n$).

In Figure 2, the neutral threshold is displayed in the $(k,\mathrm{R})$ plane for different values of $\widetilde{\mathrm{El}}$ and $n=0.5$. It is evident that the parameter $\widetilde{\mathrm{El}}$ has a destabilising effect. In fact, with any fixed wavenumber k, the neutral stability condition for increasing values of $\widetilde{\mathrm{El}}$ is achieved for decreasing values of $\mathrm{R}$.

The most important information of a neutral stability curve is the global minimum of $\mathrm{R}$, which defines the critical condition. In order to compute directly the critical points, one can derive the original Equation (20), as well as the boundary conditions (21), with respect to k and impose the condition $\partial \mathrm{R}/\partial k=0$. The derivative of the eigenfunctions $\tilde{f}$ and h with respect to k are treated as new eigenfunctions and the differential system expanded in this way is solved together with the original one. More information regarding this technique can be found in the paper by Alves et al. [18] and in the book by Barletta [3].

Figure 3 and Figure 4 show the critical stability condition by reporting ${\mathrm{R}}_c$ and $k_c$ as functions of $\widetilde{\mathrm{El}}$ for different values of n. By analysing these results, the destabilising role of $\widetilde{\mathrm{El}}$ is confirmed. In addition, one may note that, for a vanishing $\widetilde{\mathrm{El}}$, the critical data do not depend on n. Indeed, a vanishing value of $\widetilde{\mathrm{El}}$ means either $\mathrm{El}\to 0$ or $\partial p_b/\partial x\to 0$. The first case means absence of a shear–thinning effect, as suggested by the second Equation (14), while the second case means no basic shear stress and, thus, no influence of the shear–thinning effect on the dynamics. In both cases, we expect that the Ellis fluid behaves as a Newtonian fluid. In fact, for $\widetilde{\mathrm{El}}\to 0$, we obtain critical values that coincide with those reported in the literature for Newtonian fluids [1], namely ${\mathrm{R}}_c=27.1$ and $k_c=2.33$.

Figure 5 and Figure 6 show similar findings, even if ${\mathrm{R}}_c$ and $k_c$ are presented as functions of n for different values of $\widetilde{\mathrm{El}}$. Again, the conclusion that, for $\widetilde{\mathrm{El}}\to 0$, the critical data do not depend on n is confirmed. One may also note that, by increasing the value of $\widetilde{\mathrm{El}}$, the value of ${\mathrm{R}}_c$ decreases with n. For a nonzero value of $\widetilde{\mathrm{El}}$, the critical value of $\mathrm{R}$ depends on n and this dependence becoming weaker and weaker as $n\to 1$. Figure 6 shows that, for nonzero values of $\widetilde{\mathrm{El}}$, the critical wavenumber tends to $k_c=2.33$ when $n\to 1$. In the limiting case $n\to 0$, there is a significant decrease of both ${\mathrm{R}}_c$ and $k_c$.

Figure 7 displays the isolines of both pressure and temperature, showing the cellular flow pattern for $\widetilde{\mathrm{El}}=10$ and $n=0.5$. The convective cells are evidently asymmetric with respect to the horizontal midplane. This is due to the different type of conditions imposed at the horizontal boundaries.

4. Conclusions

The linear convective instability of shear–thinning fluids saturating a porous channel has been investigated. The lower channel wall has been considered as subject to a uniform heat flux, while the upper wall is maintained at a constant temperature. The rheology of the shear–thinning fluid is described by employing the Ellis model. More precisely, a modified version of Darcy’s law has been assumed to represent the shear–thinning behaviour of the fluid. A study of small–amplitude disturbances has been proposed in order to determine the threshold for the onset of instability. The resulting system of differential equations, defining the stability eigenvalue problem, has been solved numerically. The governing parameters of this analysis are the modified Ellis number, $\widetilde{\mathrm{El}}$, the power–law index n, the Darcy–Rayleigh number $\mathrm{R}$, the wavenumber k and the inclination angle $\varphi$ of the wave vector to the basic flow direction. The main results can be summarized as follows:

• The transverse rolls turned out to be the most unstable modes.

• The neutral stability curves display, qualitatively, the same shape when changing the fluid flow parameters.

• The effect of decreasing the value of n on the threshold value for the onset of the instability is similar to the effect of increasing the value of the modified Ellis number, $\widetilde{\mathrm{El}}$: both decreasing n or increasing $\widetilde{\mathrm{El}}$ yield a destabilising effect.

• When $\widetilde{\mathrm{El}}\to 0$, the onset of the instability is not affected by n and the critical values of the governing parameter match the values reported in the literature for Newtonian fluids, namely ${\mathrm{R}}_c=27.1$ and $k_c=2.33$.

Author Contributions

Conceptualization P.V.B., M.C. and A.B.; methodology P.V.B., M.C. and A.B.; validation P.V.B., M.C. and A.B.; formal analysis P.V.B., M.C. and A.B.; investigation P.V.B., M.C. and A.B.; writing–original draft preparation P.V.B., M.C. and A.B.; writing–review and editing P.V.B., M.C. and A.B.; funding acquisition P.V.B., M.C. and A.B. All authors have read and agreed to the published version of the manuscript.

Funding

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior–Brazil (CAPES)–Grant n$${}^{\circ }$$ 88881.174085/2018–01. Financial support was also provided by Ministero dell’Istruzione, dell’Università e della Ricerca (Italy)–Grant n$${}^{\circ }$$ PRIN2017F7KZWS.

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

Figure 1. A sketch of the porous layer, of the boundary conditions and of the coordinate axes.

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

Figure 3. Critical value of [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] for different values of n.

Figure 4. Critical wavenumber [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] for different values of n.

Figure 5. Critical value of [Forumla omitted. See PDF.] as a function of n for different values of [Forumla omitted. See PDF.].

Figure 6. Critical wavenumber k as a function of [Forumla omitted. See PDF.].

Figure 7. Isolines of pressure (left) and temperature (right) drawn for [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.].

Appendix A. A Proof That $$\tilde{\mathit{\omega }}=\mathbf{0}$$

In order to prove that $$\tilde{\omega }=0$$, we multiply the first Equation (20) by the complex conjugate of the eigenfunction $$\tilde{f}$$, $${\tilde{f}}^{*}$$. Then, after integrating by parts and using the boundary conditions given by Equation (21), one obtains (A1)$$\begin{array}{c}\begin{array}{c}{\int }_0^1|{\tilde{f}}^{\prime }{|}^2\mathrm{d}y+k^2\tilde{n}{\int }_0^1{\left|\tilde{f}\right|}^2\mathrm{d}y-\tilde{\mathrm{R}}{\int }_0^{1}h{\tilde{f}}^{{*}^{\prime }}\mathrm{d}y=0.\end{array}\end{array}$$

From Equation (A1), by recognising that the first two terms are both real and positive, we reach the conclusion that (A2)$$\begin{array}{c}\begin{array}{c}{\int }_0^{1}h{\tilde{f}}^{{*}^{\prime }}\mathrm{d}y\in \mathbb{R},\end{array}\end{array}$$ where $$\mathbb{R}$$ denotes the set of real numbers. By taking the complex conjugate of the integral in Equation (A2), we conclude also that (A3)$$\begin{array}{c}\begin{array}{c}{\int }_0^1{\tilde{f}}^{\prime }{h}^{*}\mathrm{d}y\in \mathbb{R}.\end{array}\end{array}$$

We multiply the second Equation (20) by the complex conjugate of h, $$h^{*}$$, and we integrate by parts over the domain $$y\in (0,1)$$ to obtain (A4)$$\begin{array}{c}\begin{array}{c}{\int }_0^1|h^{\prime }{|}^2\mathrm{d}y+\left(k^2-\tilde{\mathrm{R}}-i\tilde{\omega }\right){\int }_0^1{\left|h\right|}^2\mathrm{d}y+{\int }_0^1{\tilde{f}}^{\prime }{h}^{*}\mathrm{d}y=0.\end{array}\end{array}$$

By utilising Equation (A3), we may now split Equation (A4) into its real and imaginary parts to get (A5)$$\begin{array}{c}\begin{array}{c}\tilde{\omega }{\int }_0^1{\left|h\right|}^2\mathrm{d}y=0.\end{array}\end{array}$$

Finally, by recalling that the trivial solution ($$h=0$$) must be excluded, we conclude that (A6)$$\tilde{\omega }=0.$$

Appendix B. Dominant Modes

If we account for the dependence of $$\tilde{n}$$ on $$\varphi$$, we can say that $$\tilde{n}$$ is a monotonic decreasing function of $$\varphi$$, considering that $$0\le \varphi \le \pi /2$$. The proof is simple: we evaluate $$\mathrm{d}\tilde{n}/\mathrm{d}\varphi$$ from Equation (22), (A7)$$\begin{array}{c}\frac{\mathrm{d}\tilde{n}}{\mathrm{d}\varphi }=-{\displaystyle \frac{\widetilde{\mathrm{El}}(1-n)sin\left(2\varphi \right)}{n(\widetilde{\mathrm{El}}+1)}}.\end{array}$$

Since $$n\le 1$$ and $$sin\left(2\varphi \right)\ge 0$$ within the interval $$0\le \varphi \le \pi /2$$, one can easily conclude that $$\mathrm{d}\tilde{n}/\mathrm{d}\varphi \le 0$$, i.e., that $$\tilde{n}$$ is a monotonic decreasing function of $$\varphi$$.

Thus, by looking at Figure A1, we can observe that $${\mathrm{R}}_c$$ is a monotonic decreasing function of $$\tilde{n}$$. Hence, we can infer that $${\mathrm{R}}_c$$ is a monotonic increasing function of $$\varphi$$. This conclusion implies that the most unstable modes are transverse $$(\varphi =0)$$.

Figure A1. Critical value of [Forumla omitted. See PDF.] as a function of [Forumla omitted. See PDF.] for different values of [Forumla omitted. See PDF.].

Appendix C. Numerical Method

The numerical method used to solve the differential eigenvalue problem is the shooting method. The basis of this method is the reformulation of the original eigenvalue problem (20), subject to the boundary conditions (21), as an initial value problem by using additional initial conditions, namely (A8)$$\begin{array}{c}\begin{array}{c}{\tilde{f}}^{″}-\tilde{n}{k}^2\tilde{f}-\tilde{\mathrm{R}}{h}^{\prime }=0,\\ h^{″}-\left(k^2-\tilde{\mathrm{R}}\right)h-{\tilde{f}}^{\prime }=0,\\ {\tilde{f}}^{\prime }\left(0\right)=\tilde{\mathrm{R}}h\left(0\right),\tilde{f}\left(0\right)=1,h^{\prime }\left(0\right)=0,h\left(0\right)=\xi .\end{array}\end{array}$$

We mention that the initial condition $$\tilde{f}\left(0\right)=1$$ can be imposed because the governing differential equations in (A8) are homogeneous. On the other hand, $$\xi$$ is an unknown real parameter to be found. The system of Equation (A8) is solved numerically by means of the Runge–Kutta method. The solution for the eigenfunctions $$\tilde{f}$$ and h depends on the four governing parameters, k, $$\tilde{n}$$, $$\tilde{\mathrm{R}}$$ and $$\xi$$. Such parameters are to be determined through a root–finding algorithm in order to satisfy the target conditions (A9)$$\begin{array}{c}\begin{array}{c}{\tilde{f}}^{\prime }\left(1\right)=0,h\left(1\right)=0.\end{array}\end{array}$$

A comparison between our results and those reported in Chapter 6 of the book by Nield and Bejan [1] reveals, in the special case $$\mathrm{El}\to 0$$, an excellent agreement, with $$k_c=2.33$$ and $${\mathrm{R}}_c=27.1$$.

References

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

2. Straughan, B. Stability and Wave Motion in Porous Media; Springer: New York, NY, USA, 2008.

3. Barletta, A. Routes to Absolute Instability in Porous Media; Springer: New York, NY, USA, 2019.

4. Kim, M.C.; Lee, S.B.; Kim, S.; Chung, B.J. Thermal instability of viscoelastic fluids in porous media. Int. J. Heat Mass Transf.; 2003; 46, pp. 5065-5072. [DOI: https://dx.doi.org/10.1016/S0017-9310(03)00363-6]

5. 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]

6. Pallavi, G.; Hemanthkumar, C.; Shivakumara, I.S.; Rushikumar, B. Oscillatory Darcy–Bénard–Poiseuille mixed convection in an Oldroyd–B fluid–saturated porous layer. Advances in Fluid Dynamics; Springer: New York, NY, USA, 2021; pp. 827-837.

7. Rees, D.A.S. Darcy–Bénard–Bingham convection. Phys. Fluids; 2020; 32, 084107. [DOI: https://dx.doi.org/10.1063/5.0018775]

8. Barletta, A.; Nield, D. Linear instability of the horizontal throughflow in a plane porous layer saturated by a power–law fluid. Phys. Fluids; 2011; 23, 013102. [DOI: https://dx.doi.org/10.1063/1.3532805]

9. Alves, L.S.d.B.; Barletta, A. Convective instability of the Darcy–Bénard problem with through flow in a porous layer saturated by a power–law fluid. Int. J. Heat Mass Transf.; 2013; 62, pp. 495-506. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2013.02.058]

10. Celli, M.; Barletta, A.; Longo, S.; Chiapponi, L.; Ciriello, V.; Di Federico, V.; Valiani, A. Thermal instability of a power–law fluid flowing in a horizontal porous layer with an open boundary: A two–dimensional analysis. Transp. Porous Media; 2017; 118, pp. 1-23. [DOI: https://dx.doi.org/10.1007/s11242-017-0863-7]

11. Petrolo, D.; Chiapponi, L.; Longo, S.; Celli, M.; Barletta, A.; Di Federico, V. Onset of Darcy–Bénard convection under throughflow of a shear–thinning fluid. J. Fluid Mech.; 2020; 889, R2. [DOI: https://dx.doi.org/10.1017/jfm.2020.84]

12. Brandão, P.V.; Ouarzazi, M.N. Darcy–Carreau model and nonlinear natural convection for pseudoplastic and dilatant fluids in porous media. Transp. Porous Media; 2021; 136, pp. 521-539. [DOI: https://dx.doi.org/10.1007/s11242-020-01523-9]

13. Brandão, P.V.; Ouarzazi, M.N.; Hirata, S.C.; Barletta, A. Darcy–Carreau–Yasuda rheological model and onset of inelastic non–Newtonian mixed convection in porous media. Phys. Fluids; 2021; 33, 044111. [DOI: https://dx.doi.org/10.1063/5.0048143]

14. Celli, M.; Barletta, A.; Brandão, P.V. Rayleigh–Bénard instability of an Ellis fluid saturating a porous medium. Transp. Porous Media; 2021; 138, pp. 679-692. [DOI: https://dx.doi.org/10.1007/s11242-021-01640-z]

15. Sochi, T. Non–Newtonian flow in porous media. Polymer; 2010; 51, pp. 5007-5023. [DOI: https://dx.doi.org/10.1016/j.polymer.2010.07.047]

16. Sadowski, T.J.; Bird, R.B. Non–Newtonian flow through porous Media. I. Theoretical. Trans. Soc. Rheol.; 1965; 9, pp. 243-250. [DOI: https://dx.doi.org/10.1122/1.549000]

17. Sadowski, T.J. Non–Newtonian flow through porous media. II. Experimental. Trans. Soc. Rheol.; 1965; 9, pp. 251-271. [DOI: https://dx.doi.org/10.1122/1.549023]

18. Alves, L.D.B.; Hirata, S.D.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]