Thermal Convection of an Ellis Fluid Saturating a Porous Layer with Constant Heat Flux Boundary Conditions

Abstract

The present work analyzes the thermal instability of mixed convection in a horizontal porous channel that is saturated by a shear-thinning fluid following Ellis’ rheology. The fluid layer is heated from below by a constant heat flux and cooled from above by the same heat flux. The instability of such a system is investigated by means of a small-disturbances analysis and the resulting eigenvalue problem is solved numerically by means of a shooting method. It is demonstrated that the most unstable modes on the instability threshold are those with infinite wavelength and an analytical expression for such conditions is derived from an asymptotic analysis. Results show that the non-Newtonian character of the fluid has a destabilizing role.

1. Introduction

The study of the Rayleigh–Bénard instability for a fluid-saturated porous layer is a topic of utmost importance for the analysis of convection heat transfer in porous materials. A comprehensive review of this research topic is presented in Nield and Bejan [1]. The literature on this subject is mainly focused on Newtonian fluids. When a non-Newtonian shear-thinning fluid is considered, the Ostwald–De Waele model is frequently employed for the definition of the fluid rheology [2]. Such a model displays a peculiar feature when the shear stress applied to the fluid tends to zero. For this limiting case, the apparent viscosity of the fluid tends to infinity when dealing with pseudoplastic fluids, and to zero for dilatant fluids. This is not always an issue since usually one focuses on non-vanishing mass flow rates of shear-thinning fluids. On the other hand, when dealing with transition to instability, it can be a major problem, as a fluid at rest is a possible basic state for the onset of instability.

In order to fix the singular behavior of the Ostwald–De Waele law, one can employ different rheological models. Nield [3,4] suggested a modified drag force term for Darcy’s law. According to this suggested model, the singularity problem is solved, but only for dilatant fluids. Instead, for pseudoplastic fluids, the problem persists since the apparent viscosity still goes to infinity in the limit of the vanishing shear rate. More recently, Brandao and Ouarzazi [5] and Brandao et al. [6] proposed a modified Darcy’s law based on the Carreau and the Carreau–Yasuda models, respectively. In both cases, the apparent viscosity is reduced to the Newtonian one when the shear stress is negligible. On the other hand, Celli et al. [7] and Brandao et al. [8] proposed the Ellis model as a suitable framework for the description of the rheology of shear-thinning fluids. For the Ellis model, the apparent viscosity is also reduced to the Newtonian one in the absence of shear rate. While the Carreau and Carreau–Yasuda extensions for porous media are recent developments, the Ellis model has been used for a long time in the study of non-Newtonian flow in porous media [9] and is supported by experimental validations.

The pioneering studies on thermal instability in porous media involve Dirichlet boundary conditions for the temperature field, i.e., prescribed temperature at the walls [10,11,12]. Recently, Celli et al. [7] investigated the effect of the shear-thinning characteristics on the thermal instability of an Ellis fluid subject to an unstable temperature gradient due to prescribed temperatures at the walls. It is well-known from the literature that relaxing the Dirichlet boundary conditions for the temperature field leads to a lower critical Rayleigh number for the Newtonian rheology both for clear fluids and for saturated porous media [13,14,15,16]. In addition, such a change in the boundary conditions yields the emergent convective cell to be of larger wavelength (smaller wavenumber). According to Park and Sirovich [14], a smaller value of the Rayleigh number on the instability threshold means a smaller input energy to the fluid and, to accommodate this decrease in the energy, the fluid develops convective cells of infinite wavelength.

The present study aims to investigate the threshold conditions for the onset of thermal convection when a shear-thinning fluid saturates a horizontal porous layer heated from below and cooled from above with the same uniform heat flux. The horizontal boundaries are impermeable and a horizontal throughflow is prescribed. The Ellis model is employed to describe the rheology of the fluid. The stability of a stationary basic state characterized by a uniform flow and a vertical downward temperature gradient is investigated. A linear stability analysis is performed by superposing to the basic state small-amplitude perturbations in the form of normal modes and by solving the eigenvalue problem thus obtained. The aim of this paper is to analyze how the shear-thinning character of the fluid affects the onset of thermal convection.

