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

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


Introduction
The emergence of convection cell patterns in a fluid saturated porous medium 13 heated from below is a cornerstone topic in the research on convection heat transfer 14 in porous solids. There is an impressively wide literature on this topic, often called 15 Darcy-Bénard instability, which is surveyed in Chapter 6 of Nield and Bejan [1]. Other 16 surveys relative to this area of research are provided by Straughan [2], Barletta [3]. Except  The Darcy-Bénard instability for saturating fluids with a non-Newtonian rheology 21 has been investigated with reference to viscoelastic fluids [4][5][6], to fluids with yield-stress 22 [7] and to purely viscous fluids [8][9][10][11][12][13][14]. 23 The usual setup envisaged in the studies on the onset of Darcy-Bénard instability 24 in a saturated porous layer is one where the plane parallel boundaries are considered as 25 isothermal with different temperatures. In the present paper, the objective is to extend 26 the analysis carried out by Celli et al. [14]. The temperature boundary conditions are 27 modified, with respect to the study developed by Celli et al. [14], by assuming a uniform 28 heat flux on the lower boundary of the layer. Such a modification is accompanied by a 29 numerical solution of the stability eigenvalue problem, whereas an analytical dispersion 30 relation is available with isothermal boundary conditions [14]. The shear-thinning 31 rheology of the fluid is formulated by adopting the Ellis model extension of Darcy's law.

32
Results will be discussed for the instability onset in terms of the neutral stability curves

38
A porous layer of height H saturated by a non-Newtonian fluid is considered. The 39 boundary walls are considered to be impermeable. The channel is heated from below by 40 a constant heat-flux q 0 , while the upper boundary is kept at a constant temperature T 0 .

41
The Oberbeck-Boussinesq approximation is employed and local thermal equilibrium 42 between the fluid and the solid matrix is assumed. A basic uniform throughflow is 43 imposed along the horizontal x direction. 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 η, according to the Ellis model, is given by [15,16] where ζ is a positive parameter such that ζ = 1/n, where n is the power-law index.
The apparent viscosity at zero shear-stress is represented by η 0 , while τ 1/2 represents the value of τ at which η = η 0 /2. By employing the power-law index, Eq.(2) can be rewritten as The apparent viscosity predicted by the Ellis model, for some limiting cases, is given by the following list: On the other hand, in the limiting case of strong shear stresses, τ τ 1/2 , Eq. (2) reduces to

46
For a Newtonian fluid saturating a porous medium, Darcy's law can be written as where u is the seepage velocity vector of components (u, v, w), K is the permeability of the porous medium and f d is the drag force which contains the pressure gradient and the buoyancy term based on the Oberbeck-Boussinesq approximation, namely In Eq. (6), p is the local difference between the pressure and the hydrostatic pressure, ρ 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. Darcy's law can be generalised for non-Newtonian fluids by replacing the fluid viscosity with an effective viscosity [16,17]. Hence, Eq. (5) becomes where η e f f is the effective viscosity given by 1 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 ϕ [16] Thus, the modified Darcy's law for a porous medium saturated by an Ellis fluid can be rewritten in a more convenient way as where A is a coefficient determined by the properties of the porous medium and of the fluid, It is worth noting that, near the instability threshold, the drag forces have, usually, a 47 small magnitude. In fact, in the limit of vanishing drag forces, f d → 0, Eq. (10) reduces 48 to Eq. (5). The mass balance, momentum balance and energy balance equations for the present problem are given by where the bars over the symbols represent the dimensional quantities. In Eq. The following scaling allows one to express Eq. (12) in a dimensionless form where ∆T = q 0 H/χ and x is the position vector with Cartesian components (x, y, z). By substituting Eq. (13) into Eqs. (12), 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 55
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 uniform pressure gradient, where the subscript b denotes the basic state. Without any loss of generality, we can 56 assume that ∂p b /∂x is negative, so that u b assumes only positive values. In order to simplify the numerical treatment, the governing equations can be rewritten by employing a pressure-temperature formulation. By considering Eq. (14a) and by applying the divergence operator to Eq. (14b), we can get the rid of the velocity field in the momentum balance equation. For the energy balance equation, we employ Eq. (14b) to express the velocity as function of pressure. Thus, the pressure-temperature formulation of the governing equations is given by ∂T ∂t where the boundary conditions based on the pressure-temperature formulation are 59 obtained by employing the governing equations (14). 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 where λ is the growth rate, k x and k z are the wavenumbers in the x and z direction, ω is  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 λ = 0, we ensure the neutral stability condition. Hence, the differential dispersion relation for the disturbances is given bỹ subject to the boundary conditions where the following quantities have been defined: Here, φ the angle between the wave vector of the normal mode and the x axis. In

66
Appendix A, a proof that the eigenvalue problem is to be solved withω = 0 is provided.

67
In Appendix B, it is shown that the most unstable modes are the transverse rolls, whose where R denotes the set of real numbers. By taking the complex conjugate of the integral in Eq. (A2), we conclude also that 1 0f h * dy ∈ R.
We multiply Eq. (20b) by the complex conjugate of h, h * , and we integrate by parts over the domain y ∈ (0, 1) to obtain 1 0 |h | 2 dy + k 2 −R − iω Finally, by recalling that the trivial solution (h = 0) must be excluded, we conclude that If we account for the dependence ofñ on φ, we can say thatñ is a monotonic decreasing function of φ, considering that 0 ≤ φ ≤ π/2. The proof is simple: we evaluate dñ/dφ from Eq. (22), Since n ≤ 1 and sin(2φ) ≥ 0 within the interval 0 ≤ φ ≤ π/2, one can easily conclude 144 that dñ/dφ ≤ 0, i.e. thatñ is a monotonic decreasing function of φ. 145 Thus, by looking at Fig. A1, we can observe that R c is a monotonic decreasing 146 function ofñ. Hence, we can infer that R c is a monotonic increasing function of φ. This 147 conclusion implies that the most unstable modes are transverse (φ = 0).

148
Appendix C. Numerical method 149 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, namelỹ We mention that the initial conditionf (0) = 1 can be imposed because the governing differential equations in (A8) are homogeneous. On the other hand, ξ is an unknown real parameter to be found. The system of equations (A8) is solved numerically by means of the Runge-Kutta method. The solution for the eigenfunctionsf and h depends on the four governing parameters, k,ñ,R and ξ. Such parameters are to be determined through a root-finding algorithm in order to satisfy the target conditions f (1) = 0, h(1) = 0. (A9) A comparison between our results and those reported in Chapter 6 of the book by Nield 150 and Bejan [1] reveals, in the special case El → 0, an excellent agreement, with k c = 2.33 151 and R c = 27.1.