Rayleigh-B\'enard instability of an Ellis fluid saturating a porous medium

Unlike the power-law model, the Ellis model describes the apparent viscosity of a shear-thinning fluid with no singularity in the limit of a vanishingly small shear stress. In particular, this model matches the Newtonian behaviour when the shear stresses are very small. The emergence of the Rayleigh-B\'enard instability is studied when a horizontal pressure gradient, yielding a basic throughflow, is prescribed in a horizontal porous layer. The threshold conditions for the linear instability of this system are obtained both analytically and numerically. In the case of a negligible flow rate, the onset of the instability occurs for the same parametric conditions reported in the literature for a Newtonian fluid saturating a porous medium. On the other hand, when high flow rates are considered, a negligibly small temperature difference imposed across the horizontal boundaries is sufficient to trigger the convective instability.


Introduction
The investigation of the threshold conditions for the onset of buoyancy-driven convection of non-Newtonian fluids is a research topic that displayed a significant development in the last decades [1,2,3,4,5,6,7]. Within this area of fluid dynamics, the shear-thinning fluids, also well-known as pseudoplastic fluids, are extremely common. Pseudoplastic fluids are important in different research areas. For instance, polymer solutions display shear-thinning behaviour. The same happens for some biological fluids like blood and a significant number of liquid foods [1,8,9].
The viscosity of pseudoplastic fluids is often described by employing the Ostwald-De Weale (powerlaw) model. The drawback of this model is in its singular behaviour for negligibly small shear stresses. In fact, for this particular case, the power-law model predicts that pseudoplastic fluids display an infinite apparent viscosity [10]. The Ellis model is employed to overcome this issue. This rheological model yields the Newtonian viscosity when the shear stresses applied to the shear-thinning fluid are extremely small [11].
The analysis presented in this paper is aimed to study the threshold conditions for the onset of buoyancydriven convection in shear-thinning fluids saturating a porous medium. Since the stresses involved at onset of thermal instability may be negligibly small, the Ellis model will be employed. More precisely, the Rayleigh-Bénard instability will be analysed when an Ellis fluid saturates a horizontal porous layer. Isothermal impermeable boundaries kept at different temperatures are envisaged providing a heating-from-below condition. In perspective, the results of this study are important as they can be suitable for an experimental validation by using, for instance, a Hele-Shaw cell system [12]. In fact, the most unstable rolls for shear-thinning fluids were predicted to be transverse [13,14], by employing a power-law model. The height of the horizontal porous layer is H and the boundaries of the layer are impermeable and isothermal such that a heating-from-below configuration is present. The lower boundary is held at temperature T 0 +∆T (with ∆T > 0), while the upper boundary is held at temperature T 0 , as displayed in Fig. 1. A basic throughflow is imposed by prescribing a horizontal pressure gradient.

The Ellis model
The rheological Ellis model defines the apparent viscosity η of the non-Newtonian shear-thinning fluid as reported in [15], namely where n is a positive parameter such that 0 < n < 1, τ 0 represents the value of τ at which the apparent viscosity drops by half its reference value η 0 , τ is the scalar quantity Here τ is the shear stress tensor and τ : τ = τ ij τ ij , where the Einstein notation for the sum over repeated indices is implied. The behaviour of the viscosity ratio η/η 0 versus the shear stress ratio τ /τ 0 for different values of n is reported in Fig. 2 together with the behaviour of η/η 0 versus n for different values of τ /τ 0 . In the limiting cases of τ 0 → 0, τ 0 → ∞, n → 0, n → 1, Ellis model reduces to the Newtonian model as shown in Table 1.
On the other hand, when the fluid undergoes intense shear stresses, τ τ 0 , Eq. (1) simplifies to

Ellis model and power-law model
The power-law fluid model prescribes that the apparent viscosity of the fluid be the following function of the shear stress: where χ is the consistency factor and n is the power-law index. The limiting case described in Eq.

Modified Darcy's law for an Ellis fluid
The momentum balance equation for a Newtonian fluid saturating a porous medium is Darcy's law, namely where u is the filtration velocity vector of components (u, v), K is the permeability of the porous medium and f d is the drag force defined as follows: In Eq. (6) the Oberbeck-Boussinesq approximation is invoked, p is the pressure head, ρ 0 is the fluid density evaluated at the reference temperature T 0 , g is the gravity acceleration vector and β is the thermal expansion coefficient of the fluid. A modified Darcy's law that describes a porous medium saturated by an Ellis fluid has been proposed by [16], as well as by [17], namely where A is a fluid property [(Pa/m) 1−1/n ]. In the limiting case of A|f d | 1/n−1 1, that is when negligible drag forces are acting on the fluid, Eq. (7) matches Darcy's law (5). It is worth noting that at the onset of natural convection the intensity of the drag forces may be negligibly small.