2. Mathematical Formulation

A horizontal porous channel saturated by a shear-thinning fluid is considered. The boundary walls, separated by a height H, are considered to be impermeable and both are subject to a constant heat flux $q_0$ (see Figure 1). The Oberbeck–Boussinesq approximation is considered and the 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 behavior of the shear-thinning fluid is described by means of the Ellis model. Such a model involves three parameters and describes the rheology of a time-independent, shear-thinning and non-yield-stress fluid [17]. According to the Ellis model, the apparent viscosity $\eta$ is given by [17,18]:

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

(1)

where a is the Ellis index, ${\eta }_0$ is the apparent viscosity at zero shear-stress and ${\tau }_{1/2}$ is the value of $\tau$ when $\eta ={\eta }_0/2$. The model is formulated with $0<a<1$, although typical values of the Ellis index are within the range $0.14<a<0.91$ [19].

2.2. Modified Darcy’s Law

The classical version of the Darcy’s law for a Newtonian fluid is given by:

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

(2)

where $\mathbf{u}$ represents the seepage velocity vector, whose components are denoted as $(u,v,w)$, K is the permeability of the porous medium and ${\mathbf{f}}_d$ is the drag force, which is given by:

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

(3)

where p is the local difference between the pressure and the hydrostatic pressure, ${\rho }_0$ is the fluid density at the reference temperature $T_0$, $\mathbf{g}$ is the gravity acceleration vector and $\beta$ is the thermal expansion coefficient of the fluid.

A generalized form of Darcy’s law is represented by:

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

(4)

where ${\eta }_{eff}$ is the effective viscosity given by [18]:

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

(5)

with $r_h$ being the mean hydraulic radius, which is directly proportional to the square root of the ratio between the permeability K and the porosity $\phi$ of the solid medium.

The modified Darcy’s law for a porous medium saturated by a shear-thinning fluid of Ellis type can thus be rewritten as:

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

(6)

where A is a coefficient that depends on the properties of the fluid and the porous medium.

2.3. Governing Equations

The governing equations for mass, momentum and energy describing the present problem are thus given by:

$$\nabla ·\mathbf{u}=0,$$

(7)

$$\frac{{\eta }_0}{K}\mathbf{u}=\left(1+A|{\mathbf{f}}_d{|}^{\frac{1-a}{a}}\right){\mathbf{f}}_d,$$

(8)

$${\mathbf{f}}_d=-\nabla p-{\rho }_0\mathbf{g}\beta (T-T_0),$$

(9)

$$\sigma \frac{\partial T}{\partial t}+\mathbf{u}·\nabla T=\alpha {\nabla }^{2}T,$$

(10)

$$y=0,H:v=0,-\chi \frac{\partial T}{\partial y}=q_0.$$

(11)

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

The following dimensionless quantities, denoted with an overline, are introduced:

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

(12)

where $\Delta T=q_{0}H/\chi$ and $\mathbf{x}=(x,y,z)$ is the position vector. Furthermore, we define $\mathrm{El}$ as the Darcy–Ellis number and $\mathrm{R}$ as the Darcy–Rayleigh number, given by:

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

(13)

By substituting Equation (12) into Equations (7)–(11) and by utilizing Equation (13), the governing equations can be written in the following dimensionless form:

$$\nabla ·\mathbf{u}=0,$$

(14)

$$\mathbf{u}=\left(1+\mathrm{El}|{\mathbf{f}}_d{|}^{\frac{1-a}{a}}\right){\mathbf{f}}_d,$$

(15)

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

(16)

$$y=0,1:v=0,-\frac{\partial T}{\partial y}=1,$$

(17)

where:

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

(18)

In Equations (14)–(18), as well as in the forthcoming analysis, the overline symbol for the dimensionless quantities has been omitted for the sake of brevity.

We note that the limit of Newtonian rheology is recovered when $\mathrm{El}\to 0$. On the other hand, a very large value of $\mathrm{El}$ drives the Ellis model to the power–law model for pseudoplastic fluids. In the latter case, the Ellis index a can be considered as equivalent to the power–law index.

