- Email: [email protected]
Linear instability of the vertical throughflow in a horizontal porous layer saturated by a power-law fluid
International Journal of Heat and Mass Transfer 99 (2016) 293–302
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Linear instability of the vertical throughflow in a horizontal porous layer saturated by a power-law fluid A. Barletta a,⇑, L. Storesletten b a b
Department of Industrial Engineering, Alma Mater Studiorum Università di Bologna, Viale Risorgimento 2, 40136 Bologna, Italy Department of Mathematics, University of Agder, Postboks 422, 4604 Kristiansand, Norway
a r t i c l e
i n f o
Article history: Received 24 February 2016 Received in revised form 29 March 2016 Accepted 29 March 2016 Available online 16 April 2016 Keywords: Porous medium Non-Newtonian fluid Vertical throughflow Convective instability Rayleigh number
a b s t r a c t The effects of the vertical throughflow of a non-Newtonian power-law fluid on the onset of convective instability in a horizontal porous layer are investigated. The extended Darcy’s model of momentum diffusion is employed together with the Oberbeck–Boussinesq approximation. A stationary basic solution for the vertical throughflow is determined analytically. The basic velocity and temperature fields turn out to be independent of the non-Newtonian rheology. A linear stability analysis is carried out, leading to a fourth-order eigenvalue problem. A numerical solution of the eigenvalue problem is employed to obtain the neutral stability curves and the critical Rayleigh number for the onset of instability. The governing parameters of the transition to instability are the Péclet number associated with the throughflow, and the power-law index of the fluid. These parameters influence the position of the neutral stability curve and also the critical Rayleigh number. The asymptotic cases of absolute pseudoplasticity, absolute dilatancy, and small Péclet number are discussed in detail. The latter case leads to a simple analytical solution that approximates fairly well the numerical data when the Péclet number is smaller than unity. Ó 2016 Published by Elsevier Ltd.
1. Introduction The vertical temperature gradient in Earth’s crust, if large enough, may produce fluid flows with local free convection cells. As a consequence, the porous flow circulation of groundwater can drive soluble minerals in different regions. This may lead to the so-called mineralization problem [1]. Buried nuclear debris or chemical contaminants dispersion can be driven by free convection flow of groundwater in Earth’s crust and thus yield a severe contamination of aquifers. A system involving a porous layer heated from below and saturated by a Newtonian or non-Newtonian fluid is the basic model for the study of these phenomena. Free convection cells produced in such porous layer can be combined with a forced vertical throughflow caused by pressure gradients. For instance, in geophysics, a net vertical through flow may be caused by lithostatic pressure gradients in the lower crust [1]. The modelling of convective instability in a horizontal porous layer heated from below and subject to a vertical throughflow was developed in the pioneering papers by Wooding [2], by Sutton [3], and by Homsy and Sherwood [4]. The flow in the saturated porous medium was studied by adopting a scheme based on ⇑ Corresponding author. E-mail addresses: [email protected] (A. Barletta), leiv.storesletten@ uia.no (L. Storesletten). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.03.115 0017-9310/Ó 2016 Published by Elsevier Ltd.
Darcy’s law and on the Oberbeck–Boussinesq approximation. Wooding [2] assumed a layer with infinite thickness, while both Sutton [3]and Homsy and Sherwood [4] considered a layer with a finite thickness. The stability of the vertical throughflow was studied by a linear analysis [2–4], as well as by an energy method nonlinear analysis [4]. Later on, further studies [5–14] were carried out to investigate different effects such as anisotropy, nonlinearity, internal heat generation, viscous dissipation, heterogeneity and local thermal nonequilibrium. An exhaustive survey of the literature on this topic can be found in the book by Nield and Bejan [15]. The main feature of the throughflow is that its effect is stabilising, if compared with the pure Rayleigh–Bénard instability occurring in the porous layer when the throughflow is absent. A slight departure from this behaviour was observed by Jones and Persichetti [5] for the case when the bounding planes of the porous layer are assumed to be isobaric and, thus, permeable to the unstable convective rolls. Many problems in geophysics feature non-Newtonian fluids as, for instance, those involving hydrocarbons and oil reservoirs. In these cases, the analysis of the onset of convection cells in a porous layer implies a departure from the classical Darcy’s law, due to the non-Newtonian rheology of the fluid. A typical non-Newtonian behaviour is the shear-thinning or shear-thickening nature of some fluids. This behaviour can be modelled by the Ostwald–de Waele power law [16]. Christopher and Middleman [17] formulated a
294
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
Nomenclature a A0 ; B0 b ~ b C C ez f F g H Jm K n n0 p
dimensionless wave number, Eq. (18) dimensionless integration constants, Eq. (45) dimensionless parameter, Eq. (25) dimensionless parameter, Eq. (34) dimensionless integration constant, Eq. (36) tortuosity factor, Eq. (4) unit vector in the z direction solution of the Helmholtz equation, Eq. (18) dimensionless function, Eq. (20) gravitational acceleration, with modulus it g thickness of the porous layer Bessel function of first kind and order m permeability power-law index dimensionless parameter, Eq. (56) difference between the pressure and the hydrostatic pressure Pe P´eclet number, Eq. (7) q integer number, Eq. (45) Ra Rayleigh number, Eq. (7) Rac;max maximum critical Rayleigh number, Eq. (55) s dimensionless variable, Eqs. (25) and (34) t time T temperature T0 reference temperature u velocity u; v ; w Cartesian velocity components W dimensionless amplitude of the vertical velocity disturbance, Eq. (17) W 0 ; H0 ; k0 zero order terms of the expansion (42)
modified Darcy’s law to be applied when the saturating fluid has a power-law behaviour. Several further studies on the modified Darcy’s law for power-law fluids were carried out since this groundbreaking paper, providing both experimental and theoretical results to validate this model. Thorough reviews and bibliographies on this topic can be found in the papers by Shenoy [18] and by Di Federico et al. [19]. Several studies have been recently published on flow and convection of power-law fluids in porous media [20–30]. The aim of this paper is to extend the classical analysis of the onset of convective instability in a horizontal porous layer with vertical throughflow to a non-Newtonian fluid of the power-law type. We mention that a similar analysis has been carried out for a viscoelastic fluid in a couple of papers by Shivakumara and Sureshkumar [31,32]. Our linear stability analysis is meant to determine the neutral stability curves for different power-law indexes, n, and different Péclet numbers, Pe, associated with the vertical throughflow. This analysis will be carried out numerically, even if an analytical evaluation of the neutral stability data is done for the limiting cases of absolute pseudoplasticity ðn ! 0Þ, absolute dilatancy ðn ! 1Þ, or very small Péclet numbers. The special case of absent throughflow ðPe ¼ 0Þ yields a perfect agreement with the conclusions drawn by Barletta and Nield [21].
W 1 ; H1 ; k1 first order terms of the expansion (42) x; y; z Cartesian coordinates Ym Bessel function of second kind and order m Greek symbols b thermal expansion coefficient c_ shear rate dz numerical step in the z direction DT temperature difference e dimensionless perturbation parameter, Eq. (13) g dimensionless complex parameter H dimensionless amplitude of the temperature disturbance, Eq. (17) , effective thermal diffusivity k dimensionless parameter, Eq. (15) ~k dimensionless parameter, k=n l consistency factor l effective consistency factor m dimensionless parameter, Eq. (37) q0 reference fluid density r heat capacity ratio s shear stress / porosity Superscripts, subscripts ^ perturbation of the basic flow 0 derivative with respect to z b basic solution c critical value
reference temperature and DT is a positive temperature difference. The velocity field u has Cartesian components ðu; v ; wÞ. 2.1. Governing equations For this kind of non-Newtonian fluid, the shear stress s (SI units are Pa) is a function of the shear rate c_ (SI units are s1 ) given by
s ¼ l ðc_ Þn ;
ð1Þ
where l is the consistency factor and n is the power-law index. In the case of a Newtonian fluid, we have n ¼ 1 and the consistency factor coincides with the dynamic viscosity. The SI measurement units of the consistency factor are Pa sn and, for a Newtonian fluid, the SI units of l are Pa s. Pseudoplastic fluid behaviour is obtained with n < 1, while dilatant fluids are described by n > 1. When a power-law fluid saturates a porous medium and the flow is driven by the thermal buoyancy force, modelled through
2. Mathematical model We consider a plane horizontal porous layer with thickness H, saturated by an Ostwald–de Waele (power law) fluid. Let the vertical axis z be parallel to the gravitational acceleration g, but with opposite direction (see Fig. 1). The boundary planes, z ¼ 0 and z ¼ H, are permeable and kept isothermal at T ¼ T 0 þ DT and T ¼ T 0 , respectively. Here T is the temperature field, T 0 is a
Fig. 1. Sketch of the fluid saturated porous layer and of the boundary conditions.
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
the Oberbeck–Boussinesq approximation, the classical Darcy’s law is generalised as [18,19]
$ u ¼ 0;
l
@T þ u $T ¼ r2 T; @t z ¼ 0 : w ¼ Pe; T ¼ 1;
K
jujn1 u ¼ $p q0 g bðT T 0 Þ;
ð2Þ
where l is the effective consistency factor (SI units are Pa sn m1n ), K is the permeability (SI units are m2 ) and p is the dynamic pressure (SI units are Pa), viz. the local difference between the pressure and the hydrostatic pressure. Moreover, q0 is the fluid density at the reference temperature T 0 ; b is the thermal expansion coefficient of the fluid, and g is the gravity acceleration whose modulus is denoted as g. The ratio l =K depends on the consistency factor l, on the power-law index n, on the permeability K, on the porosity /, and on the tortuosity factor C through the relationship [18,19],
l K
¼ 2Cl
3/ 50K
ðnþ1Þ=2 n 3n þ 1 : n/
ð3Þ
The tortuosity factor C was regarded by Christopher and Middleman [17] as a numerical constant independent of n and equal to 25=12, so that l =K coincides with l=K in the Newtonian limit, n ! 1. Let us assume the validity of the generalised Darcy’s law, Eq. (2). With the assumptions above, the local mass, momentum, and energy balance yield the governing equations
$ u ¼ 0;
l
n1
juj
ð4aÞ
u ¼ $p q0 g bðT T 0 Þ;
ð4bÞ
K @T r þ u $T ¼ , r2 T; @t
ð4cÞ
where t is the time, r is the ratio between the average volumetric heat capacity of the saturated porous medium and the volumetric heat capacity of the fluid, and , is the effective thermal diffusivity of the porous medium. A vertical throughflow is considered, so that the boundary conditions are given by
z¼0:
w ¼ w0 ;
T ¼ T 0 þ DT;
z¼H:
w ¼ w0 ;
T ¼ T 0;
ð5Þ
where w0 is the prescribed vertical throughflow velocity. 2.2. Dimensionless analysis We introduce a dimensionless formulation by scaling the dimensional quantities and operators as follows:
1 H H ðx; y; zÞ ! ðx; y; zÞ; u ¼ ðu; v ; wÞ ! ðu; v ; wÞ ¼ u; H , , T T0 2 ! T; H $ ! $; H r2 ! r2 : DT
, t ! t; rH2 ð6Þ
We also introduce the dimensionless parameters,
q g b DT K H n w0 H Ra ¼ 0 n ; Pe ¼ : l , ,
ð7Þ
The dimensionless parameter Ra is the non-Newtonian version of the Rayleigh number for a saturated porous medium, and Pe is the Péclet number. While Ra can only be positive, the Péclet number can either be positive (upward throughflow) or negative (downward throughflow). By taking the curl of Eq. (4b), the dimensionless formulation of the governing equation (4) and of the boundary conditions (5) is given by
$ jujn1 u ¼ Ra $ ðT ez Þ;
z¼1:
w ¼ Pe;
T ¼ 0;
295
ð8aÞ ð8bÞ ð8cÞ ð8dÞ
where ez is the unit vector in the z direction. 2.3. Basic solution A basic stationary solution of Eq. (8) is given by a uniform throughflow,
v b ¼ 0;
ub ¼ 0;
wb ¼ Pe;
ð9Þ
with a purely vertical temperature gradient dT b ðzÞ=dz. Here, the subscript b denotes the ‘‘basic solution”. The basic temperature distribution T b ðzÞ must satisfy Eqs. (8) and (9), thus yielding the differential equation 2
d Tb dT b Pe ¼ 0; dz2 dz
ð10Þ
together with the boundary conditions
T b ð0Þ ¼ 1;
T b ð1Þ ¼ 0:
ð11Þ
The solution is given by
T b ðzÞ ¼
ePe ePe z : ePe 1
ð12Þ
3. Linear stability analysis 3.1. Disturbance equations We perturb the basic solution, Eqs. (9) and (12), by defining the velocity and temperature disturbances
^; u ¼ ub þ e u
T ¼ T b þ e Tb ;
ð13Þ
where e 1 is a positive perturbation parameter. We substitute Eq. (13) into Eq. (8) and we neglect terms of order higher than e. Then, we obtain
^ @ v^ @ w ^ @u þ ¼ 0; þ @x @y @z
ð14aÞ
n
^ @ v^ @w @ Tb ; ¼k @y @z @y
ð14bÞ
n
^ @u ^ @w @ Tb ¼k ; @x @z @x
ð14cÞ
^ @ v^ @u ¼ 0; @y @x
ð14dÞ
@ Tb dT b @ Tb @ 2 Tb @ 2 Tb @ 2 Tb ^ þw ¼ 2þ 2þ 2; þ Pe @x @y @z @t @z dz
ð14eÞ
^ ¼ 0; w
z ¼ 0; 1 :
Tb ¼ 0;
ð14fÞ
where
k¼
Ra jPejn1
:
ð15Þ
We derive Eq. (14b) with respect to y, and Eq. (14c) with respect to x. Then, by summing the two resulting equations and by taking into b Þ formulation of the stability ^ T account Eq. (14a), we obtain a ðw; problem, namely
296
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
! ! ^ @2w ^ ^ @2w @2w @ 2 Tb @ 2 Tb ; n þ 2 þ 2 ¼k þ @x2 @y2 @x2 @y @z @ Tb dT b @ Tb @ 2 Tb @ 2 Tb @ 2 Tb ^ þw ¼ 2þ 2þ 2; þ Pe @x @y @z @t @z dz b ^ ¼ 0; T ¼ 0; z ¼ 0; 1 : w
ð16aÞ ð16bÞ ð16cÞ
Due to the linearity of our analysis, an arbitrary disturbance can be properly constructed by superposition of normal modes,
^ ¼ WðzÞ egt f ðx; yÞ; w
Tb ¼ HðzÞ egt f ðx; yÞ;
ð17Þ
where g is a complex parameter whose real part expresses the growth rate of the disturbance, while its imaginary part yields the angular frequency. A zero growth rate defines the condition of neutral stability, i.e., the threshold between stability and instability. We assume the validity of the principle of exchange of stabilities discussed, for the Newtonian case ðn ¼ 1Þ, by Homsy and Sherwood [4]. We mention that, for a power-law fluid saturating a porous medium, the principle of exchange of stabilities has been proved rigorously by Alves and Barletta [28] with reference to a case where the throughflow is horizontal. Thus, hereafter, the parameter g will be regarded as real, meaning that neutrally stable modes will be considered as time-independent, i.e., with g ¼ 0. In Eq. (17), f ðx; yÞ is any solution of the two-dimensional Helmholtz equation,
@2f @2f þ þ a2 f ¼ 0; @x2 @y2
ð18Þ
where a > 0 is the wave number. By substituting Eqs. (17) and (18) into Eq. (16) and by setting g ¼ 0, we obtain the eigenvalue problem for the neutrally stable modes, namely
W 00 a2 ðn W k HÞ ¼ 0; 00
0
H Pe H a H þ Pe FðzÞ W ¼ 0; z ¼ 0; 1 : W ¼ 0; H ¼ 0; 2
ð19aÞ ð19bÞ ð19cÞ
where
FðzÞ ¼
1 dT b ePe z ¼ Pe Pe dz e 1
ð20Þ
and primes are used to denote derivatives with respect to z. The fundamental symmetry of the eigenvalue problem (19) is an invariance under the transformation,
z ! 1 z;
Pe ! Pe:
ð21Þ
As a consequence, for a fixed n, changing the sign of Pe means that WðzÞ and HðzÞ are transformed into Wð1 zÞ and Hð1 zÞ, respectively. On the other hand, the eigenvalue Ra as a function of a depends on n and on the absolute value of Pe, but not on its sign. In fact, function RaðaÞ determined by solving Eq. (19) defines the neutral stability curves in the ða; RaÞ plane. Thus one may draw, without any loss of generality, the neutral stability curves only for positive values of Pe, as they depend just on jPej and n. Also the critical values ðac ; Rac Þ, namely the values of ða; RaÞ yielding the point of minimum Ra along the neutral stability curve, depend only on jPej and n. We point out that the symmetry defined by Eq. (21) had been previously revealed for Newtonian fluids ðn ¼ 1Þ [4]. 3.2. Numerical solution The eigenvalue problem controlling the onset of instability, Eq. (19), can be solved numerically. An accurate and efficient procedure is based on the combined use of a Runge–Kutta solver and of the shooting method. One has to formulate Eq. (19) as an initial value problem with an extra unknown parameter, namely
Table 1 Neutral stability values of Ra for n ¼ 0:5 and Pe ¼ 1, with either a ¼ 4 or a ¼ 6: comparison between the Runge–Kutta solver with fixed step-size, dz, and the Runge– Kutta solver with adaptive step-size algorithm. dz
a¼4
a¼6
0.00100 0.00050 0.00010 0.00005 0.00001 Adaptive
29.8029489571230 29.8030190532842 29.8030414820325 29.8030422001112 29.8030424329556 29.8030424388895
36.2054197360977 36.2054769320744 36.2054952405836 36.2054958143489 36.2054959956716 36.2054960084208
W 00 a2 ðn W k HÞ ¼ 0;
ð22aÞ
H00 Pe H0 a2 H þ Pe FðzÞ W ¼ 0; z ¼ 0 : W ¼ 0; W 0 ¼ 1; H ¼ 0;
ð22bÞ
H0 ¼ n:
ð22cÞ
In fact, n is the extra unknown parameter to be determined so that the target conditions Wð1Þ ¼ 0 and Hð1Þ ¼ 0 are satisfied. The initial condition W 0 ð0Þ ¼ 1, in Eq. (22c), can be imposed without any loss of generality as a consequence of the scale invariance of the solutions ðW; HÞ of Eq. (19). The Runge–Kutta method [33] serves to solve Eq. (22) for assigned ða; n; k; Pe; nÞ. For prescribed input parameters ða; n; PeÞ, one can determine the eigenvalue pair ðk; nÞ by imposing the target conditions Wð1Þ ¼ 0 and Hð1Þ ¼ 0. The latter operation can be accomplished by the shooting method and implemented through any root-finding algorithm, such as the Newton–Raphson technique [33]. Function NDSolve of Mathematica 10 (Ó Wolfram Research, Inc.) is employed to implement the Runge–Kutta solver, while the shooting method is applied by using function FindRoot [34]. The default option of NDSolve is based on an adaptive step-size algorithm so that the interval 0 6 z 6 1 is subdivided into steps of variable size, dz. In order to test the accuracy of the solver, the Runge–Kutta method with adaptive step-size is compared with the fixed step-size Runge–Kutta method for decreasing values of dz. In Table 1, the comparison is relative to the neutral stability values of Ra with either a ¼ 4 or a ¼ 6, and the test case is n ¼ 0:5 and Pe ¼ 1. One observes a monotonic convergence of the fixed step-size algorithm to the data obtained by the adaptive step-size algorithm as dz gradually decreases. The comparison between the data obtained with dz ¼ 105 and those with adaptive step-size yields an agreement within no fewer than eight significant figures. 4. Discussion of the results 4.1. The neutral stability curves Fig. 2 displays neutral stability curves on the parametric plane ða; RaÞ. In the cases considered, the Péclet number Pe ranges from 0:1 to 2 and the power-law index n from 0:5 to 2. The neutral stability curve, in each case, confines the region of linear stability laying below the curve from the region of instability above the curve. The neutral stability curve is affected by both parameters Pe and n. The effect of an increasing n for a given Pe is different for either smaller or larger values of Pe. As shown in Fig. 2, when Pe ¼ 0:1 and 0:2, an increasing value of n has a destabilising effect. This behaviour changes gradually when Pe ¼ 0:5 and 0:7, as the transition to instability depends non-monotonically on n. For larger Péclet numbers, such as Pe ¼ 1 and 2, the effect of an increasing n is stabilising. Thus, pseudoplastic fluids are more stable than dilatant fluids when Pe is small, while they become more unstable if Pe is larger. This behaviour is partly a consequence of the scaling with jPejn1 of the Rayleigh number Ra, implicit in the definition (15) of the eigenvalue k. However, this is just a rough
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
297
Fig. 2. Neutral stability curves in the plane ða; RaÞ with different values of Pe and n. The dashed curves are relative to the asymptotic solution for small Péclet numbers given by Eq. (52).
aspect of the actual trend, due to the sensible dependence of k on n at neutral stability. That things are more complicated than the law
collapses into the single fourth-order differential equation in the unknown W given by
Ra jPejn1 is made clear by the frame with Pe ¼ 1 in Fig. 2, where Ra displays an evident change with n. The dashed curves, drawn in Fig. 2, are relative to the asymptotic case of small Péclet numbers which will be discussed in Section 4.1.3.
W 0000 Pe W 000 a2 W 00 a2 k Pe FðzÞ W ¼ 0;
4.1.1. The asymptotic case of absolute pseudoplasticity, n ! 0 A mathematically interesting asymptotic condition is achieved when the power-law index becomes vanishingly small. When n ! 0, the system of two second-order differential equation (19)
z ¼ 0; 1 :
W ¼ 0;
W 00 ¼ 0:
ð23aÞ ð23bÞ
The trend of the neutral stability curves plotted in Fig. 2 suggests that, with n > 0; Ra tends to infinity when a ! 1. Incidentally, this is a usual feature, in the domain of small wave lengths, for the stability analysis of a Rayleigh–Bénard setup. On the other hand, if one considers Eq. (23), the behaviour for a ! 1 is described by
298
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
W 00 þ k Pe FðzÞ W ¼ 0;
ð24aÞ
z ¼ 0; 1 :
ð24bÞ
W ¼ 0:
By defining, 2
b ¼
k ; PeðePe 1Þ
s ¼ ePe z ;
ð25Þ
Eq. (24) can be rewritten as
ð26aÞ ð26bÞ
The general solution of Eq. (26a) can be expressed in terms of Bessel functions as
pffiffi pffiffi W ¼ J0 2 b s þ C Y 0 2 b s :
ð27Þ
The constants b and C are determined so that the boundary conditions (26b) are satisfied. Thus, on account of Eq. (26b), one can evaluate b and C for a given Péclet number by solving the system
J0 ð2 bÞ þ C Y0 ð2 bÞ ¼ 0; J0 2 b ePe=2 þ C Y0 2 b ePe=2 ¼ 0;
ð28Þ
Y0 ð2 bÞ J0 2 b ePe=2 J0 ð2 bÞ Y0 2 b ePe=2 ¼ 0:
ð30Þ
Ra ¼ 1:445796;
ð31Þ
which represents the smallest asymptotic value displayed for a neutral stability curve. Thus, one may imagine that the neutral stability curves of Fig. 3 become closer and closer as Pe increases, and their asymptotes tend to the value reported in Eq. (31) when Pe ! 1. One may well say that the value of Ra given by Eq. (31) is the lower bound to the critical values of Ra when n ! 0. 4.1.2. The asymptotic case of absolute dilatancy, n ! 1 Another asymptotic regime is obtained when n is extremely large. No doubt that this condition is unlikely for real fluids, but nonetheless it deserves some mathematical interest. One may take the limit n ! 1 of Eq. (19) by defining
k ¼ n ~k
or, equivalently, one can obtain b by solving the equation
pffiffiffiffiffiffi J0 2 Ra ¼ 0:
On evaluating the lowest root of Eq. (30), one obtains
2
d W dW 2 þ b W ¼ 0; þ ds2 ds s ¼ 1; ePe : W ¼ 0:
s
One may wonder whether the neutral stability curves cross the axis Ra ¼ 0 for a sufficiently high Péclet number. This is impossible physically as it would mean that thermal instability can arise even when the upper boundary is hotter than the lower boundary. Moreover, this is also impossible mathematically as, on account of Eq. (25) in the limit Pe ! 1, Eq. (29) tends to
ð29Þ
Eqs. (25) and (29) can be used to evaluate the asymptotic value of Ra, when a ! 1, for a prescribed Pe. Fig. 3 shows the neutral stability curves in the plane ða; RaÞ for the limiting case n ! 0. For a given Péclet number, this figure reveals that Ra is a monotonic decreasing function of a along the neutral stability curve. When a ! 1, an asymptote is reached where the value of Ra is that obtained by solving Eqs. (25) and (29). These asymptotes are displayed in Fig. 3 as dashed lines. Thus, one may identify the asymptotic values of Ra when a ! 1 with the critical values of Ra for the onset of instability.
ð32Þ
and by assuming ~ k Oð1Þ. Thus, Eq. (19a) yields W ¼ ~ k H and Eqs. (19b) and (19c) yield
W 00 Pe W 0 a2 W þ ~k Pe FðzÞ W ¼ 0;
ð33aÞ
z ¼ 0; 1 :
ð33bÞ
W ¼ 0:
In analogy to the procedure for the analytical solution of Eq. (24), we define
~2 ¼ b
k~ ; PeðePe 1Þ
s ¼ ePe z :
ð34Þ
Eq. (33) can be rewritten as
s2
2 2 d W ~2 s a W ¼ 0; þ b ds2 Pe2
s ¼ 1; ePe :
W ¼ 0:
ð35aÞ ð35bÞ
The general solution of Eq. (35a) can be expressed again in terms of Bessel functions,
W¼
pffiffii pffiffih ~ pffiffi ~ s ; s Jm 2 b s þ C Y m 2 b
ð36Þ
where C is an integration constant and
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 a2 m ¼ 1 þ 2: Pe
ð37Þ
The general solution (36) satisfies the boundary conditions (35b) if the following dispersion relation holds:
~ J 2b ~ ePe=2 J 2 b ~ Ym 2 b ~ ePe=2 ¼ 0: Ym 2 b m m
ð38Þ
Eqs. (37) and (38) can be employed to draw the neutral stability curves in the plane ða; ~ kÞ for different values of Pe. Fig. 4 displays these neutral stability curves. This figure reveals that the cases having Pe < 0:2 cannot be distinguished graphically, while slight differences can be seen for Pe ¼ 0:5 and 0:7 becoming more significant ~ are finite so when Pe ¼ 1 or 2. One may note that the values of k Fig. 3. Asymptotic case n ! 0: neutral stability curves in the plane ða; RaÞ with different values of Pe. The dashed lines denote the critical value of Ra attained for a ! 1 and evaluated by solving Eq. (29).
that, for n 1, one may approximate the neutral stability values of Ra by invoking Eqs. (15) and (32), namely
Ra ¼ n ~k jPejn1 :
ð39Þ
299
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
W 0 ðzÞ ¼ A0 sinðq pzÞ;
H0 ðzÞ ¼ B0 sinðq pzÞ;
q ¼ 1; 2; 3; . . . ; ð45Þ
where A0 and B0 are constants. By substituting Eq. (45) into Eq. (44), we obtain a linear system in the unknowns A0 and B0 ,
A0 q2 p2 þ a2 n a2 k0 B0 ¼ 0; 2 2 A0 B0 q p þ a2 ¼ 0:
ð46aÞ ð46bÞ
Nonzero solutions ðA0 ; B0 Þ are allowed only if the determinant of the coefficient matrix is zero, namely if
k0 ¼
q2 p2 þ a2 n q2 p2 þ a2 : a2
ð47Þ
Eq. (47) yields the neutral stability condition in the limit Pe ! 0, i.e. when the vertical throughflow is switched off. The lowest branch of neutral stability is that with q ¼ 1,
k0 ¼
p2 þ a2 n p2 þ a2 a2
;
ð48Þ
yielding a minimum of k0 for the critical values ac and k0c given by
p
ffiffiffi ; ac ¼ p 4 n ~ with Fig. 4. Asymptotic case n ! 1: neutral stability curves in the plane ða; kÞ different values of Pe. The curve for Pe ¼ 0 is drawn by employing Eq. (53).
An immediate consequence of Eq. (39) is that, at neutral stability,
lim Ra ¼ 0
ð40Þ
n!1
for every a if jPej < 1, while
lim Ra ¼ 1
ð41Þ
n!1
for every a if jPej P 1. Eqs. (40) and (41) allow us to extrapolate the trend of the neutral stability curves displayed in Fig. 2 when n increases to large values far above n ¼ 2. If 0 < Pe < 1, these curves are expected to move downward in the plane ða; RaÞ until they eventually collapse to the horizontal line Ra ¼ 0 when n ! 1. If Pe P 1, the neutral stability curves with n > 2 are expected to gradually move upward, as n increases, and eventually tend to infinity when n ! 1. In other words if n ! 1, Eqs. (40) and (41) predict instability for every value of Ra when jPej < 1, or stability for every value of Ra when jPej P 1. 4.1.3. An asymptotic solution for Pe 1 An asymptotic solution can be obtained in the limiting case of a very small Péclet number. We start from Eq. (19) and we expand the eigenfunctions W and H, as well as the parameter k, in powers of Pe, namely
WðzÞ ¼ W 0 ðzÞ þ W 1 ðzÞ Pe þ OðPe2 Þ; ð42Þ
From Eq. (20) we also obtain
1 þ OðPe2 Þ: Pe FðzÞ ¼ 1 þ Pe z 2
ð43Þ
Therefore, Eq. (19) yield to order zero,
W 000 a2 ðn W 0 k0 H0 Þ ¼ 0; 00 0
H a H0 þ W 0 ¼ 0; z ¼ 0; 1 :
W 0 ¼ 0;
ð44aÞ ð44bÞ
2
H0 ¼ 0:
Eq. (44) can be solved analytically with
pffiffiffi 2 nþ1 :
ð49Þ
On account of Eq. (15), we point out that Eq. (49) means that, when Pe ! 0, the critical value of Ra is either 0, if n > 1, or infinite, if n < 1. Physically, this means that dilatant fluids are always unstable when the vertical throughflow is switched off, while pseudoplastic fluids are always stable. Eqs. (15) and (49) imply that Newtonian fluids are a special case where Rac ¼ 4 p2 , when Pe ! 0. These results are perfectly consistent with the conclusions drawn by Barletta and Nield [21] with reference to the Rayleigh–Bénard instability of a power-law fluid saturating a horizontal porous layer. These results are also consistent with the trend of the neutral stability curves for Pe ¼ 0:2 and Pe ¼ 0:1 reported in Fig. 2. In fact, Fig. 2 suggests that the neutral stability curves with n < 1 move upward for smaller and smaller values of Pe, while those for n > 1 move downward. Eventually, when Pe ! 0, the neutral stability curves are pushed either to infinity (if n < 1) or to zero (if n > 1). In a real non-Newtonian fluid, these features are likely to be altered due to the loss of precision in the power law of viscosity as the shear rate tends to zero. Due to this effect, the neutral stability curves are expected to approach a very large (but not infinite) upper bound when Pe ! 0 and n < 1, or a nonzero lower bound when Pe ! 0 and n > 1. The constants A0 and B0 in Eq. (45) can be chosen so that the normalization condition W 00 ð0Þ ¼ 1, given by Eqs. (22c) and (46b) are satisfied. Thus, with q ¼ 1, one has
A0 ¼
1
p
;
B0 ¼
1
pðp2 þ a2 Þ
:
ð50Þ
Eq. (19) yield to first order,
HðzÞ ¼ H0 ðzÞ þ H1 ðzÞ Pe þ OðPe2 Þ; k ¼ k0 þ k1 Pe þ OðPe2 Þ:
k0c ¼ p2
ð44cÞ
W 001 a2 ðn W 1 k0 H1 k1 H0 Þ ¼ 0; 1 H001 H00 a2 H1 þ W 1 þ z W 0 ¼ 0; 2 z ¼ 0; 1 :
W 1 ¼ 0;
H1 ¼ 0:
ð51aÞ ð51bÞ ð51cÞ
Here, W 0 and H0 are known functions given by Eq. (45), with q ¼ 1, and (50). Since the normalization condition W 0 ð0Þ ¼ 1, Eq. (22c), is satisfied by W 0 , one may easily conclude that W 01 ð0Þ ¼ 0. Thus, the boundary value problem (51) is over-determined. One may use, for instance, the boundary condition H1 ð1Þ ¼ 0 to determine the unknown constant k1 . We omit the detailed calculations for the sake of brevity, but the result is k1 ¼ 0. This conclusion is expected since the neutral stability curves in the ða; kÞ plane are
300
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
invariant under a sign change of Pe, as already pointed out in Section 3.1. Hence, the power series expansion of k, Eq. (42), is actually expected to include only even powers of Pe. In order to approximate the neutral stability condition at very small values of Pe one may well use Eqs. (48) and (49), together with Eq. (15), namely
Ra ffi jPejn1
p
ffiffiffi ; ac ffi p 4 n
p2 þ a2 n p2 þ a2
; p ffiffiffi 2 Rac ffi jPejn1 p2 n þ 1 ; a2
kc ffi p2
pffiffiffi 2 nþ1 :
ð52Þ
Fig. 2 shows a comparison between the neutral stability curves drawn by employing Eq. (52). These dashed curves are barely distinguishable from the solid ones, obtained numerically by employing the numerical procedure described in Section 3.2, for Péclet numbers smaller than 0:7. A weak discrepancy is observed when Pe ¼ 1 which becomes quite strong when Pe ¼ 2. Eq. (52) can be employed also for the special case n ! 1, discussed in Section 4.1.2. One may evaluate parameter ~ k ¼ k=n by using Eq. (48) with k ¼ k0 and then by taking the limit n ! 1. Hence, one obtains the extremely simple neutral stability condition
~k ffi p2 þ a2 :
ð53Þ
Fig. 4 shows the neutral stability curve for Pe ¼ 0, Eq. (53), compared with those for small values of Pe; 0:1 or 0:2. These three neutral stability curves appear as almost perfectly overlapped.
Fig. 6. Critical value of k versus n for different Péclet numbers. The line relative to the case Pe ¼ 0 is drawn by employing the asymptotic solution for small Péclet numbers given by Eq. (52).
4.2. Critical values The trends of ac ; Rac and kc versus n for different values of Pe are reported in Figs. 5–7, respectively. The solid lines in Fig. 5 are drawn by employing the numerical solution of the eigenvalue stability problem, Eq. (19), described in Section 3.2. On the other hand, the dashed lines are relative to the asymptotic solution, Eq. (52), for the limiting case Pe 1. Fig. 5 displays an excellent agreement between the analytical asymptotic solution for Pe 1 and the numerical data for either Pe ¼ 0:1; 0:2 or 0:5. An increasingly large discrepancy is displayed in Fig. 5 when Pe ¼ 1 or larger. A similar conclusion is drawn from Figs. 6 and 7. The asymptotic solution (52) shows that either ac and kc do not depend on Pe, so that assessing the reliability of Eq. (52) is just a matter of
Fig. 7. Critical value of a versus n for different Péclet numbers. The line relative to the case Pe ¼ 0 is drawn by employing the asymptotic solution for small Péclet numbers given by Eq. (52).
comparing the curves for Pe ¼ 0 with those for Pe – 0 in either Fig. 6 or Fig. 7. We conclude again that the curves with Pe ¼ 0:1; 0:2 and 0:5 are almost coincident, meaning that these cases are well represented by the asymptotic solution (52). Figs. 6 and 7 display a monotonic increasing trend of kc versus n, and a monotonic decreasing trend of ac versus n. On the other hand, Rac is a non-monotonic function of n displaying a maximum, as it is shown in Fig. 5. This nonmonotonic behaviour is easily explained with Eq. (52). In fact, we obtain that the derivative of Rac with respect to n vanishes for Fig. 5. Critical value of Ra versus n for different Péclet numbers. The dashed lines are relative to the asymptotic solution for small Péclet numbers given by Eq. (52).
jPej ffi exp
1 pffiffiffi : nþ n
ð54Þ
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
301
Fig. 8. Maximum value of Rac : plots of n0 and Rac;max versus jPej. Comparison between the data obtained from the numerical solution, Section 3.2, and the approximate analytical expressions, Eqs. (55) and (56).
Eq. (54) defines the approximate position of the maxima displayed in Fig. 5. This equation also makes it evident that these maxima are allowed only as far as jPej < 1. The maximum value of Rac for a given jPej is approximately obtained by employing Eqs. (52) and (54), namely
pffiffiffiffiffi n0 1 2 pffiffiffiffiffi ; Rac;max ffi p2 ð n0 þ 1Þ exp n0 þ n0
ð55Þ
where n0 is given by
1 n0 ffi 4
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2 4 1 1 ln jPej
ð56Þ
and it is obtained by solving Eq. (54) for n. Fig. 8 shows the plots of n0 and Rac;max versus jPej. The numerical data for these quantities evaluated by the method described in Section 3.2 are compared with the approximate expressions given by Eqs. (55) and (56). The agreement over the range jPej < 0:95 is very good for n0 , with a discrepancy lower than 0.4%, and fairly good for Rac;max , with a discrepancy lower than 4%. The trend of Rac versus n, displayed in Fig. 5 for a fixed Pe and increasingly large values of n, reflects the conclusions stated in Section 4.1.2. When n ! 1; Rac is either 0 if jPej < 1 or infinite if jPej P 1. One may also note that the dependence of Rac on the Péclet number is significantly amplified when the fluid departs from the Newtonian behaviour ðn ¼ 1Þ. 5. Conclusions The instability of vertical throughflow in a horizontal porous layer saturated by a non-Newtonian power-law fluid has been studied by employing a linear analysis. The lower and upper boundary planes are considered as permeable and kept at constant temperatures. Besides the Rayleigh number Ra, there are two governing dimensionless parameters in this problem, the powerlaw index n and the Péclet number Pe. The latter parameter is proportional to the prescribed vertical throughflow velocity. We may have either positive or negative values of Pe meaning upward throughflow or downward throughflow, respectively. The linear stability eigenvalue problem has been solved by employing a numerical solution based on the shooting method. The main features characterising the onset of convective instability are resumed as follows: (1) The neutral stability curve is affected by both n and jPej, while it is unaffected by the sign of Pe. The effect of an increasing n for a given Pe is different for either smaller or larger values
of Pe. Generally speaking, pseudoplastic fluids ðn < 1Þ are more stable than dilatant fluids ðn > 1Þ when Pe is small, while they become more unstable if Pe is large. (2) A mathematically interesting asymptotic condition is achieved when the power-law index n becomes vanishing small. In this limit, the Rayleigh number Ra is a monotonic decreasing function of the wavenumber a along the neutral stability curve. When a ! 1, horizontal asymptotes are reached depending on the Péclet number. Thus, these asymptotes represent the critical value of Ra for the onset of instability. (3) Another asymptotic regime is obtained when n is extremely large. In this regime, the eigenvalue problem can be solved analytically, leading to an expression of the dispersion relation for neutral stability written in terms of Bessel functions. Based on this analytical solution, we concluded that the basic flow is either always stable if jPej P 1 or always unstable if jPej < 1. (4) An analytical asymptotic solution has been obtained for the regime of small Péclet numbers. It has been shown that simple analytical formulations for the neutral stability curve, RaðaÞ, and for the critical values, ac and Rac , are obtained by employing this asymptotic solution. These analytical formulations fit well the numerical data for jPej < 0:5, while they give a fairly acceptable approximation all over the range jPej < 1. (5) If jPej < 1, the critical Rayleigh number displays a maximum when plotted versus n. On the other hand, if jPej P 1; Rac is a monotonic increasing function of n. These items outline the main facets of the non-Newtonian behaviour of vertical throughflow in a horizontal porous layer. We stress that departures from the Newtonian behaviour turned out to enhance the instability effects of the vertical throughflow, by causing an increasingly strong dependence of the critical Rayleigh number on the Péclet number. References [1] C. Zhao, B.E. Hobbs, A. Ord, Convective and Advective Heat Transfer in Geological Systems, Springer-Verlag, Berlin, Germany, 2008. [2] R.A. Wooding, Rayleigh instability of a thermal boundary layer in flow through a porous medium, J. Fluid Mech. 9 (1960) 183–192. [3] F.M. Sutton, Onset of convection in a porous channel with net through flow, Phys. Fluids 13 (1970) 1931–1934. [4] G.M. Homsy, A.E. Sherwood, Convective instabilities in porous media with through flow, AIChE J. 22 (1976) 168–174. [5] M.C. Jones, J.M. Persichetti, Convective instability in packed beds with throughflow, AIChE J. 32 (1986) 1555–1557. [6] D.A. Nield, Convective instability in porous media with throughflow, AIChE J. 33 (1987) 1222–1224. [7] A. Khalili, I.S. Shivakumara, Onset of convection in a porous layer with net through-flow and internal heat generation, Phys. Fluids 10 (1998) 315–317.
302
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
[8] D.A.S. Rees, L. Storesletten, The linear instability of a thermal boundary layer with suction in an anisotropic porous medium, Fluid Dyn. Res. 30 (2002) 155– 168. [9] G.J.M. Pieters, H.M. Schuttelaars, On the nonlinear dynamics of a saline boundary layer formed by throughflow near the surface of a porous medium, Physica D: Nonlinear Phenom. 237 (2008) 3075–3088. [10] D.A.S. Rees, The onset and nonlinear development of vortex instabilities in a horizontal forced convection boundary layer with uniform surface suction, Transp. Porous Media 77 (2009) 243–265. [11] A. Barletta, E. Rossi di Schio, L. Storesletten, Convective roll instabilities of vertical throughflow with viscous dissipation in a horizontal porous layer, Transp. Porous Media 81 (2010) 461–477. [12] D.A. Nield, A.V. Kuznetsov, The onset of convection in a heterogeneous porous medium with vertical throughflow, Transp. Porous Media 88 (2011) 347–355. [13] P.M. Patil, D.A.S. Rees, Linear instability of a horizontal thermal boundary layer formed by vertical throughflow in a porous medium: the effect of local thermal nonequilibrium, Transp. Porous Media 99 (2013) 207–227. [14] C.J. Van Duijn, G.J.M. Pieters, R.A. Wooding, A. Van Der Ploeg, Stability criteria for the vertical boundary layer formed by throughflow near the surface of a porous medium, in: Environmental Mechanics: Water, Mass and Energy Transfer in the Biosphere: The Philip Volume, American Geophysical Union, 2013, pp. 155–169. [15] D.A. Nield, A. Bejan, Convection in Porous Media, fourth ed., Springer-Verlag, New York, 2013. [16] R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, second ed., Wiley, New York, 2006. [17] R.H. Christopher, S. Middleman, Power-law flow through a packed tube, Ind. Eng. Chem. Fundam. 4 (1965) 422–426. [18] A.V. Shenoy, Non-Newtonian fluid heat transfer in porous media, Adv. Heat Transfer 24 (1994) 102–191. [19] V. Di Federico, M. Pinelli, R. Ugarelli, Estimates of effective permeability for non-Newtonian fluid flow in randomly heterogeneous porous media, Stochastic Environ. Res. Risk Assess. 24 (2010) 1067–1076. [20] D.A. Nield, A note on the onset of convection in a layer of a porous medium saturated by a non-Newtonian nanofluid of power-law type, Transp. Porous Media 87 (2011) 121–123. [21] A. Barletta, D.A. Nield, Linear instability of the horizontal throughflow in a plane porous layer saturated by a power-law fluid, Phys. Fluids 23 (2011) 013102.
[22] D.A. Nield, A further note on the onset of convection in a layer of a porous medium saturated by a non-Newtonian fluid of power-law type, Transp. Porous Media 88 (2011) 187–191. [23] Z. Alloui, N. Ben Khelifa, H. Beji, P. Vasseur, Onset of convection in a horizontal porous layer saturated by a power-law fluid, ASME J. Heat Transfer 134 (2012) 092502. [24] V. Ciriello, V. Di Federico, Similarity solutions for flow of non-Newtonian fluids in porous media revisited under parameter uncertainty, Adv. Water Resour. 43 (2012) 38–51. [25] N. Ben Khelifa, Z. Alloui, H. Beji, P. Vasseur, Natural convection in a horizontal porous cavity filled with a non-Newtonian binary fluid of power-law type, J. Non-Newtonian Fluid Mech. 169–170 (2012) 15–25. [26] V. Di Federico, R. Archetti, S. Longo, Spreading of axisymmetric non-Newtonian power-law gravity currents in porous media, J. Non-Newtonian Fluid Mech. 189–190 (2012) 31–39. [27] S. Longo, V. Di Federico, L. Chiapponi, R. Archetti, Experimental verification of power-law non-Newtonian axisymmetric porous gravity currents, J. Fluid Mech. 731 (2013). [28] L.S. de B. Alves, A. Barletta, 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. 62 (2013) 495–506. [29] A. Barletta, L.S. de B. Alves, On Gill’s stability problem for non-Newtonian Darcy’s flow, Int. J. Heat Mass Transf. 79 (2014) 759–768. [30] L.S. de B. Alves, A. Barletta, Convective to absolute instability transition in the Prats flow of a power-law fluid, Int. J. Therm. Sci. 94 (2015) 270–282. [31] I.S. Shivakumara, S. Sureshkumar, Convective instabilities in a viscoelasticfluid-saturated porous medium with throughflow, J. Geophys. Eng. 4 (2007) 104–115. [32] I.S. Shivakumara, S. Sureshkumar, Effects of throughflow and quadratic drag on the stability of a doubly diffusive Oldroyd-B fluid-saturated porous layer, J. Geophys. Eng. 5 (2008) 268–280. [33] K.F. Riley, M.P. Hobson, S.J. Bence, Mathematical Methods for Physics and Engineering: A Comprehensive Guide, Cambridge University Press, 2006. [34] Wolfram Research, Wolfram Language & System Documentation Center, Wolfram Research Inc, Mathematica, Version 10.3, Champaign, IL, 2015..
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Linear instability of the vertical throughflow in a horizontal porous layer saturated by a power-law fluid A. Barletta a,⇑, L. Storesletten b a b
Department of Industrial Engineering, Alma Mater Studiorum Università di Bologna, Viale Risorgimento 2, 40136 Bologna, Italy Department of Mathematics, University of Agder, Postboks 422, 4604 Kristiansand, Norway
a r t i c l e
i n f o
Article history: Received 24 February 2016 Received in revised form 29 March 2016 Accepted 29 March 2016 Available online 16 April 2016 Keywords: Porous medium Non-Newtonian fluid Vertical throughflow Convective instability Rayleigh number
a b s t r a c t The effects of the vertical throughflow of a non-Newtonian power-law fluid on the onset of convective instability in a horizontal porous layer are investigated. The extended Darcy’s model of momentum diffusion is employed together with the Oberbeck–Boussinesq approximation. A stationary basic solution for the vertical throughflow is determined analytically. The basic velocity and temperature fields turn out to be independent of the non-Newtonian rheology. A linear stability analysis is carried out, leading to a fourth-order eigenvalue problem. A numerical solution of the eigenvalue problem is employed to obtain the neutral stability curves and the critical Rayleigh number for the onset of instability. The governing parameters of the transition to instability are the Péclet number associated with the throughflow, and the power-law index of the fluid. These parameters influence the position of the neutral stability curve and also the critical Rayleigh number. The asymptotic cases of absolute pseudoplasticity, absolute dilatancy, and small Péclet number are discussed in detail. The latter case leads to a simple analytical solution that approximates fairly well the numerical data when the Péclet number is smaller than unity. Ó 2016 Published by Elsevier Ltd.
1. Introduction The vertical temperature gradient in Earth’s crust, if large enough, may produce fluid flows with local free convection cells. As a consequence, the porous flow circulation of groundwater can drive soluble minerals in different regions. This may lead to the so-called mineralization problem [1]. Buried nuclear debris or chemical contaminants dispersion can be driven by free convection flow of groundwater in Earth’s crust and thus yield a severe contamination of aquifers. A system involving a porous layer heated from below and saturated by a Newtonian or non-Newtonian fluid is the basic model for the study of these phenomena. Free convection cells produced in such porous layer can be combined with a forced vertical throughflow caused by pressure gradients. For instance, in geophysics, a net vertical through flow may be caused by lithostatic pressure gradients in the lower crust [1]. The modelling of convective instability in a horizontal porous layer heated from below and subject to a vertical throughflow was developed in the pioneering papers by Wooding [2], by Sutton [3], and by Homsy and Sherwood [4]. The flow in the saturated porous medium was studied by adopting a scheme based on ⇑ Corresponding author. E-mail addresses: [email protected] (A. Barletta), leiv.storesletten@ uia.no (L. Storesletten). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2016.03.115 0017-9310/Ó 2016 Published by Elsevier Ltd.
Darcy’s law and on the Oberbeck–Boussinesq approximation. Wooding [2] assumed a layer with infinite thickness, while both Sutton [3]and Homsy and Sherwood [4] considered a layer with a finite thickness. The stability of the vertical throughflow was studied by a linear analysis [2–4], as well as by an energy method nonlinear analysis [4]. Later on, further studies [5–14] were carried out to investigate different effects such as anisotropy, nonlinearity, internal heat generation, viscous dissipation, heterogeneity and local thermal nonequilibrium. An exhaustive survey of the literature on this topic can be found in the book by Nield and Bejan [15]. The main feature of the throughflow is that its effect is stabilising, if compared with the pure Rayleigh–Bénard instability occurring in the porous layer when the throughflow is absent. A slight departure from this behaviour was observed by Jones and Persichetti [5] for the case when the bounding planes of the porous layer are assumed to be isobaric and, thus, permeable to the unstable convective rolls. Many problems in geophysics feature non-Newtonian fluids as, for instance, those involving hydrocarbons and oil reservoirs. In these cases, the analysis of the onset of convection cells in a porous layer implies a departure from the classical Darcy’s law, due to the non-Newtonian rheology of the fluid. A typical non-Newtonian behaviour is the shear-thinning or shear-thickening nature of some fluids. This behaviour can be modelled by the Ostwald–de Waele power law [16]. Christopher and Middleman [17] formulated a
294
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
Nomenclature a A0 ; B0 b ~ b C C ez f F g H Jm K n n0 p
dimensionless wave number, Eq. (18) dimensionless integration constants, Eq. (45) dimensionless parameter, Eq. (25) dimensionless parameter, Eq. (34) dimensionless integration constant, Eq. (36) tortuosity factor, Eq. (4) unit vector in the z direction solution of the Helmholtz equation, Eq. (18) dimensionless function, Eq. (20) gravitational acceleration, with modulus it g thickness of the porous layer Bessel function of first kind and order m permeability power-law index dimensionless parameter, Eq. (56) difference between the pressure and the hydrostatic pressure Pe P´eclet number, Eq. (7) q integer number, Eq. (45) Ra Rayleigh number, Eq. (7) Rac;max maximum critical Rayleigh number, Eq. (55) s dimensionless variable, Eqs. (25) and (34) t time T temperature T0 reference temperature u velocity u; v ; w Cartesian velocity components W dimensionless amplitude of the vertical velocity disturbance, Eq. (17) W 0 ; H0 ; k0 zero order terms of the expansion (42)
modified Darcy’s law to be applied when the saturating fluid has a power-law behaviour. Several further studies on the modified Darcy’s law for power-law fluids were carried out since this groundbreaking paper, providing both experimental and theoretical results to validate this model. Thorough reviews and bibliographies on this topic can be found in the papers by Shenoy [18] and by Di Federico et al. [19]. Several studies have been recently published on flow and convection of power-law fluids in porous media [20–30]. The aim of this paper is to extend the classical analysis of the onset of convective instability in a horizontal porous layer with vertical throughflow to a non-Newtonian fluid of the power-law type. We mention that a similar analysis has been carried out for a viscoelastic fluid in a couple of papers by Shivakumara and Sureshkumar [31,32]. Our linear stability analysis is meant to determine the neutral stability curves for different power-law indexes, n, and different Péclet numbers, Pe, associated with the vertical throughflow. This analysis will be carried out numerically, even if an analytical evaluation of the neutral stability data is done for the limiting cases of absolute pseudoplasticity ðn ! 0Þ, absolute dilatancy ðn ! 1Þ, or very small Péclet numbers. The special case of absent throughflow ðPe ¼ 0Þ yields a perfect agreement with the conclusions drawn by Barletta and Nield [21].
W 1 ; H1 ; k1 first order terms of the expansion (42) x; y; z Cartesian coordinates Ym Bessel function of second kind and order m Greek symbols b thermal expansion coefficient c_ shear rate dz numerical step in the z direction DT temperature difference e dimensionless perturbation parameter, Eq. (13) g dimensionless complex parameter H dimensionless amplitude of the temperature disturbance, Eq. (17) , effective thermal diffusivity k dimensionless parameter, Eq. (15) ~k dimensionless parameter, k=n l consistency factor l effective consistency factor m dimensionless parameter, Eq. (37) q0 reference fluid density r heat capacity ratio s shear stress / porosity Superscripts, subscripts ^ perturbation of the basic flow 0 derivative with respect to z b basic solution c critical value
reference temperature and DT is a positive temperature difference. The velocity field u has Cartesian components ðu; v ; wÞ. 2.1. Governing equations For this kind of non-Newtonian fluid, the shear stress s (SI units are Pa) is a function of the shear rate c_ (SI units are s1 ) given by
s ¼ l ðc_ Þn ;
ð1Þ
where l is the consistency factor and n is the power-law index. In the case of a Newtonian fluid, we have n ¼ 1 and the consistency factor coincides with the dynamic viscosity. The SI measurement units of the consistency factor are Pa sn and, for a Newtonian fluid, the SI units of l are Pa s. Pseudoplastic fluid behaviour is obtained with n < 1, while dilatant fluids are described by n > 1. When a power-law fluid saturates a porous medium and the flow is driven by the thermal buoyancy force, modelled through
2. Mathematical model We consider a plane horizontal porous layer with thickness H, saturated by an Ostwald–de Waele (power law) fluid. Let the vertical axis z be parallel to the gravitational acceleration g, but with opposite direction (see Fig. 1). The boundary planes, z ¼ 0 and z ¼ H, are permeable and kept isothermal at T ¼ T 0 þ DT and T ¼ T 0 , respectively. Here T is the temperature field, T 0 is a
Fig. 1. Sketch of the fluid saturated porous layer and of the boundary conditions.
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
the Oberbeck–Boussinesq approximation, the classical Darcy’s law is generalised as [18,19]
$ u ¼ 0;
l
@T þ u $T ¼ r2 T; @t z ¼ 0 : w ¼ Pe; T ¼ 1;
K
jujn1 u ¼ $p q0 g bðT T 0 Þ;
ð2Þ
where l is the effective consistency factor (SI units are Pa sn m1n ), K is the permeability (SI units are m2 ) and p is the dynamic pressure (SI units are Pa), viz. the local difference between the pressure and the hydrostatic pressure. Moreover, q0 is the fluid density at the reference temperature T 0 ; b is the thermal expansion coefficient of the fluid, and g is the gravity acceleration whose modulus is denoted as g. The ratio l =K depends on the consistency factor l, on the power-law index n, on the permeability K, on the porosity /, and on the tortuosity factor C through the relationship [18,19],
l K
¼ 2Cl
3/ 50K
ðnþ1Þ=2 n 3n þ 1 : n/
ð3Þ
The tortuosity factor C was regarded by Christopher and Middleman [17] as a numerical constant independent of n and equal to 25=12, so that l =K coincides with l=K in the Newtonian limit, n ! 1. Let us assume the validity of the generalised Darcy’s law, Eq. (2). With the assumptions above, the local mass, momentum, and energy balance yield the governing equations
$ u ¼ 0;
l
n1
juj
ð4aÞ
u ¼ $p q0 g bðT T 0 Þ;
ð4bÞ
K @T r þ u $T ¼ , r2 T; @t
ð4cÞ
where t is the time, r is the ratio between the average volumetric heat capacity of the saturated porous medium and the volumetric heat capacity of the fluid, and , is the effective thermal diffusivity of the porous medium. A vertical throughflow is considered, so that the boundary conditions are given by
z¼0:
w ¼ w0 ;
T ¼ T 0 þ DT;
z¼H:
w ¼ w0 ;
T ¼ T 0;
ð5Þ
where w0 is the prescribed vertical throughflow velocity. 2.2. Dimensionless analysis We introduce a dimensionless formulation by scaling the dimensional quantities and operators as follows:
1 H H ðx; y; zÞ ! ðx; y; zÞ; u ¼ ðu; v ; wÞ ! ðu; v ; wÞ ¼ u; H , , T T0 2 ! T; H $ ! $; H r2 ! r2 : DT
, t ! t; rH2 ð6Þ
We also introduce the dimensionless parameters,
q g b DT K H n w0 H Ra ¼ 0 n ; Pe ¼ : l , ,
ð7Þ
The dimensionless parameter Ra is the non-Newtonian version of the Rayleigh number for a saturated porous medium, and Pe is the Péclet number. While Ra can only be positive, the Péclet number can either be positive (upward throughflow) or negative (downward throughflow). By taking the curl of Eq. (4b), the dimensionless formulation of the governing equation (4) and of the boundary conditions (5) is given by
$ jujn1 u ¼ Ra $ ðT ez Þ;
z¼1:
w ¼ Pe;
T ¼ 0;
295
ð8aÞ ð8bÞ ð8cÞ ð8dÞ
where ez is the unit vector in the z direction. 2.3. Basic solution A basic stationary solution of Eq. (8) is given by a uniform throughflow,
v b ¼ 0;
ub ¼ 0;
wb ¼ Pe;
ð9Þ
with a purely vertical temperature gradient dT b ðzÞ=dz. Here, the subscript b denotes the ‘‘basic solution”. The basic temperature distribution T b ðzÞ must satisfy Eqs. (8) and (9), thus yielding the differential equation 2
d Tb dT b Pe ¼ 0; dz2 dz
ð10Þ
together with the boundary conditions
T b ð0Þ ¼ 1;
T b ð1Þ ¼ 0:
ð11Þ
The solution is given by
T b ðzÞ ¼
ePe ePe z : ePe 1
ð12Þ
3. Linear stability analysis 3.1. Disturbance equations We perturb the basic solution, Eqs. (9) and (12), by defining the velocity and temperature disturbances
^; u ¼ ub þ e u
T ¼ T b þ e Tb ;
ð13Þ
where e 1 is a positive perturbation parameter. We substitute Eq. (13) into Eq. (8) and we neglect terms of order higher than e. Then, we obtain
^ @ v^ @ w ^ @u þ ¼ 0; þ @x @y @z
ð14aÞ
n
^ @ v^ @w @ Tb ; ¼k @y @z @y
ð14bÞ
n
^ @u ^ @w @ Tb ¼k ; @x @z @x
ð14cÞ
^ @ v^ @u ¼ 0; @y @x
ð14dÞ
@ Tb dT b @ Tb @ 2 Tb @ 2 Tb @ 2 Tb ^ þw ¼ 2þ 2þ 2; þ Pe @x @y @z @t @z dz
ð14eÞ
^ ¼ 0; w
z ¼ 0; 1 :
Tb ¼ 0;
ð14fÞ
where
k¼
Ra jPejn1
:
ð15Þ
We derive Eq. (14b) with respect to y, and Eq. (14c) with respect to x. Then, by summing the two resulting equations and by taking into b Þ formulation of the stability ^ T account Eq. (14a), we obtain a ðw; problem, namely
296
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
! ! ^ @2w ^ ^ @2w @2w @ 2 Tb @ 2 Tb ; n þ 2 þ 2 ¼k þ @x2 @y2 @x2 @y @z @ Tb dT b @ Tb @ 2 Tb @ 2 Tb @ 2 Tb ^ þw ¼ 2þ 2þ 2; þ Pe @x @y @z @t @z dz b ^ ¼ 0; T ¼ 0; z ¼ 0; 1 : w
ð16aÞ ð16bÞ ð16cÞ
Due to the linearity of our analysis, an arbitrary disturbance can be properly constructed by superposition of normal modes,
^ ¼ WðzÞ egt f ðx; yÞ; w
Tb ¼ HðzÞ egt f ðx; yÞ;
ð17Þ
where g is a complex parameter whose real part expresses the growth rate of the disturbance, while its imaginary part yields the angular frequency. A zero growth rate defines the condition of neutral stability, i.e., the threshold between stability and instability. We assume the validity of the principle of exchange of stabilities discussed, for the Newtonian case ðn ¼ 1Þ, by Homsy and Sherwood [4]. We mention that, for a power-law fluid saturating a porous medium, the principle of exchange of stabilities has been proved rigorously by Alves and Barletta [28] with reference to a case where the throughflow is horizontal. Thus, hereafter, the parameter g will be regarded as real, meaning that neutrally stable modes will be considered as time-independent, i.e., with g ¼ 0. In Eq. (17), f ðx; yÞ is any solution of the two-dimensional Helmholtz equation,
@2f @2f þ þ a2 f ¼ 0; @x2 @y2
ð18Þ
where a > 0 is the wave number. By substituting Eqs. (17) and (18) into Eq. (16) and by setting g ¼ 0, we obtain the eigenvalue problem for the neutrally stable modes, namely
W 00 a2 ðn W k HÞ ¼ 0; 00
0
H Pe H a H þ Pe FðzÞ W ¼ 0; z ¼ 0; 1 : W ¼ 0; H ¼ 0; 2
ð19aÞ ð19bÞ ð19cÞ
where
FðzÞ ¼
1 dT b ePe z ¼ Pe Pe dz e 1
ð20Þ
and primes are used to denote derivatives with respect to z. The fundamental symmetry of the eigenvalue problem (19) is an invariance under the transformation,
z ! 1 z;
Pe ! Pe:
ð21Þ
As a consequence, for a fixed n, changing the sign of Pe means that WðzÞ and HðzÞ are transformed into Wð1 zÞ and Hð1 zÞ, respectively. On the other hand, the eigenvalue Ra as a function of a depends on n and on the absolute value of Pe, but not on its sign. In fact, function RaðaÞ determined by solving Eq. (19) defines the neutral stability curves in the ða; RaÞ plane. Thus one may draw, without any loss of generality, the neutral stability curves only for positive values of Pe, as they depend just on jPej and n. Also the critical values ðac ; Rac Þ, namely the values of ða; RaÞ yielding the point of minimum Ra along the neutral stability curve, depend only on jPej and n. We point out that the symmetry defined by Eq. (21) had been previously revealed for Newtonian fluids ðn ¼ 1Þ [4]. 3.2. Numerical solution The eigenvalue problem controlling the onset of instability, Eq. (19), can be solved numerically. An accurate and efficient procedure is based on the combined use of a Runge–Kutta solver and of the shooting method. One has to formulate Eq. (19) as an initial value problem with an extra unknown parameter, namely
Table 1 Neutral stability values of Ra for n ¼ 0:5 and Pe ¼ 1, with either a ¼ 4 or a ¼ 6: comparison between the Runge–Kutta solver with fixed step-size, dz, and the Runge– Kutta solver with adaptive step-size algorithm. dz
a¼4
a¼6
0.00100 0.00050 0.00010 0.00005 0.00001 Adaptive
29.8029489571230 29.8030190532842 29.8030414820325 29.8030422001112 29.8030424329556 29.8030424388895
36.2054197360977 36.2054769320744 36.2054952405836 36.2054958143489 36.2054959956716 36.2054960084208
W 00 a2 ðn W k HÞ ¼ 0;
ð22aÞ
H00 Pe H0 a2 H þ Pe FðzÞ W ¼ 0; z ¼ 0 : W ¼ 0; W 0 ¼ 1; H ¼ 0;
ð22bÞ
H0 ¼ n:
ð22cÞ
In fact, n is the extra unknown parameter to be determined so that the target conditions Wð1Þ ¼ 0 and Hð1Þ ¼ 0 are satisfied. The initial condition W 0 ð0Þ ¼ 1, in Eq. (22c), can be imposed without any loss of generality as a consequence of the scale invariance of the solutions ðW; HÞ of Eq. (19). The Runge–Kutta method [33] serves to solve Eq. (22) for assigned ða; n; k; Pe; nÞ. For prescribed input parameters ða; n; PeÞ, one can determine the eigenvalue pair ðk; nÞ by imposing the target conditions Wð1Þ ¼ 0 and Hð1Þ ¼ 0. The latter operation can be accomplished by the shooting method and implemented through any root-finding algorithm, such as the Newton–Raphson technique [33]. Function NDSolve of Mathematica 10 (Ó Wolfram Research, Inc.) is employed to implement the Runge–Kutta solver, while the shooting method is applied by using function FindRoot [34]. The default option of NDSolve is based on an adaptive step-size algorithm so that the interval 0 6 z 6 1 is subdivided into steps of variable size, dz. In order to test the accuracy of the solver, the Runge–Kutta method with adaptive step-size is compared with the fixed step-size Runge–Kutta method for decreasing values of dz. In Table 1, the comparison is relative to the neutral stability values of Ra with either a ¼ 4 or a ¼ 6, and the test case is n ¼ 0:5 and Pe ¼ 1. One observes a monotonic convergence of the fixed step-size algorithm to the data obtained by the adaptive step-size algorithm as dz gradually decreases. The comparison between the data obtained with dz ¼ 105 and those with adaptive step-size yields an agreement within no fewer than eight significant figures. 4. Discussion of the results 4.1. The neutral stability curves Fig. 2 displays neutral stability curves on the parametric plane ða; RaÞ. In the cases considered, the Péclet number Pe ranges from 0:1 to 2 and the power-law index n from 0:5 to 2. The neutral stability curve, in each case, confines the region of linear stability laying below the curve from the region of instability above the curve. The neutral stability curve is affected by both parameters Pe and n. The effect of an increasing n for a given Pe is different for either smaller or larger values of Pe. As shown in Fig. 2, when Pe ¼ 0:1 and 0:2, an increasing value of n has a destabilising effect. This behaviour changes gradually when Pe ¼ 0:5 and 0:7, as the transition to instability depends non-monotonically on n. For larger Péclet numbers, such as Pe ¼ 1 and 2, the effect of an increasing n is stabilising. Thus, pseudoplastic fluids are more stable than dilatant fluids when Pe is small, while they become more unstable if Pe is larger. This behaviour is partly a consequence of the scaling with jPejn1 of the Rayleigh number Ra, implicit in the definition (15) of the eigenvalue k. However, this is just a rough
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
297
Fig. 2. Neutral stability curves in the plane ða; RaÞ with different values of Pe and n. The dashed curves are relative to the asymptotic solution for small Péclet numbers given by Eq. (52).
aspect of the actual trend, due to the sensible dependence of k on n at neutral stability. That things are more complicated than the law
collapses into the single fourth-order differential equation in the unknown W given by
Ra jPejn1 is made clear by the frame with Pe ¼ 1 in Fig. 2, where Ra displays an evident change with n. The dashed curves, drawn in Fig. 2, are relative to the asymptotic case of small Péclet numbers which will be discussed in Section 4.1.3.
W 0000 Pe W 000 a2 W 00 a2 k Pe FðzÞ W ¼ 0;
4.1.1. The asymptotic case of absolute pseudoplasticity, n ! 0 A mathematically interesting asymptotic condition is achieved when the power-law index becomes vanishingly small. When n ! 0, the system of two second-order differential equation (19)
z ¼ 0; 1 :
W ¼ 0;
W 00 ¼ 0:
ð23aÞ ð23bÞ
The trend of the neutral stability curves plotted in Fig. 2 suggests that, with n > 0; Ra tends to infinity when a ! 1. Incidentally, this is a usual feature, in the domain of small wave lengths, for the stability analysis of a Rayleigh–Bénard setup. On the other hand, if one considers Eq. (23), the behaviour for a ! 1 is described by
298
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
W 00 þ k Pe FðzÞ W ¼ 0;
ð24aÞ
z ¼ 0; 1 :
ð24bÞ
W ¼ 0:
By defining, 2
b ¼
k ; PeðePe 1Þ
s ¼ ePe z ;
ð25Þ
Eq. (24) can be rewritten as
ð26aÞ ð26bÞ
The general solution of Eq. (26a) can be expressed in terms of Bessel functions as
pffiffi pffiffi W ¼ J0 2 b s þ C Y 0 2 b s :
ð27Þ
The constants b and C are determined so that the boundary conditions (26b) are satisfied. Thus, on account of Eq. (26b), one can evaluate b and C for a given Péclet number by solving the system
J0 ð2 bÞ þ C Y0 ð2 bÞ ¼ 0; J0 2 b ePe=2 þ C Y0 2 b ePe=2 ¼ 0;
ð28Þ
Y0 ð2 bÞ J0 2 b ePe=2 J0 ð2 bÞ Y0 2 b ePe=2 ¼ 0:
ð30Þ
Ra ¼ 1:445796;
ð31Þ
which represents the smallest asymptotic value displayed for a neutral stability curve. Thus, one may imagine that the neutral stability curves of Fig. 3 become closer and closer as Pe increases, and their asymptotes tend to the value reported in Eq. (31) when Pe ! 1. One may well say that the value of Ra given by Eq. (31) is the lower bound to the critical values of Ra when n ! 0. 4.1.2. The asymptotic case of absolute dilatancy, n ! 1 Another asymptotic regime is obtained when n is extremely large. No doubt that this condition is unlikely for real fluids, but nonetheless it deserves some mathematical interest. One may take the limit n ! 1 of Eq. (19) by defining
k ¼ n ~k
or, equivalently, one can obtain b by solving the equation
pffiffiffiffiffiffi J0 2 Ra ¼ 0:
On evaluating the lowest root of Eq. (30), one obtains
2
d W dW 2 þ b W ¼ 0; þ ds2 ds s ¼ 1; ePe : W ¼ 0:
s
One may wonder whether the neutral stability curves cross the axis Ra ¼ 0 for a sufficiently high Péclet number. This is impossible physically as it would mean that thermal instability can arise even when the upper boundary is hotter than the lower boundary. Moreover, this is also impossible mathematically as, on account of Eq. (25) in the limit Pe ! 1, Eq. (29) tends to
ð29Þ
Eqs. (25) and (29) can be used to evaluate the asymptotic value of Ra, when a ! 1, for a prescribed Pe. Fig. 3 shows the neutral stability curves in the plane ða; RaÞ for the limiting case n ! 0. For a given Péclet number, this figure reveals that Ra is a monotonic decreasing function of a along the neutral stability curve. When a ! 1, an asymptote is reached where the value of Ra is that obtained by solving Eqs. (25) and (29). These asymptotes are displayed in Fig. 3 as dashed lines. Thus, one may identify the asymptotic values of Ra when a ! 1 with the critical values of Ra for the onset of instability.
ð32Þ
and by assuming ~ k Oð1Þ. Thus, Eq. (19a) yields W ¼ ~ k H and Eqs. (19b) and (19c) yield
W 00 Pe W 0 a2 W þ ~k Pe FðzÞ W ¼ 0;
ð33aÞ
z ¼ 0; 1 :
ð33bÞ
W ¼ 0:
In analogy to the procedure for the analytical solution of Eq. (24), we define
~2 ¼ b
k~ ; PeðePe 1Þ
s ¼ ePe z :
ð34Þ
Eq. (33) can be rewritten as
s2
2 2 d W ~2 s a W ¼ 0; þ b ds2 Pe2
s ¼ 1; ePe :
W ¼ 0:
ð35aÞ ð35bÞ
The general solution of Eq. (35a) can be expressed again in terms of Bessel functions,
W¼
pffiffii pffiffih ~ pffiffi ~ s ; s Jm 2 b s þ C Y m 2 b
ð36Þ
where C is an integration constant and
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 4 a2 m ¼ 1 þ 2: Pe
ð37Þ
The general solution (36) satisfies the boundary conditions (35b) if the following dispersion relation holds:
~ J 2b ~ ePe=2 J 2 b ~ Ym 2 b ~ ePe=2 ¼ 0: Ym 2 b m m
ð38Þ
Eqs. (37) and (38) can be employed to draw the neutral stability curves in the plane ða; ~ kÞ for different values of Pe. Fig. 4 displays these neutral stability curves. This figure reveals that the cases having Pe < 0:2 cannot be distinguished graphically, while slight differences can be seen for Pe ¼ 0:5 and 0:7 becoming more significant ~ are finite so when Pe ¼ 1 or 2. One may note that the values of k Fig. 3. Asymptotic case n ! 0: neutral stability curves in the plane ða; RaÞ with different values of Pe. The dashed lines denote the critical value of Ra attained for a ! 1 and evaluated by solving Eq. (29).
that, for n 1, one may approximate the neutral stability values of Ra by invoking Eqs. (15) and (32), namely
Ra ¼ n ~k jPejn1 :
ð39Þ
299
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
W 0 ðzÞ ¼ A0 sinðq pzÞ;
H0 ðzÞ ¼ B0 sinðq pzÞ;
q ¼ 1; 2; 3; . . . ; ð45Þ
where A0 and B0 are constants. By substituting Eq. (45) into Eq. (44), we obtain a linear system in the unknowns A0 and B0 ,
A0 q2 p2 þ a2 n a2 k0 B0 ¼ 0; 2 2 A0 B0 q p þ a2 ¼ 0:
ð46aÞ ð46bÞ
Nonzero solutions ðA0 ; B0 Þ are allowed only if the determinant of the coefficient matrix is zero, namely if
k0 ¼
q2 p2 þ a2 n q2 p2 þ a2 : a2
ð47Þ
Eq. (47) yields the neutral stability condition in the limit Pe ! 0, i.e. when the vertical throughflow is switched off. The lowest branch of neutral stability is that with q ¼ 1,
k0 ¼
p2 þ a2 n p2 þ a2 a2
;
ð48Þ
yielding a minimum of k0 for the critical values ac and k0c given by
p
ffiffiffi ; ac ¼ p 4 n ~ with Fig. 4. Asymptotic case n ! 1: neutral stability curves in the plane ða; kÞ different values of Pe. The curve for Pe ¼ 0 is drawn by employing Eq. (53).
An immediate consequence of Eq. (39) is that, at neutral stability,
lim Ra ¼ 0
ð40Þ
n!1
for every a if jPej < 1, while
lim Ra ¼ 1
ð41Þ
n!1
for every a if jPej P 1. Eqs. (40) and (41) allow us to extrapolate the trend of the neutral stability curves displayed in Fig. 2 when n increases to large values far above n ¼ 2. If 0 < Pe < 1, these curves are expected to move downward in the plane ða; RaÞ until they eventually collapse to the horizontal line Ra ¼ 0 when n ! 1. If Pe P 1, the neutral stability curves with n > 2 are expected to gradually move upward, as n increases, and eventually tend to infinity when n ! 1. In other words if n ! 1, Eqs. (40) and (41) predict instability for every value of Ra when jPej < 1, or stability for every value of Ra when jPej P 1. 4.1.3. An asymptotic solution for Pe 1 An asymptotic solution can be obtained in the limiting case of a very small Péclet number. We start from Eq. (19) and we expand the eigenfunctions W and H, as well as the parameter k, in powers of Pe, namely
WðzÞ ¼ W 0 ðzÞ þ W 1 ðzÞ Pe þ OðPe2 Þ; ð42Þ
From Eq. (20) we also obtain
1 þ OðPe2 Þ: Pe FðzÞ ¼ 1 þ Pe z 2
ð43Þ
Therefore, Eq. (19) yield to order zero,
W 000 a2 ðn W 0 k0 H0 Þ ¼ 0; 00 0
H a H0 þ W 0 ¼ 0; z ¼ 0; 1 :
W 0 ¼ 0;
ð44aÞ ð44bÞ
2
H0 ¼ 0:
Eq. (44) can be solved analytically with
pffiffiffi 2 nþ1 :
ð49Þ
On account of Eq. (15), we point out that Eq. (49) means that, when Pe ! 0, the critical value of Ra is either 0, if n > 1, or infinite, if n < 1. Physically, this means that dilatant fluids are always unstable when the vertical throughflow is switched off, while pseudoplastic fluids are always stable. Eqs. (15) and (49) imply that Newtonian fluids are a special case where Rac ¼ 4 p2 , when Pe ! 0. These results are perfectly consistent with the conclusions drawn by Barletta and Nield [21] with reference to the Rayleigh–Bénard instability of a power-law fluid saturating a horizontal porous layer. These results are also consistent with the trend of the neutral stability curves for Pe ¼ 0:2 and Pe ¼ 0:1 reported in Fig. 2. In fact, Fig. 2 suggests that the neutral stability curves with n < 1 move upward for smaller and smaller values of Pe, while those for n > 1 move downward. Eventually, when Pe ! 0, the neutral stability curves are pushed either to infinity (if n < 1) or to zero (if n > 1). In a real non-Newtonian fluid, these features are likely to be altered due to the loss of precision in the power law of viscosity as the shear rate tends to zero. Due to this effect, the neutral stability curves are expected to approach a very large (but not infinite) upper bound when Pe ! 0 and n < 1, or a nonzero lower bound when Pe ! 0 and n > 1. The constants A0 and B0 in Eq. (45) can be chosen so that the normalization condition W 00 ð0Þ ¼ 1, given by Eqs. (22c) and (46b) are satisfied. Thus, with q ¼ 1, one has
A0 ¼
1
p
;
B0 ¼
1
pðp2 þ a2 Þ
:
ð50Þ
Eq. (19) yield to first order,
HðzÞ ¼ H0 ðzÞ þ H1 ðzÞ Pe þ OðPe2 Þ; k ¼ k0 þ k1 Pe þ OðPe2 Þ:
k0c ¼ p2
ð44cÞ
W 001 a2 ðn W 1 k0 H1 k1 H0 Þ ¼ 0; 1 H001 H00 a2 H1 þ W 1 þ z W 0 ¼ 0; 2 z ¼ 0; 1 :
W 1 ¼ 0;
H1 ¼ 0:
ð51aÞ ð51bÞ ð51cÞ
Here, W 0 and H0 are known functions given by Eq. (45), with q ¼ 1, and (50). Since the normalization condition W 0 ð0Þ ¼ 1, Eq. (22c), is satisfied by W 0 , one may easily conclude that W 01 ð0Þ ¼ 0. Thus, the boundary value problem (51) is over-determined. One may use, for instance, the boundary condition H1 ð1Þ ¼ 0 to determine the unknown constant k1 . We omit the detailed calculations for the sake of brevity, but the result is k1 ¼ 0. This conclusion is expected since the neutral stability curves in the ða; kÞ plane are
300
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
invariant under a sign change of Pe, as already pointed out in Section 3.1. Hence, the power series expansion of k, Eq. (42), is actually expected to include only even powers of Pe. In order to approximate the neutral stability condition at very small values of Pe one may well use Eqs. (48) and (49), together with Eq. (15), namely
Ra ffi jPejn1
p
ffiffiffi ; ac ffi p 4 n
p2 þ a2 n p2 þ a2
; p ffiffiffi 2 Rac ffi jPejn1 p2 n þ 1 ; a2
kc ffi p2
pffiffiffi 2 nþ1 :
ð52Þ
Fig. 2 shows a comparison between the neutral stability curves drawn by employing Eq. (52). These dashed curves are barely distinguishable from the solid ones, obtained numerically by employing the numerical procedure described in Section 3.2, for Péclet numbers smaller than 0:7. A weak discrepancy is observed when Pe ¼ 1 which becomes quite strong when Pe ¼ 2. Eq. (52) can be employed also for the special case n ! 1, discussed in Section 4.1.2. One may evaluate parameter ~ k ¼ k=n by using Eq. (48) with k ¼ k0 and then by taking the limit n ! 1. Hence, one obtains the extremely simple neutral stability condition
~k ffi p2 þ a2 :
ð53Þ
Fig. 4 shows the neutral stability curve for Pe ¼ 0, Eq. (53), compared with those for small values of Pe; 0:1 or 0:2. These three neutral stability curves appear as almost perfectly overlapped.
Fig. 6. Critical value of k versus n for different Péclet numbers. The line relative to the case Pe ¼ 0 is drawn by employing the asymptotic solution for small Péclet numbers given by Eq. (52).
4.2. Critical values The trends of ac ; Rac and kc versus n for different values of Pe are reported in Figs. 5–7, respectively. The solid lines in Fig. 5 are drawn by employing the numerical solution of the eigenvalue stability problem, Eq. (19), described in Section 3.2. On the other hand, the dashed lines are relative to the asymptotic solution, Eq. (52), for the limiting case Pe 1. Fig. 5 displays an excellent agreement between the analytical asymptotic solution for Pe 1 and the numerical data for either Pe ¼ 0:1; 0:2 or 0:5. An increasingly large discrepancy is displayed in Fig. 5 when Pe ¼ 1 or larger. A similar conclusion is drawn from Figs. 6 and 7. The asymptotic solution (52) shows that either ac and kc do not depend on Pe, so that assessing the reliability of Eq. (52) is just a matter of
Fig. 7. Critical value of a versus n for different Péclet numbers. The line relative to the case Pe ¼ 0 is drawn by employing the asymptotic solution for small Péclet numbers given by Eq. (52).
comparing the curves for Pe ¼ 0 with those for Pe – 0 in either Fig. 6 or Fig. 7. We conclude again that the curves with Pe ¼ 0:1; 0:2 and 0:5 are almost coincident, meaning that these cases are well represented by the asymptotic solution (52). Figs. 6 and 7 display a monotonic increasing trend of kc versus n, and a monotonic decreasing trend of ac versus n. On the other hand, Rac is a non-monotonic function of n displaying a maximum, as it is shown in Fig. 5. This nonmonotonic behaviour is easily explained with Eq. (52). In fact, we obtain that the derivative of Rac with respect to n vanishes for Fig. 5. Critical value of Ra versus n for different Péclet numbers. The dashed lines are relative to the asymptotic solution for small Péclet numbers given by Eq. (52).
jPej ffi exp
1 pffiffiffi : nþ n
ð54Þ
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
301
Fig. 8. Maximum value of Rac : plots of n0 and Rac;max versus jPej. Comparison between the data obtained from the numerical solution, Section 3.2, and the approximate analytical expressions, Eqs. (55) and (56).
Eq. (54) defines the approximate position of the maxima displayed in Fig. 5. This equation also makes it evident that these maxima are allowed only as far as jPej < 1. The maximum value of Rac for a given jPej is approximately obtained by employing Eqs. (52) and (54), namely
pffiffiffiffiffi n0 1 2 pffiffiffiffiffi ; Rac;max ffi p2 ð n0 þ 1Þ exp n0 þ n0
ð55Þ
where n0 is given by
1 n0 ffi 4
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2 4 1 1 ln jPej
ð56Þ
and it is obtained by solving Eq. (54) for n. Fig. 8 shows the plots of n0 and Rac;max versus jPej. The numerical data for these quantities evaluated by the method described in Section 3.2 are compared with the approximate expressions given by Eqs. (55) and (56). The agreement over the range jPej < 0:95 is very good for n0 , with a discrepancy lower than 0.4%, and fairly good for Rac;max , with a discrepancy lower than 4%. The trend of Rac versus n, displayed in Fig. 5 for a fixed Pe and increasingly large values of n, reflects the conclusions stated in Section 4.1.2. When n ! 1; Rac is either 0 if jPej < 1 or infinite if jPej P 1. One may also note that the dependence of Rac on the Péclet number is significantly amplified when the fluid departs from the Newtonian behaviour ðn ¼ 1Þ. 5. Conclusions The instability of vertical throughflow in a horizontal porous layer saturated by a non-Newtonian power-law fluid has been studied by employing a linear analysis. The lower and upper boundary planes are considered as permeable and kept at constant temperatures. Besides the Rayleigh number Ra, there are two governing dimensionless parameters in this problem, the powerlaw index n and the Péclet number Pe. The latter parameter is proportional to the prescribed vertical throughflow velocity. We may have either positive or negative values of Pe meaning upward throughflow or downward throughflow, respectively. The linear stability eigenvalue problem has been solved by employing a numerical solution based on the shooting method. The main features characterising the onset of convective instability are resumed as follows: (1) The neutral stability curve is affected by both n and jPej, while it is unaffected by the sign of Pe. The effect of an increasing n for a given Pe is different for either smaller or larger values
of Pe. Generally speaking, pseudoplastic fluids ðn < 1Þ are more stable than dilatant fluids ðn > 1Þ when Pe is small, while they become more unstable if Pe is large. (2) A mathematically interesting asymptotic condition is achieved when the power-law index n becomes vanishing small. In this limit, the Rayleigh number Ra is a monotonic decreasing function of the wavenumber a along the neutral stability curve. When a ! 1, horizontal asymptotes are reached depending on the Péclet number. Thus, these asymptotes represent the critical value of Ra for the onset of instability. (3) Another asymptotic regime is obtained when n is extremely large. In this regime, the eigenvalue problem can be solved analytically, leading to an expression of the dispersion relation for neutral stability written in terms of Bessel functions. Based on this analytical solution, we concluded that the basic flow is either always stable if jPej P 1 or always unstable if jPej < 1. (4) An analytical asymptotic solution has been obtained for the regime of small Péclet numbers. It has been shown that simple analytical formulations for the neutral stability curve, RaðaÞ, and for the critical values, ac and Rac , are obtained by employing this asymptotic solution. These analytical formulations fit well the numerical data for jPej < 0:5, while they give a fairly acceptable approximation all over the range jPej < 1. (5) If jPej < 1, the critical Rayleigh number displays a maximum when plotted versus n. On the other hand, if jPej P 1; Rac is a monotonic increasing function of n. These items outline the main facets of the non-Newtonian behaviour of vertical throughflow in a horizontal porous layer. We stress that departures from the Newtonian behaviour turned out to enhance the instability effects of the vertical throughflow, by causing an increasingly strong dependence of the critical Rayleigh number on the Péclet number. References [1] C. Zhao, B.E. Hobbs, A. Ord, Convective and Advective Heat Transfer in Geological Systems, Springer-Verlag, Berlin, Germany, 2008. [2] R.A. Wooding, Rayleigh instability of a thermal boundary layer in flow through a porous medium, J. Fluid Mech. 9 (1960) 183–192. [3] F.M. Sutton, Onset of convection in a porous channel with net through flow, Phys. Fluids 13 (1970) 1931–1934. [4] G.M. Homsy, A.E. Sherwood, Convective instabilities in porous media with through flow, AIChE J. 22 (1976) 168–174. [5] M.C. Jones, J.M. Persichetti, Convective instability in packed beds with throughflow, AIChE J. 32 (1986) 1555–1557. [6] D.A. Nield, Convective instability in porous media with throughflow, AIChE J. 33 (1987) 1222–1224. [7] A. Khalili, I.S. Shivakumara, Onset of convection in a porous layer with net through-flow and internal heat generation, Phys. Fluids 10 (1998) 315–317.
302
A. Barletta, L. Storesletten / International Journal of Heat and Mass Transfer 99 (2016) 293–302
[8] D.A.S. Rees, L. Storesletten, The linear instability of a thermal boundary layer with suction in an anisotropic porous medium, Fluid Dyn. Res. 30 (2002) 155– 168. [9] G.J.M. Pieters, H.M. Schuttelaars, On the nonlinear dynamics of a saline boundary layer formed by throughflow near the surface of a porous medium, Physica D: Nonlinear Phenom. 237 (2008) 3075–3088. [10] D.A.S. Rees, The onset and nonlinear development of vortex instabilities in a horizontal forced convection boundary layer with uniform surface suction, Transp. Porous Media 77 (2009) 243–265. [11] A. Barletta, E. Rossi di Schio, L. Storesletten, Convective roll instabilities of vertical throughflow with viscous dissipation in a horizontal porous layer, Transp. Porous Media 81 (2010) 461–477. [12] D.A. Nield, A.V. Kuznetsov, The onset of convection in a heterogeneous porous medium with vertical throughflow, Transp. Porous Media 88 (2011) 347–355. [13] P.M. Patil, D.A.S. Rees, Linear instability of a horizontal thermal boundary layer formed by vertical throughflow in a porous medium: the effect of local thermal nonequilibrium, Transp. Porous Media 99 (2013) 207–227. [14] C.J. Van Duijn, G.J.M. Pieters, R.A. Wooding, A. Van Der Ploeg, Stability criteria for the vertical boundary layer formed by throughflow near the surface of a porous medium, in: Environmental Mechanics: Water, Mass and Energy Transfer in the Biosphere: The Philip Volume, American Geophysical Union, 2013, pp. 155–169. [15] D.A. Nield, A. Bejan, Convection in Porous Media, fourth ed., Springer-Verlag, New York, 2013. [16] R.B. Bird, W.E. Stewart, E.N. Lightfoot, Transport Phenomena, second ed., Wiley, New York, 2006. [17] R.H. Christopher, S. Middleman, Power-law flow through a packed tube, Ind. Eng. Chem. Fundam. 4 (1965) 422–426. [18] A.V. Shenoy, Non-Newtonian fluid heat transfer in porous media, Adv. Heat Transfer 24 (1994) 102–191. [19] V. Di Federico, M. Pinelli, R. Ugarelli, Estimates of effective permeability for non-Newtonian fluid flow in randomly heterogeneous porous media, Stochastic Environ. Res. Risk Assess. 24 (2010) 1067–1076. [20] D.A. Nield, A note on the onset of convection in a layer of a porous medium saturated by a non-Newtonian nanofluid of power-law type, Transp. Porous Media 87 (2011) 121–123. [21] A. Barletta, D.A. Nield, Linear instability of the horizontal throughflow in a plane porous layer saturated by a power-law fluid, Phys. Fluids 23 (2011) 013102.
[22] D.A. Nield, A further note on the onset of convection in a layer of a porous medium saturated by a non-Newtonian fluid of power-law type, Transp. Porous Media 88 (2011) 187–191. [23] Z. Alloui, N. Ben Khelifa, H. Beji, P. Vasseur, Onset of convection in a horizontal porous layer saturated by a power-law fluid, ASME J. Heat Transfer 134 (2012) 092502. [24] V. Ciriello, V. Di Federico, Similarity solutions for flow of non-Newtonian fluids in porous media revisited under parameter uncertainty, Adv. Water Resour. 43 (2012) 38–51. [25] N. Ben Khelifa, Z. Alloui, H. Beji, P. Vasseur, Natural convection in a horizontal porous cavity filled with a non-Newtonian binary fluid of power-law type, J. Non-Newtonian Fluid Mech. 169–170 (2012) 15–25. [26] V. Di Federico, R. Archetti, S. Longo, Spreading of axisymmetric non-Newtonian power-law gravity currents in porous media, J. Non-Newtonian Fluid Mech. 189–190 (2012) 31–39. [27] S. Longo, V. Di Federico, L. Chiapponi, R. Archetti, Experimental verification of power-law non-Newtonian axisymmetric porous gravity currents, J. Fluid Mech. 731 (2013). [28] L.S. de B. Alves, A. Barletta, 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. 62 (2013) 495–506. [29] A. Barletta, L.S. de B. Alves, On Gill’s stability problem for non-Newtonian Darcy’s flow, Int. J. Heat Mass Transf. 79 (2014) 759–768. [30] L.S. de B. Alves, A. Barletta, Convective to absolute instability transition in the Prats flow of a power-law fluid, Int. J. Therm. Sci. 94 (2015) 270–282. [31] I.S. Shivakumara, S. Sureshkumar, Convective instabilities in a viscoelasticfluid-saturated porous medium with throughflow, J. Geophys. Eng. 4 (2007) 104–115. [32] I.S. Shivakumara, S. Sureshkumar, Effects of throughflow and quadratic drag on the stability of a doubly diffusive Oldroyd-B fluid-saturated porous layer, J. Geophys. Eng. 5 (2008) 268–280. [33] K.F. Riley, M.P. Hobson, S.J. Bence, Mathematical Methods for Physics and Engineering: A Comprehensive Guide, Cambridge University Press, 2006. [34] Wolfram Research, Wolfram Language & System Documentation Center, Wolfram Research Inc, Mathematica, Version 10.3, Champaign, IL, 2015.