Transport in Porous Media 53: 245–263, 2003. © 2003 Kluwer Academic Publishers. Printed in the Netherlands.
245
Non-Darcian Effects on the Onset of Convection in a Porous Layer with Throughflow A. KHALILI1 and I. S. SHIVAKUMARA2
1 Max Planck Institute for Marine Microbiology, Celsiusstr. 1, 28359 Bremen, Germany. e-mail:
[email protected] 2 UGC-DSA Centre in Fluid Mechanics, Department of Mathematics, Bangalore University, Bangalore 560001, India
(Received: 21 January 2002; in final form: 22 July 2002) Abstract. The Brinkman extended Darcy model including Lapwood and Forchheimer inertia terms with fluid viscosity being different from effective viscosity is employed to investigate the effect of vertical throughflow on thermal convective instabilities in a porous layer. Three different types of boundary conditions (free–free, rigid–rigid and rigid–free) are considered which are either conducting or insulating to temperature perturbations. The Galerkin method is used to calculate the critical Rayleigh numbers for conducting boundaries, while closed form solutions are achieved for insulating boundaries. The relative importance of inertial resistance on convective instabilities is investigated in detail. In the case of rigid–free boundaries, it is found that throughflow is destabilizing depending on the choice of physical parameters and the model used. Further, it is noted that an increase in viscosity ratio delays the onset of convection. Standard results are also obtained as particular cases from the general model presented here. Key words: throughflow, convective instability, viscosity ratio, hydrothermal vents.
1. Introduction Thermal convective instability in porous media has been studied widely in recent decades because of its occurrence in nature and in many technological applications, such as, thermal insulation engineering, enhanced recovery of petroleum reservoirs and deep sea hydrodynamics to name a few. In the classical problem, no base flow through the porous environment is considered. However, the most naturally occurring problems, for example, hydrothermal vents (Holm, 1992; Dando et al., 2000) and underground spreading of chemical wastes involve non-isothermal flow of fluids through porous media, referred to as throughflow. The pioneering works considering the throughflow effects on convective instability in porous media go back to Wooding (1960), Homsy and Sherwood (1976), Jones and Persichetti (1986) and Nield (1987). An ample documentation of the research progress made concerning this topic till 1990 can be found in the book of Nield and Bejan (1992). Recently, Khalili and Shivakumara (1998) have discussed the effect of throughflow
246
A. KHALILI AND I. S. SHIVAKUMARA
and internal heat generation on the onset of convection in a porous layer with throughflow. It was found that throughflow destabilizes the system even if the boundaries are of the same type; a result which is not true in the absence of internal heat source. All the studies cited above have considered Darcy flow model which neglects the inertia and viscous effects. Nevertheless, these effects are important, especially when a strong throughflow (high Péclet number) exists (Vafai and Tien, 1981; Vafai and Kim, 1990). Further, Givler and Altobelli (1994) have observed in their experiments that the effective viscosity is different from fluid viscosity. Therefore, in the present study, a general flow model is considered to account for the non-Darcian effects on convective instability in a porous layer with transverse throughflow. In addition, the ratio of viscosities appears as a separate parameter. The model considered allows us to deduce results for ordinary viscous fluids and Darcy model as limiting cases. 2. Mathematical Formulation We consider a fluid-saturated horizontal porous layer of thickness d with an imposed flow of magnitude W0 which is either gravity aligned or anti-gravity in its direction (see Figure 1). The top and the bottom boundaries are maintained at different constant temperatures, which are low at the top (T1 ) and high at the bottom (T0 ). A rectangular coordinate system is taken with the origin at the bottom of the layer and the z-axis vertically upward. The governing equations for continuity, momentum and energy with Boussinesq approximation are, respectively ∇ · q¯ = 0, ρ0
1 µ ρ0 cf 1 ∂ q¯ ¯ q, ¯ + 2 (q¯ · ∇)q¯ = −∇p − ρg kˆ − q¯ + µ∇ ˜ 2 q¯ − √ |q| ∂t k k
Figure 1. Physical configuration.
(1) (2)
247
NON-DARCIAN EFFECTS
A
∂T + (q¯ · ∇)T = κ∇ 2 T , ∂t
ρ = ρ0 [1 − α(T − T0 )],
(3) (4)
where q¯ = (u, v, w), T and p represent the flow velocity, temperature and pressure respectively, while µ, ˜ µ, k, , cf , g, α and A correspondingly denote the effective viscosity, fluid viscosity, permeability, porosity, form drag constant, gravitational acceleration, thermal expansion coefficient and ratio of the effective thermal capacitance to the fluid thermal capacitance. Further, kˆ is the unit vector in z-direction and ρ0 is the reference density. Note that the basic steady state is not quiescent and the temperature distribution is given by
1 − eW0 z/κ . Tb (z) = T0 + (T0 − T1 ) 1 − eW0 d/κ
(5)
The nonlinearity in basic temperature distribution is due to throughflow and it has a profound effect on the stability of the system. In order to investigate the stability of ¯ p, T ) are introduced as the basic state, infinitesimal perturbations (q¯ , p , T ) to (q, q¯ = qb + q¯ , p = pb + p and T = Tb + T . Following the standard linear stability analysis procedure, and non-dimensionalizing the equations as (after dropping the ¯ T∗ = primes) t ∗ = t/(d 2 /κ), (x ∗ , y ∗ , z∗ ) = (x/d, y/d, z/d), q¯ ∗ = q/(κ/d), T /(T0 − T1 ), we get (asterisks are dropped) √ Pe ∂ 1 ∂ + 2 + cf |Pe| Da−1 + Da−1 − ∇ 2 ∇ 2 w = Ra ∇h2 T , Pr ∂t ∂z
A
∂T ∂T + Pe + f (z)w = ∇ 2 T . ∂t ∂z
(6)
(7)
In the above equations Ra = αg(T0 − T1 )d 3 /νκ is the Rayleigh number, Pr = ν/κ is the Prandtl number, Da = k/d 2 is the Darcy number, Pe = W0 d/κ is the Péclet number, = µ/µ ˜ is the ratio of viscosities, ∇h2 = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 is the horizontal Laplacian operator, ∇ 2 = ∇h2 + ∂ 2 /∂z2 is the Laplacian operator and f (z) =
ePez ∂Tb = −Pe Pe . ∂z e −1
(8)
These equations are supplemented by boundary conditions at z = 0 and z = 1. The boundaries are assumed to be either conducting or insulating to temperature perturbations and either rigid but permeable or free of tangential stresses. The
248
A. KHALILI AND I. S. SHIVAKUMARA
following types of boundary conditions are considered: (i) free–free boundaries (a) conducting: ∂ 2w w= = T = 0 at z = 0, 1, (9) ∂z2 (b) insulating: ∂T ∂ 2w = 0 at z = 0, 1, (10) = w= ∂z2 ∂z (ii) rigid–rigid boundaries (a) conducting: ∂w = T = 0 at z = 0, 1, (11) w= ∂z (b) insulating: ∂T ∂w = = 0 at z = 0, 1, (12) w= ∂z ∂z (iii) lower rigid and upper free boundaries (a) conducting: ∂ 2w ∂w = T = 0 at z = 0, w= = T = 0 at z = 1, w= ∂z ∂z2 (13) (b) insulating: ∂T ∂ 2w ∂T ∂w = = 0 at z = 0, w= = 0 at z = 1. = w= 2 ∂z ∂z ∂z ∂z (14) Noting that the principle of exchange of instability holds (Homsy and Sherwood, 1976), a normal mode representation of the form (w, T ) = [W (z), T (z)] ei(lx+my)
(15)
is assumed, where W (z) and T (z) are the related amplitudes. Substituting Equation (15) in Equations (6) and (7) results in the system of equations Pe (16) (D 2 − a 2 ) (D 2 − a 2 ) − G − 2 D W = Ra a 2 T , Pr (D 2 − a 2 )T − Pe DT = f (z)W,
(17)
√ where D = d/dz and G = Da−1 + cf |Pe| Da−1 /Pr. The boundary conditions now take the form: (i) free–free boundaries (a) conducting: W = D 2 W = T = 0 at z = 0, 1, (b) insulating: W = D 2 W = DT = 0 at z = 0, 1,
(18) (19)
NON-DARCIAN EFFECTS
(ii) rigid–rigid boundaries (a) conducting: W = DW = T = 0 at z = 0, 1, (b) insulating: W = DW = DT = 0 at z = 0, 1, (iii) lower rigid and upper free boundaries (a) conducting: W = DW = T = 0 at z = 0, W = D2W = T = 0 (b) insulating: W = DW = DT = 0 at z = 0, W = D 2 W = DT = 0 at z = 1.
249
(20) (21)
at z = 1, (22)
(23)
Equations (16) and (17) with the conducting or insulating boundary conditions constitute an eigenvalue problem. The solution techniques used to find the eigenvalues for these two types of temperature boundary conditions are different, and are given below separately.
3. Solution for the Conducting Boundaries When the boundaries are conducting, no analytical solution is possible. Hence the Galerkin method is used to find the eigenvalue Ra. Accordingly, the variables are written in a series of basis functions W =
N
Ai Wi ,
(24)
Bi Ti ,
(25)
i=1
T =
N i=1
in which Ai and Bi are constants and the basis functions Wi and Ti are represented by the power series satisfying the respective boundary conditions. Substituting back Equations (24) and (25) into Equations (16) and (17); multiplying Equation (16) by Wj and Equation (17) by Tj , integrating between z = 0 and 1, and simplifying them using the boundary conditions, we obtain the following system of homogeneous algebraic equations Cj i Ai + Dj i Bi = 0,
(26)
Ej i Ai + Fj i Bi = 0.
(27)
250
A. KHALILI AND I. S. SHIVAKUMARA
The coefficients Cj i to Fj i involve the inner products of the basis functions and are given by Cj i = D 2 Wj D 2 Wi + (2a 2 + G)DWj DWi + a 2 (a 2 + G)Wj Wi + Pe (DWj D 2 Wi − a 2 DWj Wi ), (28) + Pr 2 Dj i = −Ra2 Wj Ti ,
(29)
Ej i = f (z)Tj Wi ,
(30)
Fj i = DTj DTi + a 2 Tj Ti − PeDTj Ti .
(31)
To obtain the eigenvalue for different sets of boundary conditions, the following polynomial type of basis functions are used: (i) (ii) (iii) (iv)
Wi = zi − 2zi+2 + zi+3 (free–free boundaries), Wi = zi+1 − 2zi+2 + zi+3 (rigid–rigid boundaries), Wi = 3zi+1 − 5zi+2 + 2zi+3 (lower rigid and upper free boundaries), Ti = zi+1 − zi .
The basis functions chosen above satisfy the boundary conditions automatically. Using these, the inner products were evaluated analytically rather than integrating them numerically. The coupled equations (26) and (27) can now be expressed in the form of a generalized eigenvalue problem with Rayleigh number, Ra, as the eigenvalue. The subroutine GVLRG of IMSL Library is used to find the critical Rayleigh number Rac with respect to wavenumber when other parameters are specified. 4. Solution for the Insulating Boundaries When the boundaries are insulating to temperature perturbations, a closed form solution can be obtained for the eigenvalue problem. In this case, the critical wavenumber is small, and a regular perturbation technique with wavenumber a as the perturbation parameter may be employed. Hence, the variables are expressed in a series solution of the form (W, T ) =
∞
a 2m [Wm , Tm ].
(32)
m=0
Substituting Equation (32) in Equations (16) and (17) and in the boundary conditions (24), (26) and (28), one obtains at zeroth order D 4 W0 −
Pe 3 D W0 − GD 2 W0 = 0, Pr 2
D 2 T0 − Pe DT0 = f (z)W0
(33) (34)
251
NON-DARCIAN EFFECTS
with W0 = 0 = D 2 W0 on the free, and W0 = 0 = DW0 on the rigid boundaries. The solution of Equations (33) and (34) for all the velocity boundary conditions (i.e. free–free, rigid–rigid and rigid–free) is W0 = 0 and T0 = 1. The equations at the first order of a 2 are, then D 4 W1 −
Pe 3 D W1 − GD 2 W1 = Ra, Pr 2
D 2 T1 − Pe DT1 = 1 + f (z)W1 .
(35) (36)
The general solution of Equation (35) is W1 = K1 + K2 z + K3 eaz + K4 ebz −
Ra z2 , 2G
(37)
where
Pe Pe + + 4 G, a = 2 Pr Pr 2 Pe Pe + + 4 G, b = Pr 2 Pr 2
and Ki (i = 1, 4) are constants of integration to be determined by using the boundary conditions. Applying the solvability condition to Equations (35) and (36) (i.e. zeroth order solutions are orthogonal to the first order solutions), one gets Ra W0 = 0 and [1 + f (z)W1 ]T0 = 0, where the angle brackets denote the integration from z = 0 to 1. Note that the first condition is satisfied automatically, and the second requires 1 + f (z)W1 = 0.
(38)
Using this equation, analytical expressions for critical Rayleigh number are obtained for free–free, rigid–rigid and rigid–free boundaries by determining the corresponding W1 , from Equation (37). The constants Ki (i = 1, 4) obtained for these boundary conditions are given in Appendix A. 4.1. FREE – FREE BOUNDARIES For this case, Rac is found to be 4 Pe2 G(q 2 − r 2 )2 (cosh s2 − cosh s1 ) , c1 cosh s1 + c2 cosh s2 + c3 cosh s3 + c4 sinh s1 + c5 sinh s2 + e1 (39) where r = Pe/Pr 2 + 4 G/2 , q = Pe/2 Pr 2 , s1,2 = Pe/2 ± r, s3,4 = Pe/2 ± q and the coefficients c1 –c5 and e1 are given in the Appendix. Rac =
252
A. KHALILI AND I. S. SHIVAKUMARA
As can be seen easily from Equation (39), Rac is an even function of Pe, and hence, it remains independent of the throughflow direction. In other words, the direction of throughflow does not have any influence on the stability of the system. Though a similar observation is made for conducting boundaries through numerical calculations, it should be noted that this behaviour could not be seen explicitly. As Da−1 → ∞, Equation (39) takes the form Rac ∼
2 Pe2 G , [Pe coth(Pe/2) − 2]
(40)
which corresponds to the Darcy–Forchheimer case. Further, it should be noted that if cf = 0, the above result reduces to that of Nield (1987), for which Rac is independent of Pr. In addition, Rac ∼ 12G{1 + Pe2 /60} as Pe → 0. Note that, in the absence of throughflow (Pe = 0), the inertia and hence, the Prandtl number effects do not appear, and Rac ∼ 12 Da−1 , which is the known exact value. Also, Equation (39) reduces to a 12 D as Pe → 0, (41) 3 a − 12 D a + 24 tanh(D a/2)] [D a = Da−1 / . From Equation (41), it follows that Rac ∼ 120 as Da−1 → where D 0, = 1, and Rac ∼ 12 Da−1 as Da−1 → ∞, which are the known exact values. 5
Rac ∼
4.2. RIGID – RIGID BOUNDARIES In this case Rac =
2 Pe G[−2r(sinh s4 + sinh s3 − sinh s1 − sinh s2 ) + e1 ] , c1 cosh s1 + c2 cosh s2 + c3 cosh s3 + c4 sinh s1 + e2
(42)
where r, q and s1 –s4 have their pre-defined meaning, while other parameters are given in Appendix A. It is seen that, also in this case, Rac is an even function of Pe because of symmetric boundary conditions. For the Darcy–Forchheimer case (Da−1 → ∞), the above equation reduces to Rac ∼
2 Pe2 G , [Pe coth(Pe/2) − 2]
(43)
which is identical to Equation (40). Further, as Pe → 0, Equation (42) becomes a4 (2 + D a sinh D a − 2 cosh D a) 12 D . (44) Rac ∼ a2 ) + Da(24 + D a2 ) sinh D a − 8(D a2 + 3) cosh D a] [4(6 − D −1 −1 This √ in turn reduces to Rac ∼ 720 as Da → 0, = 1 and Rac ∼ 12 Da + 72 Da−1 as Da−1 → ∞, = 1, which represent the well-known results.
253
NON-DARCIAN EFFECTS
4.3. LOWER BOUNDARY RIGID AND UPPER BOUNDARY FREE For this case, Rac is found to be Rac =
c1 ePe+q+r + c2 ePe+q
2 Pe2 G(ePe − 1)H . + c3 ePe+r + c4 ePe + c5 eq + c6 er + c7 eq+r + c8 (45)
The coefficients c1 –c8 are given in Appendix A, and other parameters are defined earlier. In contrast to the previous two types of boundary conditions, here, Rac is not an even function of Pe. As a consequence, the direction of throughflow also has a say on the stability of the system. However, Rac for the Darcy–Forchheimer case turns out to be the same as that of free–free and rigid–rigid boundaries. This implies that in the absence of viscous shear, Rac is independent of the boundary conditions considered. Further, Rac ∼
a5 (sinh D a − D a cosh D a) 12 D , 2 4 a(D a − 6) sinh D a + (24 − D a ) cosh D a + 12(D a2 − 2)] [4 D (46)
−1 as the throughflow √ vanishes. Note that Rac ∼ 320 as Da → 0, = 1 and Rac ∼ −1 −1 −1 12 Da + 36 Da as Da → ∞, = 1, which are the known standard results.
5. Results and Discussion The effects of macroscopic and microscopic inertial forces, along with the nonlinear basic temperature gradient arising due to throughflow are investigated on convective instabilities in a porous layer. The Lapwood–Forchheimer–Brinkmanextended Darcy model is used to describe the fluid motion. The boundaries are taken to be either rigid–permeable or stress-free and either conducting or insulating to temperature perturbations. The Galerkin method is employed to find the critical Rayleigh number, Rac for free–free, rigid–rigid and rigid–free conducting boundaries. Prior to applying the method to the complete problem, the method of solution was verified by comparing our results with the results known for the case of no throughflow and also without a porous medium by setting Pe → 0, Da−1 → 0, = 1, = 1 and cf = 0. Our results of the critical wavenumber and Rayleigh number (ac , Rac ) = (2.221, 657.512), (3.116, 1707.761) and (2.682, 1100.656) for free– free, rigid–rigid and rigid–free boundaries, respectively, with i = 5 in Equations (24) and (25) agree exactly with those given by Chandrasekhar (1961). However, the presence of throughflow brings in Prandtl number effect through inertia terms, and the coefficient matrix formed from Equations (26) and (27) become asymmetric. As a consequence, it is observed that the convergence of results crucially
254
A. KHALILI AND I. S. SHIVAKUMARA
depends on the value of |Pe|/Pr. For the parametric values considered, it is found that the method gives convergent results for i = 8. In this study, Da−1 of 10 and 50, of 1 and 5, cf of 0.1 and 0.55 and Pr of 0.8 and 7 for a fixed value of 2 = 0.5 are chosen for calculations. The critical Rayleigh numbers were calculated for all types of boundary conditions. It was found that the results for free–free boundaries differ from those of rigid–rigid boundaries quantitatively, and hence only the latter case is considered
Figure 2. Variation of Rac versus Pe for rigid–rigid conducting boundaries and different models; (a): cf = 0.1, Pr = 0.8, (b): = 5, Pr = 0.8, (c): cf = 0.55, Da−1 = 50.
NON-DARCIAN EFFECTS
255
for discussion to avoid repetition. Figure 2(a)–2(c) show the variation of Rac as a function of Pe for rigid–rigid conducting boundaries. In these figures, the relative importance of inertial resistance on the onset of convection is demonstrated separately by retaining either only Lapwood term (L) or only Forchheimer term (F) or both (L + F), in the momentum equation. It can be seen that the effect of inertia on the stability of the system gains significance with an increase in the throughflow strength. In general, the presence of both L and F makes the system more stable as compared to their presence in isolation. However, among L and F, the system is less stable with the presence of only F when cf = 0.1 (see Figure 2(a)), while an opposite or mixed type of behaviour is noticed when cf = 0.55 (see Figure 2(b)). It may be noted that the deviation in the Rac values between L and F is more for Da−1 = 10 than Da−1 = 50. Whereas this deviation is found to be not so significant between L and L + F. Further, an increase in Da−1 , and cf increases Rac and thus makes the system more stable. Figure 2(c) shows the effect of Pr on the stability of the system for = 1 and 5 with cf = 0.55 and Da−1 = 50. It is evident that a decrease in Pr enhances the stability of the system, and the effect of L or F or L + F on Rac could be discerned only for Pr = 0.8. In all the cases, Rac values increase with an increase in Pe and they are independent of throughflow direction. This may be due to the fact that throughflow is to confine significant thermal gradients to a thermal boundary layer at the boundary towards which the throughflow is directed. The effective length scale is thus smaller than the thickness of the porous layer d. Hence the Rayleigh number (being proportional to cube of the porous layer thickness) will be much less than the actual value of Rac . Therefore higher values of Rac are needed for the onset of convection when the throughflow strength increases. For rigid–free boundaries, an entirely different type of behaviour was noticed and the same is depicted in Figures. 3(a)–(e). As can be seen, the direction of throughflow changes the critical Rayleigh number. Moreover, a small amount of throughflow in one direction produces some destabilization depending on the choice of physical parameters and also the models employed. For Da−1 = 10, Pr = 0.8, = 1 and cf = 0.55, it is noted that the downward (Pe < 0) and upward (Pe > 0) flows destabilize the system initially for L and F, respectively, while L + F always stabilizes the system (see Figure 3(a)). The Rac curve for F falls below L and L + F for upflow, while a dual behaviour could be seen for downflow. However, when Pr, and Da−1 are increased correspondingly to 7, 5 and 50, it is observed that the system initially becomes unstable for only upward flow irrespective of the model considered (see Figures 3(b)–(e)). In all these cases, the Rac curves for F always lie above L and L + F for downflow. Furthermore, the range of Pe up to which the system destabilizes is increased with growing Pr, and reducing Da−1 . The destabilization may be due to the distortion of the steady state basic temperature distribution by the throughflow. A measure of this is given by the term f (z)W T and this can be interpreted as a rate of transfer of energy into the disturbance by interaction of the perturbation convective motion with basic temperature gradient. The maximum temperature occurs at a place where the perturbed vertical
256
A. KHALILI AND I. S. SHIVAKUMARA
Figure 3. Variation of Rac versus Pe for rigid–free conducting boundaries and different models with cf = 0.55; (a) Pr = 0.8, = 1, Da−1 = 10, (b) Pr = 7, = 1, Da−1 = 10, (c) Pr = 0.8, = 5, Da−1 = 10, (d) Pr = 7, = 5, Da−1 = 10, (e) Pr = 7, = 5, Da−1 = 50.
NON-DARCIAN EFFECTS
257
velocity is large, and this leads to an increase in energy supply for destabilization. The destabilization may also be due to other mechanisms viz., the momentum and thermal energy transports. When the boundaries are insulating, only the results for rigid–free case is presented because it provides more insight into the physics of the problem considered. Figure 4(a) and (b) shows the variation of critical Rayleigh numbers with Pe for
Figure 3. (continued).
258
A. KHALILI AND I. S. SHIVAKUMARA
Figure 3. (continued).
Da−1 = 10 and 50, respectively, for fixed Pr = 7, = 5 and cf = 0.55. It should be noted that the results are similar to those of conducting boundaries (see Figure 3(d) and (e)), except for a quantitative change. As expected, the critical
Figure 4. Variation of Rac Pe for rigid–free insulating boundaries and different models with cf = 0.55, Pr = 7, = 5; (a) Da−1 = 10, (b) Da−1 = 50.
NON-DARCIAN EFFECTS
259
Figure 4. (continued).
Rayleigh numbers for the conducting boundaries are higher than those of insulating boundaries. Figure 5(a) and (b) illustrates the variation of the critical wavenumber ac as a function of Pe for rigid–rigid and rigid–free conducting boundaries, respectively.
Figure 5. Variation of ac Pe for different models with cf = 0.55, Pr = 0.8 and Da−1 = 50; (a) rigid–rigid conducting boundaries, (b) rigid–free conducting boundaries.
260
A. KHALILI AND I. S. SHIVAKUMARA
Figure 5. (continued).
The results presented here are for = 1 and 5 with cf = 0.55, Pr = 0.8 and Da−1 = 50. From these figures, it is obvious that an increase in and the throughflow strength reduces the size of the convection cells. Acknowledgements One of the authors (I.S.S.) is pleased to express his thanks to German Research Foundation (DFG) for the financial support (grant 446 IND 17/1/00) and Max Planck Institute for Marine Microbiology, Bremen, Germany for the hospitality. Appendix A A.1. FREE – FREE BOUNDARIES The boundary conditions are: W1 = D 2 W1 = 0 at z = 0, 1. The constants of integration appearing in Equation (37) are Ra[a 2 (ea − 1) − b2 (eb − 1)] , HG Ra[a 2 b2 (eb − ea ) + 2(a 2 − b2 )(1 − ea )(1 − eb )] , K2 = 2H G Ra b2 (eb − 1) , K3 = HG
K1 =
NON-DARCIAN EFFECTS
K4 =
−Ra a 2 (ea − 1) HG
with H = a 2 b2 (eb − ea ). The coefficients appearing in Equation (39) are given by c1 = c2 = c3 = c4 = c5 = e1 =
4 Pe2 (q 2 − r 2 )(q − r) + 4(q 2 − r 2 )2 , Pe + q + r −4 Pe2 (q 2 − r 2 )(q + r) − 4(q 2 − r 2 )2 , Pe + q − r [q 2 + r 2 + 2q(Pe + q)] 16qr Pe2 − 8r Pe3 , (Pe + q)2 − r 2 −2 Pe(q 2 − r 2 )2 + 16 Pe qr, 2 Pe(q 2 + r 2 )2 + 16 Pe qr, −16 Pe qr(sinh s3 + sinh s4 ).
A.2. RIGID – RIGID BOUNDARIES The boundary conditions are: W1 = DW1 = 0 at z = 0, 1. Hence, Ra[a(1 + ea ) − b(1 + eb ) + 2(eb − ea )] , HG Ra[2a(eb − 1) − 2b(ea − 1) + ab(eb − ea )] , K2 = HG Ra[b(1 + eb ) − 2(eb − 1)] , K3 = HG Ra[2(ea − 1) − a(1 + ea )] K4 = HG
K1 =
with H = 2(b − a)(1 − ea − eb + ea+b ) + 2ab(ea − eb ). The coefficients appearing in Equation (42) are given by 2 Pe(q + r)2 − (q 2 − r 2 ), Pe + q + r 2 Pe(q − r)2 + (q 2 − r 2 ), c2 = − Pe + q − r 4 Pe r − 4 Pe r, c3 = (Pe + q)2 − r 2
c1 =
261
262
A. KHALILI AND I. S. SHIVAKUMARA
c4 = c5 = c6 = c7 = e1 = e2 =
Pe(q + r)2 (q − r) − 4r, Pe + q + r Pe(q − r)2 (q + r) − 4r, Pe + q − r Pe2 (q + r)(Pe + 2r) + 4r, Pe + q + r Pe3 (q + r) + 4r, − Pe + q − r (q 2 − r 2 )(cosh s1 − cosh s2 ), c5 sinh s2 + c6 sinh s3 + c7 sinh s4 .
A.3. LOWER BOUNDARY RIGID AND UPPER BOUNDARY FREE The boundary conditions are: W1 = DW1 = 0 at z = 0,
W1 = D 2 W1 = 0 at z = 1.
Hence, Ra[a 2 ea − b2 eb + 2(eb − b) − 2(ea − a)] , 2H G Ra[a 2 b ea − ab2 eb + 2a(eb − 1) − 2b(ea − 1)] , K2 = 2H G Ra[b2 eb − 2(eb − b − 1)] , K3 = 2H G Ra[−a 2 ea + 2(ea − a − 1)] K4 = 2H G
K1 =
with H = b2 eb (ea − a − 1) − a 2 ea (eb − b − 1). The coefficients occurring in Equation (45) are given by Pe3 (q 2 − 2) + (q 2 − r 2 )(Pe2 + 2 Pe + 2), Pe + q −2 Pe2 q(1 + r) + Pe q 2 (r + 2) + 2 Pe r − 2(1 + r)q 2 , c2 = Pe + q 2 Pe2 r(1 + q) − Pe r 2 (q + 2) + 2 Pe q + 2(1 + q)r 2 , c3 = Pe + r c4 = −2 Pe(q + r), c1 =
NON-DARCIAN EFFECTS
263
Pe r 2 (q 2 − 2) + 2(1 + r)q 2 , Pe + r −Pe q 2 (r 2 − 2) − 2(1 + q)r 2 , c6 = Pe + q c7 = 2(r 2 − q 2 ), 2 Pe3 (q − r)(q + r + Pe + 1) − 2 Pe(Pe + 1)(q − r). c8 = (Pe + q)(Pe + r) c5 =
References Chandrasekhar, S.: 1961, Hydrodynamic and Hydromagnetic Stability, Dover, New York. Dando, P. R., Aliani, S., Arab, H., Bianchi, C. N., Brehmer, M., Cocito, S., Fowler, S. W., Gundersen, J., Hooper, L. E., Kölbl, R., Küver, J., Linke, P., Makropoulos, K. C., Meloni, R., Miqueo, J.-C., Morri, C., Müller, S., Robinson, C., Schlesner, H., Sievert, S., Stöhr, R., Stüben, D., Thomm, M., Varnavas, S. P. and Ziebis, W.: 2000, Hydrothermal studies in the aegean sea, Phys. Chem. Earth 25, 1–8. Givler, R. C. and Altobelli, S. A.: 1994, A determination of the effective viscosity for the Brinkman– Forchheimer flow model, J. Fluid Mech. 258, 355–370. Holm, N. G. (ed.): 1992, Marine Hydrothermal Systems and the Origin of Life, Kluwer Academic Publishers, Dordrecht. Homsy, G. M. and Sherwood, A. E.: 1976, Convective instabilities in porous media with through flow, AIChE J. 22, 168–174. Jones, M. C. and Persichetti: 1986, Convective instability in packed beds with throughflow, AIChE J. 32, 1555–1557. Khalili, A. and Shivakumara, I. S.: 1998, Onset of convection in a porous layer with net through-flow and internal heat generation, Phys. Fluids 10, 315–317. Nield, D. A.: 1987, Convective instabilities in porous media with throughflow, AIChE J. 33, 1222–1224. Nield, D. A. and Bejan, A.: 1992, Convection in Porous Media, Springer, Berlin. Vafai, K. and Kim, S. J.: 1990, Analysis of surface enhancement by a porous substrate, ASME J. Heat Tran. 112, 700–706. Vafai, K. and Tien, C. L.: 1981, Boundary and inertia effects on flow and heat transfer in porous media, Int. J. Heat Mass Tran. 24, 195–203. Wooding, R. A.: 1960, Rayleigh instability of a thermal boundary layer in a flow through a porous medium at large Reynolds number or Péclet number, J. Fluid Mech. 9, 183–192.