2.4. Basic State

The basic stationary flow is generated by a prescribed constant pressure drop along the x axis, $\partial p_b/\partial x$, and by a prescribed heat flux at the walls. This basic solution is given by:

$$\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-a}{a}}\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}$$

(19)

with b standing for basic state. Without any loss of generality, $\partial p_b/\partial x$ can be assumed to be negative to obtain a positive basic velocity $u_b$.

The pressure–temperature formulation of the dimensionless governing Equations (14)–(17) is given by:

$$\nabla ·\left[\left(1+\mathrm{El}|{\mathbf{f}}_d{|}^{\frac{1-a}{a}}\right){\mathbf{f}}_d\right]=0,$$

(20)

$$\frac{\partial T}{\partial t}+\left[\left(1+\mathrm{El}|{\mathbf{f}}_d{|}^{\frac{1-a}{a}}\right){\mathbf{f}}_d\right]·\nabla T={\nabla }^{2}T,$$

(21)

$${\mathbf{f}}_d=-\nabla p+\mathrm{R}T{\mathbf{e}}_y,$$

(22)

$$y=0,1:\frac{\partial p}{\partial y}=\mathrm{R}T,-\frac{\partial T}{\partial y}=1.$$

(23)

2.5. Linear Stability Analysis

In order to perform the linear stability analysis, the original fields are decomposed in two parts: one relative to the basic solution and the other one relative to infinitesimal disturbances. Neutral stability conditions, i.e., zero growth rate disturbances, are sought. The linear dynamics of the infinitesimal disturbances is investigated by solving the linearized governing equations. Such disturbances are assumed to behave as plane waves, i.e., Fourier modes, namely:

$$\begin{array}{c}\begin{array}{c}p(x,y,z,t)=p_b(x,y)+\epsilon f\left(y\right)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^{i\left(kx+k_{z}z-\omega t\right)},\end{array}\end{array}$$

(24)

where $\epsilon$ is a positive small parameter that accounts for the amplitude of the disturbances, f and h are the eigenfunctions of the problem, $k_x$ and $k_z$ are the wavenumbers in the streamwise and spanwise directions, respectively, and $\omega$ is the angular frequency.

After substituting Equation (24) into Equations (20)–(23), the ordinary linear differential equations governing the disturbance dynamics can be written as:

$${\tilde{f}}^{″}-\tilde{a}{k}^2\tilde{f}-\tilde{\mathrm{R}}{h}^{\prime }=0,$$

(25)

$$h^{″}-(k^2-\tilde{\mathrm{R}}-i\tilde{\omega })h-{\tilde{f}}^{\prime }=0,$$

(26)

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

(27)

where the following quantities have been defined:

$$\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-a}{a}},\\ k_x=kcos\left(\varphi \right),k_z=ksin\left(\varphi \right),\tilde{\omega }=\omega -k\mathrm{Pe},\\ \tilde{a}={\displaystyle \frac{\widetilde{\mathrm{El}}+a(\widetilde{\mathrm{El}}+2)+\widetilde{\mathrm{El}}(1-a)cos\left(2\varphi \right)}{2a(\widetilde{\mathrm{El}}+1)}}.\end{array}\end{array}$$

(28)

In Equation (28), $\mathrm{Pe}$ is the Pèclet number and the wavenumbers $k_x$ and $k_z$ are written in terms of a single wavenumber k and the inclination angle $\varphi$, in order to account for both longitudinal and transverse rolls. It is possible to verify that $\tilde{\omega }=0$, which means that $\omega =k\mathrm{Pe}$. In order to demonstrate such a result, let us multiply Equation (25) by the complex conjugate of $\tilde{f}$ and Equation (26) by the complex conjugate of h and, then, let us integrate the resulting equations over y. By invoking integration by parts and by taking into account the boundary conditions (27), it is possible to show that all terms in the resulting equations are real, which allows us to conclude that $i\tilde{\omega }=0$. In other words, all modes on the instability threshold travel with a phase velocity that is equal to the average velocity of the fluid flow, or equivalently the stability analysis is performed in the comoving reference frame. The Péclet number is not explicitly defined here, but its definition can be recovered by averaging the velocity profile given by Equations (19), namely:

$$\mathrm{Pe}={\int }_0^1{u}_b\mathrm{d}y=\left|\frac{\partial p_b}{\partial x}\right|\left(1+\mathrm{El}{\left|\frac{\partial p_b}{\partial x}\right|}^{\frac{1-a}{a}}\right).$$

(29)

3. Asymptotic Analysis for Vanishing Wavenumber

As already mentioned, disturbance modes of a vanishing wavenumber tend to be the most unstable ones on the thermal instability threshold for a fluid layer heated by a constant heat flux. For this reason, it can be useful to investigate the stability of the present problem by focusing on vanishing wavenumbers, i.e., $k\to 0$. In order to do so, an asymptotic analysis is performed here in the vicinity of $k=0$ by expanding the eigenfunctions of $\tilde{f}$ and h, as well as $\tilde{\mathrm{R}}$, in power series of k and then substituting these series into Equations (25)–(27). Since in Equations (25) and (26) there are only quadratic terms of k, the power series expansions can be performed by considering only quadratic terms, namely:

$$\tilde{f}\left(y\right)=\sum _{j=0}^{\infty }{k}^{2j}{\tilde{f}}_{2j}\left(y\right),h\left(y\right)=\sum _{j=0}^{\infty }{k}^{2j}{h}_{2j}\left(y\right),\tilde{\mathrm{R}}=\sum _{j=0}^{\infty }{k}^{2j}{\tilde{\mathrm{R}}}_{2j}.$$

(30)

After substituting Equation (30) into the original eigenvalue problem represented by Equations (25)–(27), it is possible to solve the problem for each order j. Due to the boundary conditions, $h_0$ is defined up to a constant. Such a constant is used to fix the solution scale, which allows us to write:

$$h_0=1,$$

(31)

while $f_0$ is given by:

$${\tilde{f}}_0={\tilde{\mathrm{R}}}_{0}y+{\mathcal{C}}_0,$$

(32)

where ${\mathcal{C}}_0$ is an integration constant. With reference to second-order terms, it is possible to obtain the solutions for $h_2$ and ${\tilde{f}}_2$:

$$h_2=\frac{1}{24}\tilde{a}{y}^2\left[{\tilde{\mathrm{R}}}_0(y^2-2)+{\mathcal{C}}_0(4y-6)\right]+{\mathcal{C}}_1,$$

(33)

$$\begin{array}{c}{\tilde{f}}_2={\tilde{\mathrm{R}}}_{2}y+\frac{1}{360}{y}^2\left\{{\tilde{\mathrm{R}}}_0\tilde{a}\left[60+{\tilde{\mathrm{R}}}_0(3y-5)y\right]y-60{\tilde{\mathrm{R}}}_0\tilde{a}-30\left[{\tilde{\mathrm{R}}}_0(y-2)y+12\right]\right\}\hfill \\ \hfill +{\tilde{\mathrm{R}}}_{0}y{\mathcal{C}}_1+{\mathcal{C}}_2,\end{array}$$

(34)

where ${\mathcal{C}}_1$ and ${\mathcal{C}}_2$ are other integration constants. By applying the boundary conditions for $\tilde{f}$, given by:

$$y=0,1:{\tilde{f}}_2^{\prime }={\tilde{\mathrm{R}}}_2{h}_0+{\tilde{\mathrm{R}}}_0{h}_2,$$

(35)

it is possible to find the solution for ${\tilde{\mathrm{R}}}_0$:

$${\tilde{\mathrm{R}}}_0=\frac{12}{\tilde{a}},$$

(36)

which, written in terms of ${\mathrm{R}}_0$ and a by employing Equation (28), becomes:

$${\mathrm{R}}_0=\frac{24a}{\widetilde{\mathrm{El}}+(\widetilde{\mathrm{El}}+2)a+\widetilde{\mathrm{El}}(1-a)cos\left(2\varphi \right)}.$$