Governing equations
The governing equations describing the problem here presented are where the bars over the quantities identify dimensional fields, coordinates and time, σ is the ratio between the average volumetric heat capacity of the porous medium and the volumetric heat capacity of the fluid, and α is the average thermal diffusivity of the saturated porous medium. The drag force f d is given by Eq. (6). The following scaling allows us to express Eq. (8) in a dimensionless formulation: where x is the Cartesian position vector of components (x, y, z). By substituting Eq. (9) into Eqs. (8) one may write where The parameter El is the Darcy-Ellis number and the parameter R is the Darcy-Rayleigh number. They are defined as follows:

Basic state
The stationary solution of Eqs. (10) employed for the stability analysis is composed by a fully developed basic flow along the horizontal direction and a purely vertical constant temperature gradient. The horizontal flow is assumed to be generated by a prescribed pressure gradient, which is independent of the x and z coordinates, such that where the subscript b denote the basic state fields. It is not restrictive to assume that ∂p b /∂x 0 so that u b 0. By taking the average value of the velocity profile, one obtains the definition of the Péclet number, namely For El → 0 with |∂p b /∂x| = 0 one may simplify Eqs. (10) and (13) to obtain the basic state employed by the Prats problem [18]. For El → 0 with |∂p b /∂x| = 0 Eqs. (10) and (13)

Linear stability analysis
The system (15) is perturbed by defining the pressure and temperature fields as composed by a basic state plus small-amplitude disturbances expressed in terms of normal modes, namely Here, f and h are, in general, complex functions, λ is the growth rate, k = (k x , 0, k z ) is the wave vector, ω is the angular frequency. By assuming that the disturbance amplitude is small, ε 1, we perform a linear stability analysis where we consider only terms O(ε). The aim of the forthcoming investigation is finding the threshold for the onset of thermal convection. This threshold is obtained when the neutrally stable modes are considered. These modes are characterised by null growth rate. Thus, from now on, λ is set equal to zero. By substituting Eq. (16) into Eqs. (15), and by employing one obtainsf where φ is the inclination angle between the wave vector and the x-axis. For φ = 0 the wave vector is parallel to the x-axis so that the rolls axes are perpendicular to the basic flow (transverse rolls). For φ = π/2 the wave vector is parallel to the z-axis. In this case, the rolls axes are parallel to the basic flow (longitudinal rolls).
In Appendix A we prove analytically thatω = 0 and, hence, we conclude that the eigenvalue problem (18) features real eigenfunctions and eigenvalues. It is worth noting that the Péclet number is not present, at least explicitly, in Eqs. (18). The definition of the rescaled angular frequency is a classical practice [21,22] for this kind of problems that follows the Prats choice [18] of performing the stability analysis in the comoving reference frame. On account of Eqs. (14) and (17), one may obtain Pe as a function of El,Ẽl and n, namely

Results
The eigenvalue problem (18) is solved both numerically and analytically. The numerical procedure, reported in Appendix B, is employed for comparison with the results obtained analytically. We assume thatf and h are trigonometric functions satisfying the boundary conditions in Eq. (18c), namelỹ where is a positive integer, and B is the constant. By employing Eqs. (17)- (20), one obtains the dispersion relation R = k 2 + π 2 2 (Ẽl + 1) + 2 π 2 2 n k 2 + π 2 2 k 2 [Ẽl (1 − n) cos(2φ) +Ẽl (n + 1) + 2 n] . The most relevant parametric configuration for the stability analysis is the one characterised by the lowest values of R. It is worth noting that, in order to minimise the value of R, the integer and positive parameter must be minimum, i.e. = 1. Moreover, by recalling that 0 φ π/2, the minimum values of R are obtained for transverse rolls, φ = 0, since this angle minimises the contribution of the second term in the right-hand side of Eq. (21). Thus, at the onset of instability, Eq. (21) can be simplified to R = k 2 + π 2 k 2 (Ẽl + n) + n π 2 (Ẽl + 1) Equation (22) allows one to draw the neutral stability curves presented in Fig. 3. This figure is obtained for the sample n = 0.2 and different values ofẼl. The absolute minimum of each neutral stability curve defines the parametric threshold for the onset of convective instability. The term "critical values" is employed to denote these threshold values of the governing parameters R and k. In order to obtain the critical values, we calculate through Eq. (22) the derivative of R with respect to k. Hence, the critical values are given by R c = π 2 1 +Ẽl n (1 +Ẽl) n +Ẽl The values of R c and k c given by Eq. (23) are reported versus n in Fig. 4, for different values ofẼl. These critical values are shown to be monotonic decreasing functions of the parameterẼl, while they are monotonic increasing functions of n. The four limiting casesẼl → 0,Ẽl → ∞, n → 0, and n → 1 deserve some particular attention.
These results coincide with those reported in Barletta & Nield [23] for the same limiting case. For a finite non vanishing value ofẼl, in the limiting case n → 0 Eq. (23), simplifies to For a finite non vanishing value ofẼl, in the limiting case n → 1 Eq. (23), simplifies to Equations 24 and 26 point out that the two limitsẼl → 0 and n → 0 do not commute.  16) and (20), only one case has been reported. Figure 5 refers to transverse rolls, φ = 0, and thus it is plotted on the plane (x, y).