(37)

From Equation (37), one can infer that ${\mathrm{R}}_0$ is a monotonic increasing function of the inclination angle $\varphi$. Consequently, the neutral stability condition for $k\to 0$ holds for a smaller value of $\mathrm{R}$ for transverse rolls. The sensitivity of ${\mathrm{R}}_0$ to the angle $\varphi$ can be inferred also from the derivative of ${\mathrm{R}}_0$ with respect to $\varphi$:

$$\frac{\partial {\mathrm{R}}_0}{\partial \varphi }=\frac{48\widetilde{\mathrm{El}}(1-a)asin\left(2\varphi \right)}{{\left[\widetilde{\mathrm{El}}+(2+\widetilde{\mathrm{El}})a+\widetilde{\mathrm{El}}(1-a)cos\left(2\varphi \right)\right]}^2}.$$

(38)

By considering $\widetilde{\mathrm{El}}\ge 0$ and $0<a<1$, Equation (38) always returns positive values in the range $0\le \varphi \le \pi /2$. Such a result confirms the monotonicity of ${\mathrm{R}}_0$ as a function of $\varphi$, and the transverse modes as the preferred ones for the onset of thermal instability. The results for the most unstable modes can be simplified as follows:

$${\mathrm{R}}_0^T=\frac{12a}{\widetilde{\mathrm{El}}+a},$$

(39)

where the superscript T stands for transverse. The analytical expression for ${\mathrm{R}}_0$ considering longitudinal modes is given by:

$${\mathrm{R}}_0^L=\frac{12}{\widetilde{\mathrm{El}}+1}.$$

(40)

From now on, the present analysis will be focused on the most unstable modes, i.e., the transverse ones. Moreover, it is assumed that $\tilde{\omega }=0$.

4. Results and Discussion

The solution of Equations (25)–(27) was obtained numerically by means of a shooting method, the same as described in Celli et al. [7]. The computations were carried out within the Mathematica environment [20]. The limit $\mathrm{El}\to 0$ yielded the Newtonian case, whose critical conditions are given by $k_c=0$ and ${\mathrm{R}}_c=12$ [1,16].

Table 1. Neutral stability values of $\mathrm{R}$ for small values of k for the Newtonian limit $\mathrm{El}\to 0$.

k	$\mathbf{R}$
$0.1$	12.0114
$0.01$	12.0001
$0.001$	12.0000
0	12

shows a comparison between the present results from the numerical solution for small values of k, and the well-known results for the Newtonian case. The results for $k=0$ are those coming from the asymptotic analysis presented in the previous section.

The results are presented in terms of neutral stability conditions, that is, the parametric regime where the disturbances have zero growth rate. Such a condition is expressed by the pair $(k,\mathrm{R})$ for prescribed values of the control parameters $(\widetilde{\mathrm{El}},a)$. Figure 2 shows the neutral stability condition for different values of $\widetilde{\mathrm{El}}$ and a. One can conclude that the parameter a has a stabilizing effect on the basic state: by increasing a, higher values of $\mathrm{R}$ are obtained. Although $\mathrm{R}$ can vary noticeably from frame to frame in Figure 2, the trend of the neutral stability curves is similar for the different values of $\widetilde{\mathrm{El}}$ shown. Since a controls the shear-thinning character of the fluid, one may affirm that the non-Newtonian character of shear-thinning fluids has a destabilizing effect for the present problem.

One may note that, for a given a, by increasing the value of $\widetilde{\mathrm{El}}$ starting from $\widetilde{\mathrm{El}}=0$, the basic flow becomes unstable for smaller values of $\mathrm{R}$. In fact, such a result is confirmed by Figure 3. This trend can be physically explained by mentioning that an increment in $\widetilde{\mathrm{El}}$ implies an increment in the non-Newtonian character of the fluid flow as well, in analogy to what has been found for the parameter a.

In real-world problems, perturbations of any wavenumber may arise. For this reason, once the smallest $\mathrm{R}$ is exceeded in Figure 2, the flow is considered to be unstable, even if it is seemingly stable for values of k larger than the critical wavenumber. By looking at Figure 2, one may deduce that the critical conditions always hold true for $k\to 0$. In fact, Figure 3 illustrates the neutral stability conditions for $k=0.001$ and compares it with the asymptotic results given by Equation (39). Thus, we can assume, without any loss of generality, that Equation (39) also expresses the critical conditions:

$${\mathrm{R}}_c={\mathrm{R}}_0^T=\frac{12a}{\widetilde{\mathrm{El}}+a}.$$

(41)

Such a conclusion is of great importance since it allows one to obtain the critical conditions entirely analytically.

Table 2. Critical values of $\mathrm{R}$.

$\widetilde{\mathbf{El}}$	$\mathit{a}=0.8$	$\mathit{a}=0.6$	$\mathit{a}=0.4$	$\mathit{a}=0.2$
0	12	12	12	12
$0.1$	10.666667	10.285714	9.6	8
$0.5$	7.3846154	6.5454545	5.3333333	3.4285714
1	5.3333333	4.5	3.4285714	2
10	0.88888889	0.67924528	0.46153846	0.23529412
100	0.095238095	0.071570577	0.047808765	0.023952096

shows the critical value of $\mathrm{R}$ for the same parametric combination chosen to express the neutral stability condition in Figure 2.

Figure 4 shows a comparison between the neutral stability results for $k=0.001$ and the asymptotic solution for $k\to 0$. Such results confirm the trend observed in Figure 2 and suggest again that $k=0.001$ can be considered as an approximation of the critical conditions. When $\widetilde{\mathrm{El}}\to 0$, the results do not depend on a and coincide with the Newtonian limit. The same can be observed in Figure 2, as the results for $\widetilde{\mathrm{El}}\to 0$ do not depend on a.

It can be useful to express the critical conditions for the onset of instability in terms of the parameters $\partial p_b/\partial x$, a and $\mathrm{El}$, as follows:

$${\mathrm{R}}_c=12a{\left(\mathrm{El}{\left|\frac{\partial p_b}{\partial x}\right|}^{\frac{1-a}{a}}+a\right)}^{-1},$$

(42)

where Equations (28) and (41) have been employed.

Figure 5 shows the dependence of the instability threshold given in terms of ${\mathrm{R}}_c$ as a function of the basic pressure gradient. The Newtonian fluid limit ($\mathrm{El}\to 0$) yields ${\mathrm{R}}_c=12$ independently of the basic pressure gradient. Thus, Figure 5 shows that ${\mathrm{R}}_c$ for a pseudoplastic fluid is always smaller than for a Newtonian fluid. This means that a departure from the Newtonian behavior destabilizes the basic flow. Figure 5 reveals that, for $a=0.8$, there is a marked sensitivity to the basic pressure gradient when $\partial p_b/\partial x$ is very small, especially for large values of $\mathrm{El}$.

5. Conclusions

The onset of thermal convection of a shear-thinning fluid saturating a horizontal porous medium subject to constant heat fluxes at its impermeable boundaries was investigated. The Ellis model was considered in order to overcome the singular behavior of the power-law model for vanishing shear rates. A linear stability analysis was performed by considering normal mode disturbances. The most important results can be summarized as follows:

• There exists a suitable variable transformation that yields a compact representation of the stability eigenvalue problem;
• The critical conditions hold always for $k=0$. The threshold values can be obtained entirely analytically due to an asymptotic analysis performed for $k\to 0$;
• The non-Newtonian character of the fluid plays a destabilizing effect on the convective flow, namely an increasing value of the Ellis number yields a destabilization of the basic flow;
• For $\mathrm{El}\to 0$, the Ellis index a does not affect the stability conditions and the results coincide with those for the limit of Newtonian fluid already available in the literature ($k_c=0$ and ${\mathrm{R}}_c=12$);
• For large values of the Ellis number, the power-law behavior is recovered. This means that the critical Rayleigh number tends to zero.