Conclusions
The onset of convective instability inside a horizontal porous layer saturated by a non-Newtonian fluid has been investigated. The fluid is shear-thinning and its apparent viscosity is defined by the Ellis model. The layer is heated from below and a basic horizontal pressure gradient is assumed. A linear stability analysis has been performed by means of the normal mode method. The governing parameters are the Darcy-Rayleigh number, R, the modified Darcy-Ellis number,Ẽl, and the Ellis power-law index, n. The modified Darcy-Ellis number is a function of the Péclet number associated with the basic flow rate, of the Ellis number and of the Ellis power-law index. The main conclusions drawn from the stability analysis are the following: • The critical values of the governing parameters can be expressed analytically as functions of n andẼl.
• The most unstable rolls are transverse, having their axes perpendicular to the direction of the basic throughflow.
• The angular frequency of the transverse rolls is equal to the product between the wavenumber and the Péclet number. Such rolls are non-travelling in the reference frame comoving with the basic throughflow.
• ForẼl → 0, the critical value of the Darcy-Rayleigh number tends to 4π 2 while the wavenumber approaches π. This limiting case identifies those configurations where the basic pressure gradient is absent and/or the fluid is Newtonian. The critical values of the governing parameters match those found in the literature for either the Prats problem or the Horton-Rogers-Lapwood problem. • ForẼl → ∞, the critical value of the Darcy-Rayleigh number tends to zero and the wavenumber tends to π n 1/4 . This limiting case identifies those configurations where the basic pressure gradient is extremely intense and/or the fluid is strongly shear-thinning. In other words, for this parametric configurations, a fluid characterized by an extremely low apparent viscosity is considered and thus a negligibly small temperature gap between the horizontal boundaries is sufficient to trigger the onset of convection.
• The parametersẼl and n play different roles:Ẽl has a stabilising effect on the basic state while n has a destabilising effect.
We finally point out that our study has been based on the Ellis model for the fluid rheology in order to encompass the singular behaviour of the simpler power-law model. In particular, as pointed out in Barletta & Nield [23], the use of the power-law model leads to the prediction of an either zero or infinite critical value of the Darcy-Rayleigh number when the flow rate in the basic state is zero. On the other hand, when the basic flow rate tends to zero, the use of the Ellis model leads to a non-singular behaviour where the same critical value of the Darcy-Rayleigh number as predicted for the case of a Newtonian fluid, namely 4 π 2 , is attained.
A Proof thatω = 0 One can multiply Eq. (18a) byf * , that is the complex conjugate of the eigenfunctionf , and integrate by parts over the domain y ∈ (0, 1) to obtain From Eq. (28), one may conclude that the last integral on the left hand side is real. By taking the complex conjugate of this integral and, on integrating it by parts, one concludes that is real. This result will be invoked later on. One can now multiply Eq. (18b) by h * , that is the complex conjugate of the eigenfunction h, and integrate by parts over the domain y ∈ (0, 1) to obtain

B Numerical method
The numerical method employed to solve the stability eigenvalue problem is the shooting method. The first step consists in defining (and solving) the initial value problem obtained from Eq. (18) simplified as a consequence of the results reported in Appendix A, namelỹ Here, the conditionf (0) = 1 can be imposed because the governing equations in Eqs. (32) are homogeneous, while ξ is an unknown real parameter. The problem (32) is solved numerically by means of the Runge-Kutta method. The obtained eigenfunctionsf and h depend on four governing parameters, (k,ñ,R, ξ). The second step of the shooting method is based on the target conditions f (1) = 0, h(1) = 0.
Such conditions serve to obtain numerically, by employing a root-finding algorithm, two out of the four governing parameters (k,ñ,R, ξ). Thus, for every givenñ, one obtains the neutral stability curveR(k).
The critical values are obtained by solving the initial value problem given by Eq. (32) and the derivative with respect to k of Eq. (32). The conditions employed in the root-finding algorithm are the two conditions given by Eq. (33) together with their derivatives with respect to k. A comparison between the results obtained analytically and those obtained numerically is reported in Table 2.
The critical values of the wavenumber k and the critical values of R, both evaluated for n = 0.2, φ = 0 and different values ofẼl, are provided in this table. The subscript a refers to the data obtained analytically, while the subscript n is relative to the numerical data. The results obtained by employing these two different approaches coincide within 12 significant figures.