Engineering
Vol.11 No.07(2019), Article ID:93606,23 pages
10.4236/eng.2019.117026

Stability of Geothermal Convection in Anisotropic River Beds

Gérard Degan1,2, Julien Yovogan1,2, Latif Fagbémi1, Zineddine Allou3

1Laboratoire d’Energétique et de Mécanique Appliquées (LEMA), Ecole Polytechnique d’Abomey-Calavi, Cotonou, Bénin

2Université Nationale des Sciences Technologies, Ingénierie et Mathématiques (UNSTIM), Abomey, Bénin

3Mechanical Engineering Department, Ecole Polytechnique, University of Montréal, Montréal, Canada

Copyright © 2019 by author(s) and Scientific Research Publishing Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

Received: February 19, 2019; Accepted: July 7, 2019; Published: July 10, 2019

ABSTRACT

The onset of thermal convection, due to heating from below in a system consisting of a fluid layer overlying a porous layer with anisotropic permeability and thermal diffusivity, is investigated analytically. The porous medium is both anisotropic in permeability whose principal axes are oriented in a direction that is oblique to the gravity vector and in thermal conductivity with principal directions coincident with the coordinate axes. The Beavers-Joseph condition is applied at the interface between the two layers. Based on parallel flow approximation theory, a linear stability analysis is conducted to study the geothermal river beds system and documented the effects of the physical parameters describing the problem. The critical Rayleigh numbers for both the fluid and porous layers corresponding, to the onset of convection arising from sudden heating and cooling at the boundaries are also predicted. The results obtained are in agreement with those found in the past for particular isotropic and anisotropic cases and for limiting cases concerning pure porous media and for pure fluid layer. It has demonstrated that the effects of anisotropic parameters are highly significant.

Keywords:

River Beds, Critical Rayleigh Numbers, Isotropic, Anisotropic

1. Introduction

Natural convection in composite fluid and porous layers heated from below can be encountered in many engineering and environmental problems. Applications include solidification of castings, aerosol production, groundwater pollution, thermal insulation, geophysical systems, etc. Therefore a considerable amount of investigations on this subject has been performed in the last decades [1].

Most investigations of buoyancy-driven flows in two-dimensional fluid/porous systems are concerned with the case of a shallow rectangular cavity partially filled with a porous medium. The interaction that occurs at the interface between a fluid and a porous layer has been formulated in the past according to different approaches. Somerton and Catton [2] studied the stability of a system consisting of a volumetrically heated porous bed overlaid with a fluid layer, heated or cooled isothermally from below. The case of a rectangular cavity, divided into fluid and porous regions, has been investigated numerically by Nishimura et al. [3] using a Brinkman model. The boundary conditions at the interface between fluid and porous medium were written explicitly in stream function-vorticity form. A good agreement was observed between their numerical results and experimental data. The case of a porous bed under a fluid, in a shallow cavity heated from the bottom by a constant heat flux, has been investigated analytically by Vasseur et al. [4]. Closed form solutions were obtained using NavierStokes and Brinkmans equations in fluid and porous regions, respectively. Chen et al. [5] used the BJ condition to predict the onset of motion in a system consisting of a fluid layer overlying an anisotropic porous medium. The onset of surface-tension driven convection has been investigated by Shivakumara et al. [6] in a two layer-system. Both BeaversJoseph and the Jones conditions [7] were applied at the contact surface between the fluid saturated porous medium and the adjacent bulk fluid. A negligible difference in the critical Marangoni numbers was obtained whether BJ or Jones conditions are being used. A study of buoyancy-driven flow in a confined fluid overlying a porous layer has been investigated by Valenca-Lopez and Ochoa-Tapia [8], using two different models. The first approach considered the NavierStokes equation and Darcys law coupled with the BJ interfacial boundary condition. The second approach used the Brinkman extended Darcys law together with the continuity of shear and velocity at the interface. Significant differences between the overall Nusselt numbers were found when the Rayleigh and Darcy numbers were large enough. Five different types of interfacial conditions between a porous medium and a fluid layer have been analyzed by Alazmi and Vafai [9]. It was found that the variances within different models, for most practical applications, have a negligible effect on the results. Recently, the influence of the interfacial jump boundary condition on the onset of the convection in superposed fluid and porous layers has been studied by Hirata et al. [10] [11]. The effective jump coefficient was found to strongly influence the marginal stability curves. However the results indicate that the inclusion of the Brinkman term, to model the porous medium, plays a secondary role on the stability results. More recently a linear stability analysis of the onset of thermosolutal convection in horizontal superposed fluid and porous layers was performed by Hirata et al. [12] using the one-domain approach. It was demonstrated that for positive thermal Rayleigh numbers, the convective flow occurs both in the fluid and porous regions. However, for negative thermal Rayleigh numbers the onset of motion is characterized by a multi-cellular flow in the fluid region. The same problem was reconsidered by Hirata et al. [13] by using both the one-domain approach and the two-domain formulation. A very good agreement was obtained upon comparing the results predicted by the two models. An analytical study of natural convection in a horizontal shallow cavity filled with a fluid layer overlying a porous medium, saturated with the same fluid, has been investigated by Alloui et al. [14] using a Darcy model. The critical Rayleigh numbers for the onset of convection in a composite system are predicted explicitly, by the present model, in terms of the governing parameters of the problem. For finite-amplitude convection, useful expressions have been obtained for velocity, and temperature distributions in the core of the enclosure. At the interface between the fluid and the porous layers the effect of the slip coefficient is more pronounced on the velocity profiles in the fluid layer than on those in the porous layer. The influence of the governing parameter on the Nusselt numbers are predicted and discussed. Finally, it is noted that the scope of this study is limited by the assumption of a parallel flow, and the model is thus valid in the limit of a shallow enclosure ( A 1 ) . Despite this restriction, the analytical model provides useful results which can be taken as a starting point for more detailed numerical computations.

The aim of the present study is to study analytically natural convection in a cavity consisting of a fluid layer over a saturated porous layer anisotropic in permeability. The system is heated from the bottom by a constant heat flux. At the interface of the two layers the Beavers-Joseph boundary conditions are applied. In this investigation, we use the Navier-Stokes equations for the fluid layer and the Darcy law with the presence of gravitational field for the porous layer A parallel flow approximation is used, which enables the temperature and velocity fields in the core region of the system to be determined in closed form. The critical Rayleigh for both fluid layer and porous layer corresponding the onset of convection are predicted.

2. Conceptual Model

The physical model illustrating the problem under different considerations is shown in Figure 1. A porous layer of thickness h p (zone 2) underlying a fluid layer of thickness h f (zone 1) is considered. Between the first and the second zone, there is an interface (called the nominal surface). The layers are horizontal and extend infinitely in the horizontal direction. The top of the fluid layer is free and the bottom of the porous layer is bounded by rigid walls with uniform heating. It is assumed that the flow is laminar, incompressible, and two-dimensional. The physical properties of the fluid are assumed constant, except for the density in the buoyancy term in the momentum equations (Boussinesq approximation). The porous medium is considered homogeneous and anisotropic and saturated with a fluid which is in local thermodynamic equilibrium with the solid matrix. A Cartesian coordinate system is chosen with the origin at the interface between the porous and fluid layers, the y’-axis vertically upward, and the x’-axis horizontal.

Figure 1. Schematic diagram of the physical model and coordinate system.

The porous medium is anisotropic, the permeabilities along the two principal axes of the porous matrix are denoted by K 1 and K 2 . The anisotropy of the porous layer is characterized by the permeability ratio K * = K 1 / K 2 and the orientation angle φ , defined as the angle between the horizontal direction and the principal axis with the permeability K 2 and the thermal diffusivity ratio ε T , defined as ratio of thermal diffusivity of the fluid layer α f to the thermal diffusivity of the porous layer α p .

3. Mathematical Formulation

The equations governing the conservation of mass, momentum, energy and the Boussinesq approximation (see, Alloui et al. [14] Yovogan and Degan [15] Yovogan et al. [16] ) can be written in each Zone as follows:

• Zone 1 (fluid layer):

Equation governing the conservation of mass

V f = 0 , (1)

Equation governing the conservation of momentum (Navier-Stokes model with the presence of gravitational field).

ϱ 0 ( V f ) V f = P f + μ 2 V f + ϱ g , (2)

Equation governing the conservation of energy

( ϱ C p ) f D T f D t = k f 2 T f , (3)

Equation governing the Boussinesq approximation

ϱ = ϱ 0 [ 1 b ( T f T 0 ) ] . (4)

• Zone 2 (porous layer):

Equation governing the conservation of mass

V p = 0 , (5)

Equation governing the conservation of momentum (Darcy law with the presence of gravitational field).

V p = K ¯ ¯ m { P p + ϱ g } , (6)

Equation governing the conservation of energy

( ϱ c p ) p T p t + ( ϱ c p ) f ( V p T p ) = k p 2 T p , (7)

Equation governing the Boussinesq approximation

ϱ = ϱ 0 [ 1 β ( T p T 0 ) ] . (8)

In these equations, V denotes the velocity vector, p the pressure and T the temperature. The subscript f denotes the fluid layer, p the porous medium. Moreover, μ the dynamic viscosity, ϱ the density, β the thermal-expansion coefficient, c p the specific heat of the fluid, k the thermal conductivity, α f = k / ( ϱ c p ) f the thermal diffusivity of the fluid layer and α p = k / ( ϱ c p ) p the thermal diffusivity of the porous layer, where ( ϱ c p ) is the volumetric heat capacity of the fluid. In Equation (6), the symmetrical second-order permeability tensor K ¯ ¯ is defined as

K ¯ ¯ = [ K 1 sin 2 φ + K 2 cos 2 φ ( K 2 K 1 ) sin φ cos φ ( K 2 K 1 ) sin φ cos φ K 2 sin 2 φ + K 1 cos 2 φ ] . (9)

The appropriate boundary conditions prevailing on the the lower impermeable boundary and the upper free surface of the channel are:

y = h p : v p = 0 , u p = 0 , T p y = q k p , (10)

y = h f : v f = 0 , d u f d y = 0 , T f y = q k f , (11)

At the interface of the two layers ( y = 0 ) , the conventional no-slip velocity boundary condition can be assumed to be valid even at the impermeable walls. However, Beavers and Joseph [17] postulated the existence of a streamwise slip velocity at the permeable bounding surface. So the dynamic conditions and the continuity of temperature, of the heat flux on the nominal surface are formulated as

y = 0 : v p = v f , u f y = β 1 K 1 ( u f u p ) , (12)

y = 0 : T p = T f , k f T f y | y = 0 = k p T p y | y = 0 . (13)

According to Rudraiah and Veerabhadraiah [18] [19], the parameter β 1 denotes a constant depending on the material property of the porous medium, which have can be determined only experimentally.

To render the equations non-dimensional, the characteristic length is chosen to be the total height, h = h p + h f , of the fluid and porous layers and the characristic velocity to be α / h . Different scales for temperature are used; they are ( T f T 0 ) / Δ T and ( T p T 0 ) / Δ T , for the fluid and porous layers, respectively. Then, the dimensionless formulation ( ψ , T ) of governing equations become:

• (for the fluid layer):

J ( ψ f , 2 ψ f ) = P r R a T f x + P r 4 ψ f , (14)

J ( ψ f , T f ) = 2 T f . (15)

• (and for the porous layer):

a ( 2 ψ p x 2 ) b ( 2 ψ p y 2 ) + 2 c ( 2 ψ p x y ) = R T p x , (16)

J ( y p , T p ) = ε T 2 T p . (17)

where P r = μ C p / k is the Prandtl number, μ the dynamic viscosity of the fluid. J ( f , g ) = f x g y f y g x the usual Jacobian operator, R a = ( g Δ T β h 3 K 1 ) / ν α p the Rayleigh number in the fluid layer, R = ( g Δ T β h K 1 ) / ν α p = D a R a the Rayleigh number in the porous layer, D a = K * / h 2 the Darcy number and ψ is the usual stream function defined as:

u = ψ y , v = ψ x . (18)

such that the mass conservation is satisfied. the constants a, b, and c are defined as.

{ a = K * sin 2 φ + cos 2 φ , b = sin 2 φ + K * cos 2 φ , c = ( 1 K * ) sin φ cos φ . (19)

The dimensionless boundary conditions on the lower impermeable boundary of the porous medium and on the upper free surface of the fluid layer are

y = η : ψ p = 0 , T p y = 1 / ε T , (20)

y = η ¯ : ψ f = d 2 ψ f d y 2 = 0 , T f y = 1. (21)

At the interface of the two layers (at y = 0 ), we have:

T p = T f , (22)

d 2 ψ f d y 2 = β 1 D a ( d ψ f d y d ψ p d y ) , (23)

T f y | y = 0 = ϵ T T p y | y = 0 . (24)

4. Parallel Flow Hypothesis

Assuming that when the flow is fully developed in the system, the axial (x-direction) velocity depends on the transverse coordinate y (i.e. u f = u f ( y ) for the fluid layer and u p = u p ( y ) for the porous layer), and then from the continuity equation, the transverse velocity component must be zero (i.e. v = 0 and v p = 0 ). The temperature field, in the central part, can be divided into the sum of a linear dependence on x and an unknown function of y. Thus, it is assumed that

ψ f ( y ) = R × C × H f ( y ) , ψ p ( y ) = R × C × H p ( y ) , (25)

T f ( x , y ) = C × x + θ f ( y ) , T p ( x , y ) = C × x + θ p ( y ) . (26)

where C is the dimensionless horizontal temperature gradient in the horizontal direction. C is the same in the two layers. Similar approximations have been used in the past by Cormack et al. [20], Sparrow et al. [21], Vasseur et al. [22], Sen et al. [23] [24] and Bahloul et al. [25], among others.

Substituting Equations (25)-(26) into Equations (14)-(17), the governing equations for the fluid layer can be reduced to the following ordinary differential equations:

d 4 ψ f d y 4 = R a × C , (27)

d 2 θ f d y 2 = C d ψ f d y . (28)

while for the porous layer, the ordinary differential equations obtained are given by

b d 2 ψ p d y 2 = R × C , (29)

d 2 θ p d y 2 = C ε T d ψ p d y . (30)

5. Analytical Solution

5.1. Velocity and Temperature Distrubution

Equations (27)-(30) can be solved, sobjected to boundary conditions, Equations. (20)-(24). The resulting expressions for the velocity, stream function and temperature fields for the fluid layer are given by

u f ( y ) = ψ 0 ( 1 6 y 3 + 1 2 A 1 y 2 + A 2 y + A 3 ) , (31)

ψ f ( y ) = ψ 0 ( 1 24 y 4 + 1 6 A 1 y 3 + 1 2 A 2 y 2 + A 3 y + A 4 ) , (32)

T f ( x,y )=Cx+C ψ 0 ( 1 120 y 5 + 1 24 A 1 y 4 + 1 6 A 2 y 3 + 1 2 A 3 y 2 + A 4 y )+ A 5 y. (33)

while those in the porous layer are given by

u p ( y ) = γ ψ 0 ( y + B 1 ) , (34)

ψ p ( y ) = γ ψ 0 ( 1 2 y 2 + B 1 y + B 2 ) , (35)

T p ( x , y ) = C x + γ C y 0 ε T ( 1 6 y 3 + 1 2 B 1 y 2 + B 2 y ) + B 3 y . (36)

where

γ = D a / b , (37)

A 1 = ( 5 η ¯ 4 ) / 24 + ( D a * η ¯ 3 ) / 2 ( γ η 2 ) / 2 η ¯ 3 / 3 D a * η ¯ 2 γ , (38)

A 2 = ( η ¯ 5 ) / 24 + ( γ η ¯ 3 ) / 2 + ( γ η η ¯ ) / 2 η ¯ 3 / 3 D a * η ¯ 2 γ , (39)

A 3 = ( D a * η ¯ 5 ) / 24 ( 5 γ η ¯ 4 ) / 24 + ( D a * η η ¯ ) / 2 + ( γ η ) 2 / 2 η ¯ 3 / 3 D a * η ¯ 2 γ , (40)

A 4 = γ [ ( D a * η η ¯ 2 ) / 2 ( η η ¯ 3 ) / 4 ( γ η 2 η ¯ ) / 2 + ( η 2 η ¯ 3 ) / 12 + ( η η ¯ 4 ) / 24 ] η ¯ 3 / 3 D a * η ¯ 2 γ . (41)

{ A 5 = 1 , B 1 = A 1 , B 2 = A 1 η η 2 / 2 , B 3 = 1 ε T , D a * = D a β 1 , ψ 0 = R a × C . (42)

5.2. Critical Rayleigh Number (Onset of Motion in the System)

The parallel flow approximation is only applicable in the core of the layers. Flows in the end regions are much more complicated and cannot be approximated in such a simple manner. For this reason, the thermal boundary condition in the x-direction cannot be reproduced exactly with this approximation. We can, however, impose an equivalent energy flux condition in that direction by writing as follows

η 0 ( d ψ p d y T p T p x ) x = 0 d y + 0 η ¯ ( d ψ f d y T f T f x ) x = 0 d y = 0. (43)

The integrands are a sum of convective heat fluxes in the fluid and porous medium, respectively. This is derived from the condition of uniform heat flux at the boundaries. Substituting Equations (32)-(33) and (35)-(36) into Equation (43) and integrating yields, after some straightforward but laborious algebra, an expression of the form:

[ E 1 + E 3 ] R a 2 C 3 + [ ( ε T η + η ¯ ) ( E 2 + E 4 ) R a ] C = 0. (44)

where

E 1 = η ¯ 9 5184 + A 1 η ¯ 8 576 + ( A 1 2 36 + A 2 24 ) η ¯ 7 7 + ( A 1 A 2 6 + A 3 12 ) η ¯ 6 6 + ( A 1 A 3 3 + A 4 12 + A 2 2 4 ) η ¯ 5 5 + ( A 1 A 4 3 + A 2 A 3 ) η ¯ 4 4 + ( A 3 2 + A 2 A 4 ) η ¯ 3 3 + A 3 A 4 η ¯ 2 + A 4 2 η ¯ , (45)

E 2 = η ¯ 5 120 + A 1 η ¯ 4 24 + A 2 η ¯ 3 6 + A 3 η ¯ 2 2 + A 4 η ¯ , (46)

E 3 = γ 2 ε T [ η 5 20 B 1 η 4 4 + ( B 1 2 + B 2 ) η 3 3 B 1 B 2 η 2 + B 2 2 η ] , (47)

E 4 = γ ε T [ η 3 6 + B 1 η 2 2 B 2 η ] . (48)

From Equation (44) it is seen that, apart of the trivial case C = 0 corresponding to the rest state, the value of C for finite amplitude convection is given by

C = ± ( E 2 + E 4 ) R a ( ε T η + η ¯ ) [ E 1 + E 3 ] R a 2 . (49)

Equation (49) indicates that when R a = ( ε T η + η ¯ ) / ( E 2 + E 4 ) , C = 0 is the real value of C and there is no convection. On the other hand, when R a > ( ε T η + η ¯ ) / ( E 2 + E 4 ) a convective cell rotating clockwise ( C < 0 ) or counter-clockwise ( C > 0 ) , bifurcates from the rest state. Physically, this follows from the fact that the convective cells resulting from the onset of Rayleigh-Benard convection, can rotate indifferently in a direction or in the other. The velocity, the stream function and temperature distributions can then be evaluated from Equations (31)-(36), once the value of C has been evaluated from Equation (49).

The critical Rayleigh number, R a c , for the onset of convection, is obtained from Equation (49), for the condition C = 0, as:

R a c = ε T η + η ¯ E 2 + E 4 . (50)

5.3. Comparison with the Results Predicted in the Past

The marginal stability of the composite system (anisotropic river beds) considered in this investigation is given by Equation (50). We can check this formula against known results for the following limiting cases.

5.3.1. Case of a Pure Anisotropic Porous Layer ( η = 1 )

Letting η = 1 and ε T = 1 , it is readily found from Equation (50) that

R c = D a R a c = 12 b , (51)

were the letter b takes into account the parameters of anisotropy, K * and φ . This result (Equation (51)) is in agreement the with the result predicted in the past by Nield [26], P. Vasseur et al. [22], Degan et al. [27], on the basis of the linear stability theory for the case of a horizontal Darcy porous layer heated from the bottom by a constant heat flux.

5.3.2. Case of a Pure Isotropic Porous Layer

When the permeability is the same in all directions (i.e. for an isotropic porous layer), we have: K 1 = K 2 and b = K * cos 2 ( φ ) + sin 2 ( φ ) = 1 . Furthermore, it is easily found that, Equation (51) yields

R c = D a R a c = 12 , (52)

which is reported in the past by Alloui et al. [14] on the basis of the linear stability theory for the case of a horizontal Darcy porous layer heated from the bottom by a constant heat flux.

5.3.3. Case of a Horizontal Fluid Layer ( η = 0 )

For the case η = 0 , corresponding to a layer of fluid with an upper solid boundary and a lower porous lining, it is found from the present theory that we have

R a c = 72000 ( 1 + 3 D a * + 3 D a ) 100 + 675 D a * + 1800 D a . (53)

The above result, in the limit D a 0 , reduces to

R a c = 720. (54)

which corresponds to the critical Rayleigh number for a layer of fluid bounded by two solid walls as predicted by Alloui et al. [14] and Sparrow et al. [21]. Furthermore, it is easily found that, when β 1 0 , Equation (52) yields

R a c = 320. (55)

which is also a value reported in the past by Alloui et al. [14] and Sparrow et al. [21] in the case of a fluid layer with a solid horizontal lower boundary and a free upper surface.

Finally, in the limit of a high Darcy number D a , it is found from Equation (52) that

R a c = 120. (56)

Which is in agreement the with the result predicted in the past by Alloui et al. [14] and Sparrow et al. [21].

5.3.4. Case of a Liquid Layer Topping a Solid Slab

In the limit D a 0 , it can be demonstrated from Equation (50) that the critical Rayleigh number is given by the following simplified expression:

R a c = 320 ( η ε T + 1 η ( 1 η ) 5 ) . (57)

which yields the marginal stability condition for a system consisting of a liquid layer over a solid slab.

5.4. Nusselt Number (Heat Transfer Rate)

Since the temperature of each thermally active wall varies linearly in the x-direction, the heat transfer rate can be expressed in terms of Nusselt number at the x = 0 section, defined as

N u = 1 Δ T . (58)

where the temperature difference across the section is given by Δ T = T p ( 0, η ) T f ( 0, η ¯ ) the Nusselt number, Equation (57), after normalized with ε T / ( ε T η ¯ + η ) (such that Nu = 1 in pure conduction), give us

N u = ε T η ¯ + η Δ T ε T . (59)

Thus, upon substituting Equations (33) and (36) into Equation (58) it is found that

N u = ( ε T η ¯ + η + ε T R a C 2 A 6 ε T η ¯ + η ) 1 (60)

where

A 6 = γ ε T ( η 3 6 B 1 η 2 2 + B 2 η ) ( η ¯ 5 120 + A 1 η ¯ 4 24 + A 2 η ¯ 3 6 + A 3 η ¯ 2 2 + A 4 η ¯ ) . (61)

The Nusselt number derived in the present study, Equation (59), can be simplified upon considering the following limiting cases.

5.4.1. Case of Pure Anisotropic Porous Layer: ( η = 1 )

In the limit of a pure porous layer the Nusselt number, as predicted by Equation (59), can be reduced to the form:

N u = ( 1 6 + 10 b D a R a ) 1 . (62)

5.4.2. Case of Pure Isotropic Porous Layer: ( η = 1 )

When the permeability is the same in all directions (i.e. for an isotropic porous layer), we have: K 1 = K 2 and b = K * cos 2 ( φ ) + sin 2 ( φ ) = 1 and taking, D a R a = R , we obtain

N u = ( 1 6 + 10 R ) 1 . (63)

Which is in agreement the with the result predicted in the past by Alloui et al. [14] and Vasseur et al. [22].

5.4.3. Case of Pure Fluid Layer ( η = 0 ):

The limit of a pure fluid layer is also predicted by Equation (59).

N u = ( 1 R a C 2 ( 100 + 675 D a * + 1800 D a ) 72000 ( 1 + 3 D a * + 3 D a ) ) 1 . (64)

6. Results and Discussion

The effect of varying of β 1 , the slip parameter, of K * , the permeaility ratio and of η , the dimensionless position of the interface on the critical Rayleigh number R a c is illustraded in Figure 2 for D a = 10 3 , ε T = 1 and φ = 0 . It observes that when the slip parameter increases, the critical Rayleigh number for the onset of convection increases too. For η = 0 (corresponding to a single layer of fluid), the curve evolutes smoothly between the limit R a c = 320 (corresponding to the critical Rayleigh number for a fluid layer with a solid horizontal lower boundary and a free upper surface), when β 1 is small towards the limit R a c = 720 (corresponding to the critical Rayleigh number for a fluid layer bounded by two solid horizontal walls), for large values of β 1 . These values are identical to the results reported in the past, by authors as Vasseur et al. [4], Degan et al. [14], Alloui et al. [27] The limit η = 1 (not presented in Figure 2), corresponds to the case of a horizontal porous layer for which the critical Rayleigh number for the apparition of the convection is given by R a c = 12 b / D a (Equation (51)). When the dimensionless position of the inteface, η , is increased from η = 0 to η = 1 , it is observed that the influence of β 1 becomes less and less important. In the limiting case η = 1 , corresponding to a pure porous layer of Darcy, the slip parameter β 1 doesn’t have any effect. Figure 2 also shows that the anisotropiy ratio K * has a great influence on the critical Rayleigh number. Indeed, the anisotropiy ratio K * increases with an increase of the position, η of the interface. Thus, upon increasing η , the effect of anisotropy becomes more and more important. Moreover, compared to isotroic situation ( η = 1 the critical Rayleigh number for the onset of convection is enhanced when K * > 1 , and decreased when K * < 1 . To this effect, we read in Figure 2 that when η = 0.8 and β 1 = 100 , the critical Rayleigh number is R a c = 3915 for K * = 0.8 , R a c = 4688 for K * = 1 and becomes R a c = 5453 for K * = 1.2 . It results then, when the principal axes of anisotropy are oriented

Figure 2. Effect of the slip parameter, β1, of the permeability ratio, K*, and of the dimensionless position of the interface, η, on the critical Rayleigh number Rac for Da = 10−3, εt = 1 and φ = 0˚.

in a direction coinciding with the coordinate axes (i.e. φ = 0 ), the critical Rayleigh number increases (or decreases) when the permeability in the vertical direction ( K 1 ) is higher (or smaller) than the permeability in the horizontal direction ( K 2 ).

The effect of the anisotropic angle φ and of the anisotropic permeability ratio K * = K 1 / K 2 on the onset of convection ( R a c ) is illustrated in Figure 3 when D a = 10 4 , β 1 = 1 , η = 0.3 and ε T = 0.01 . It is seen that for a given value of the anisotropy ratio K * , the critical Rayleigh number increases with an increase of the anisotropic angle φ when the permeability ratio is made smaller than 1 (i.e. K * < 1 ). It explains by the fact that, because of the dependence of the critical Rayleigh number with the anisotropic parameter b = K * cos 2 ( φ ) + sin 2 ( φ ) and of the fixation of the position of the interface between the two layers, for K * < 1 , b becomes K * cos 2 ( φ ) + sin 2 ( φ ) = 1 > 0 , independently of the anisotropic angle, φ . Consequently, when φ varies from φ = 0 to φ = 90 , the influence of the permeability ratio K * becomes less and less important, and thus the onset of convection is made progressively enhanced. For the limiting case corresponding to φ = 90 which means that the principal axis with permeability ( K 2 ) is oriented parrallel to the gravity, (i.e., K 2 > K 1 ), the permeability ratio K * doesn’t have any effect. This is explained by the fact that, for φ = 90 we have b = K * independently of φ and R a c . For φ = 90 , the critical Rayleigh number is maximal and minimal when φ = 0 . It is evident from the analysis of this result that the characteristic parameter of the onset of convection is maximal (or minimal) when the main axis having the most elevated permeability of the porous layer is parallel (or perpendicular) to the gravity.

The effect of the conductivity thermal ratio, ε T , and the permeability ratio, K * , on the critical Rayleigh number, R a c , for η = 0.5 , β 1 = 0.1 and φ = 10 and for different values of the Darcy number , Da, is illustrated in Figure 4. It is clear that, independently of the Darcy number, a growth of the conductivity thermal ratio ε T also corresponds to a linear growth of the critical Rayleigh number, R a c . It results by the fact that the growth of ε T implies a reduction of the thermal conductivity ratio in the fluid layer situated above the porous layer. This effect stabilizes the systme, and require a more elevated value of the critical Rayleigh number consequently for onset of movements convection in all the cavity. These results are in agreement with those gotten in the past by Alloui et al. [14]. By somewhere else, the gotten results (Figure 4) show a great influence of the permeability ratio, K * , on the onset of convection. Thus, for a value given of the thermal conductivity, ε T , and of the Darcy number Da, the critical Rayleigh number is maximal when the permeability ratio, K * , is superior to the unity (i.e. K * > 1 ) and is minimal when K * is lower to the unity (i.e. K * < 1 ). For D a = 10 9 , the effect of the anisotropy disappears. The Darcy number, being a parameter that measures the permeability of the porous layer, then it is normal that when the Darcy number, Da, stretches toward zero, the phenomenon of convection disappears and makes room to the phenomenon of conduction.

Figure 3. Effect of the anisotropic angle φ and of the anisotropic permeability ratio K*=K1/K2 on the criteria of onset of convection Rac when Da = 10−4, β1= 1, η = 0.3 and εT = 0.01.

Figure 4. Effect of the conductivity thermal ratio, εT, and the permeability ratio, K*, on the critical Rayleigh number, Rac , for η = 0.5, β1 = 0.1 , φ = 10˚ and for different values of the Darcy number , Da.

In this case (that means D a 0 ) we cannot speak anymore of anisotropy.

The effect of Rayleigh number, Ra, on the Nusselt number, Nu, is presented in Figure 5 for different values of the position of the interface, η , when D a = 10 6 , β 1 = 1 , ε T = 1 , K * = 0.01 and φ = 0 . The results show that for a given value of η , there is a critical Rayleigh number R a c , Equation (49), below which convection does not occur. When the Rayleigh number is above this critical value, there is onset of convection and the heat transfer increases consequently. So for η = 0 , the heat transfer begin only from R a = R a c = 320 which corresponds to the critical Rayleigh number for a fluid layer with a solid

Figure 5. Effect of Rayleigh number, Ra, on the Nusselt number, Nu, for different values of the position of the interface for Da = 10−6, β1= 1, εT = 1, K* = 0.01 and φ = 0˚.

horizontal lower boundary and a free upper surface. these results are in agreement with those gotten in the past by Alloui et al. [14]. However, Figure 5 clearly illustrates the fact that Nu does not increase monotonically, as expected in general, but tends asymptotically toward a constant value that depends upon η , when Da, β 1 , ε T , K * and φ are fixed. Thus, from Equation (58), the limit N u = 3.9 is gotten for η = 0 , N u = 2.3 for η = 0.2 , N u = 1.6 for η = 0.3 and N u = 6 for η = 1 .

The effect of the permeability ratio, K * , and of Rayleigh number Ra on the Nusselt number is presented in Figure 6 for D a = 10 3 , β 1 = 1 , ε T = 1 , and φ = 40 . These results show that in pure porous layer (i.e. η = 1 ), when the Rayleigh number is lower the critical Rayleigh number given by R a c = 12 b / D a , there is no heat transfer. On the other hand, for a value superior to R a c there is heat transfer and the Nusselt number increases when Ra increases from R a c and tends toward the limit where N u = 6 . Results show that the infuence of K * is very significant because of the fact that, for a value given of Ra the rate of heat transfer becomes more and more important when K * becomes more and more small that the unity. Besides, the rate of heat becomes more and more small when K * becomes more and more raised in relation to the unity.

The evolution of the Nusselt number according to the position of the interface between the fluid layer and the porous layer, for D a = 10 2 , R a = 10 6 , K * = 0.1 , β 1 = 1 , ε T = 1 , and φ = 25 , is illustrated in Figure 7. The variation of the position of the interface, η , do 0 to 1 entail a variation of the Nusselt number respectively of N u = 5 (for a pure fluid layer corresponding to η = 0 ) to N u = 6 (for a porous layer corresponding to η = 1 ). However, the curve indicates that the transition between these two limits is not uniform. Thus, for η = 0 , point (1), the intensity of the speed in the fluid layer varies from the value

Figure 6. Effect of the permeability ratio, K*, and of Rayleigh number Ra on the Nusselt number for Da = 10−3; β1= 1, εT = 1, and φ = 40˚.

Figure 7. Evolution of the Nussselt number according to the position of the interface between the fluid layer and the porous layer when Da = 10−2; Ra = 106; β1= 1, εT = 1, and φ = 25˚.

u = 33.72 to the value u = 25.77 . This speed of slip is owed to the condition of Joseph-Beavers applied to the interface. When the position of the interface, η , increases, the Nusselt number decreases first and reaches the minimal value in conduction, for which N u = 1.00 at η = 0.1 ; it is materialized by the point (2), for which the profile of speed indicates a weak discontinuity to the interface. Above this value of η , the Nusselt number increases to reach the value N u = 1.4 at the position η = 0.3 ; it is materialized by the point (3), for which the profile of speed indicate a weak discontinuity also to the interface. Then, the heat transfer decreases constantly while passing by the point (4) to the minimal value N u = 1 characteristic of the thermal conduction for η = 0.8 ; it is represented by the point (5), for which the profile of speed indicates a relatively important discontinuity to the interface. Finally for η = 1 , the point (7) indicates that the linear speed profile is typical and is in agreement with the results gotten in the past by Alloui et al. [14] for the case of a porous layer. For this last situation the heat transfer reaches his maximum. This behavior of the Nusselt number and the profile of speed is like the results gotten in the past by Alloui et al. [14]. The difference between results gotten by these last and the present results are that the pofile of speed, whatever the position of the interface, start with a nill value ( u = 0 ) at the impermeable upper surface for their study while for the present results the prfile of speed takes his maximal value to the free surface. This behavior of the speed profile to the free surface is explain by the fact that the free surface positioned to y = η ¯ of the origine, O, of axis of coordinates, is permeable and that tension owed to the shearing is null, i.e., τ = μ ( d u f / d y ) = 0 ; the speed is then maximal. It results some therefore that to the upper free surface, the speed is constant and normal to this surface.

The evolution of the Nussselt number accordng to the position of the interface between the fluid layer and the porous layer and for different values of K * when D a = 10 2 , R a = 10 6 , β 1 = 1 , ε T = 1 , and φ = 25 , is illustrated in Figure 8. The results show that for η = 0 (which corresponds to a pure fluid layer where the anisotropy is out subject), the curves leave from only one point N u = 5 and when η begins by increasing, the effect of anisotropy becomes more pronounced until every curve reaches for η = 1 a value that depends on K * . Variations of the Nusselt number are identical to tose described previously for a value given of K * ( K * = 0.1 ). For a small value of the anisotropic permeability ratio K * , in relation to the unity, results show that curves tends towar the limit where N u = 6 . It explains by the fact that for a raised values of the Rayleigh number ( R a = 10 6 ) and for η = 1 , the Nusslt number tends toward the value N u = 6 .

Figure 9(a), Figure 9(b) and Figure 10(a), Figure 10(b) represent, respectively, the horizontal speed and the temperature distribution on all the width of the conduct, when β 1 = 10 , ε T = 1 , η = 0.5 , D a = 25 × 10 5 , φ = 0 and different values of K * and of Da. Speeds in the fluid layer and in porous layer are normalized like follows: u f * = ( u f / ψ 0 ) × 10 3 and u p * = ( u p / ψ 0 ) × 10 3 . Alike, temperature distributions in fluid layer and in porous layer are ormaliwed like follows: θ f * = ( ( θ f + y ) / ( C ψ 0 ) ) × 10 3 and θ p * = ( ( ε T θ p + y ) / ( C ψ 0 ) ) × 10 3 . For every value of K * , Figure 9(a), Figure 9(b) indicates that the speed to the upper free surface is constant and normal to this surface, and that the constraint of shearing to this surface is equal to zero. We note that the slip parameter, β 1 , and the Darcy number measuring the permeability of the porous layer, have a strong influence on fields of speed and temperature. Thus, for D a = 2.5 × 10 5 (correspondent to a porous layer very little permeable), Figure 9(a) indicates that the intensity of flow inside the system is relatively weak. In fact, the major part of the flow is confined in the fluid layer ( 0 y 0.5 ) where we note a light effect of the anisotropy, while in the porous layer the flui is nearly to rest. The transfer of the heat in the porous layer is mainly a conduction. Besides, we

Figure 8. Evolution of the Nussselt number according to the position of the interface between the fluid layer and the porous layer and for different values of K* when Da = 10−2; Ra = 106; β1= 1, εT = 1, and φ = 25˚.

(a) (b)

Figure 9. (a) Horizontal speed distribution on all the width of the condut, when β1= 10, εT = 1, η = 0.5, Da = 2.5 × 10−5, φ = 0˚ and for different values of K*; (b) Horizontal speed distribution on all the width of the condut, when β1= 10, εT = 1, η = 0.5, Da = 10−4, φ = 0˚ and for different values of K*.

(a) (b)

Figure 10. (a) Horizontal temperature distribution on all the width of the condut, when β1= 10, εT = 1, η = 0.5, Da = 2.5× 10−5, φ = 0˚ and for different values of K*; (b) Horizontal temperature distribution on all the width of the condut, when β1= 10, εT = 1, η = 0.5, Da = 10−34, φ = 0˚ and for different values of K*.

observe that, when Da is very small, the effect of the anisotropic permeaility ratio K * is nearly negligible as indicate in Figure 9(a) and Figure 10(a). When the Darcy number increases ( D a = 10 4 ), the flow of the fluid intensifies slightly in the two layers (fluid and porous) and depends considerably on the permeability ratio, K * . Thus, to the free surface of the fluid, according to the curve Figure 9(a): | u f * | = 2.82 for K * = 0.5 , | u f * | = 2.72 for K * = 1.0 and | u f * | = 2.63 for K * = 10.0 . For the curve Figure 10(b), we have to the free surface of the fluid: | θ f * | = 0.35 for K * = 0.5 , | θ f * | = 0.28 for K * = 1.0 and | θ f * | = 0.22 for K * = 5.0 . These results show that the effect of anisotropy becomes important when the Darcy number measuring the permeability of the porous layer becomes more and more elevated. Otherwise, the profile of speed in the fluid layer (Figure 9(a)), for every value of K * , reaches his maximum in a point situated on a parallel line out of place downwards in relation to the central axis of the fluid layer ( y = 0.1234 ) and decreases drastic way toward the interface letting appear a dynamic limits layer. This shift is due to the presence of the nominal surface where the speed should change profile, since changing of middle. These results are compliant those gotten by Alloui et al. [14] who studied the effect of Darcy and the slip parameter to the interface.

In Figure 11(a), Figure 11(e) and Figure 12(a), Figure 12(e) are drawn profiles of temperature and speed so much in fluid layer that in porous layer according to the coordinate y and for different values of the permeability ratio K * and of the position of the interface, η , when D a = 10 1 , β 1 = 1 , ε T = 1 , and φ = 0 . When η = 0 , Figure 11(a) and Figure 12(a) show only one curve for profiles of speed and temperature. It justifies by the fact that in pure fluid layer the anisotropy is out subject. When η becomes more and more raised, the corresponding curves (Figure 11(b), Figure 11(e) and Figure 12(b), Figure 12(e)) show that the effect of anisotropy on the flow becomes accentuated more and more. The same interpretations are made on curves representing fonctions of current (Figure 13(a), Figure 13(e)) when D a = 10 3 , K * = 0.5 , β 1 = 10 and ε T = 1 .

7. Conclusions

In this investigation, an analytical study of heat transfer is conducted in a system consisting of a horizontal fluid layer over a saturated porous bed. Our research concerns the influence of hydrodynamic anisotropy on stability geothermal streams. The fluid layer superposed on the porous layer is heated from below. Using the Navier-Stokes model for the fluid layer and the Darcy model for the porous layer, an exact solution is found for a fully developed system of forced convective flow through the superposed layers. The Beavers-Joseph condition is applied at the interface between the two layers. The main conclusions of the present study are:

It appears that when the principal axes of anisotropy are oriented parallel to the coordinate axis ( φ = 0 ), the characteristic parameter of the criterion of onset of convection, Rac, increases (or decreases) when the permeability in the vertical direction ( K 1 ) is greater (or smaller) than the permeability in the horizontal direction ( K 2 ).

Figure 11. Profiles of temperature so much in fluid layer that in porous layer according to the coordinate y and for different values of the permeability ratio K* and of the position of the interface, η, when Da = 0.1, φ = 0˚ β1= 1, and εT = 1.

Figure 12. Profiles of speed so much in fluid layer that in porous layer according to the coordinate y and for different values of the permeability ratio K* and of the position of the interface, η, when Da = 0.1, φ = 0˚, β1= 1, and εT = 1.

When the permeability K * increases, the critical Rayleigh number increases.

The characteristic parameter of the criterion of occurrence of convection Rac is maximum (or minimum) when the orientation of the main axis of the high permeability of the anisotropic porous layer is parallel (or perpendicular) to the gravity.

Figure 13. Profiles of fonctions of current so much in fluid layer that in porous layer according to the coordinate y and for different values of φ and of the position of the interface, η, when Da = 0.001, β1 = 10, K* = 0.5 and εT = 1.

For a given value of Ra, the rate of heat transfer becomes increasingly important when K * is more than unity. In addition, the heat transferred becomes smaller when K * is increasingly large compared to unity.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Degan, G., Yovogan, J., Fagbémi, L. and Allou, Z. (2019) Stability of Geothermal Convection in Anisotropic River Beds. Engineering, 11, 343-365. https://doi.org/10.4236/eng.2019.117026

References

  1. 1. Nield, D.A. and Bejan, A. (2006) Convection in Porous Media. 3rd Edition, Springer, New York.

  2. 2. Somerton, C.W. and Catton, I. (1982) On the Thermal Instability of Superposed Porous and Fluid Layers. Journal of Heat Transfer, 104, 160-165. https://doi.org/10.1115/1.3245044

  3. 3. Nishimura, T., Takumi, T., Shiraishi, M., Kawamura, Y. and Ozoe, H. (1986) Numerical Analysis of Natural Convection in a Rectangular Horizontally Cavity Divided into Fluid and Porous Region. International Journal of Heat and Mass Transfer, 26, 889-898. https://doi.org/10.1016/0017-9310(86)90184-5

  4. 4. Vasseur, P., Wang, C.H. and Sen, M. (1989) Thermal Instability and Natural Convection in a Fluid Layer over a Porous Substrate. Waerme-Stoffuebertragung, 24, 337-347. https://doi.org/10.1007/BF01592780

  5. 5. Chen, F., Chen, C.F. and Pearlsein, A.J. (1991) Convective Instability in Superposed Fluid and Anisotropic Porous Layers. Physics of Fluids, 4, 556-565. https://doi.org/10.1063/1.858117

  6. 6. Shivakumara, I.S., Suma, S.P. and Chavaraddi, K.B. (2006) Onset of Surface-Ten-sion-Driven Convection in Superposed Layers of Fluid and Saturated Porous Medium. Archives of Mechanics, 58, 71-92.

  7. 7. Jones, I.P. (1973) Low Reynolds Number Flow Past a Porous Spherical Shell. Proceedings of the Cambridge Philosophical Society, 73, 231-238. https://doi.org/10.1017/S0305004100047642

  8. 8. Valenca-Lopez, J.J. and Ochoa-Tapia, J.A. (2001) A Study of Buoyancy-Driven Flow in a Confined Fluid Overlying a Porous Layer. International Journal of Heat and Mass Transfer, 44, 4725-4736. https://doi.org/10.1016/S0017-9310(01)00105-3

  9. 9. Alazmi, B. and Vafai, K. (2001) Analysis of Fluid Flow and Heat Transfer Interfacial Conditions between a Porous Medium and a Fluid Layer. International Journal of Heat and Mass Transfer, 44, 1735-1749. https://doi.org/10.1016/S0017-9310(00)00217-9

  10. 10. Hirata, S.C., Goyeau, B. and Gobin, D. (2007) Stability of Natural Convection in Superposed Fluid and Porous Layers: Influence of the Interfacial Jump Boundary Condition. Physics of Fluids, 19, Article ID: 058102. https://doi.org/10.1063/1.2730877

  11. 11. Hirata, S.C., Goyeau, B., Gobin, D., Carr, M. and Cotta, R.M. (2007) Linear Stability of Natural Convection in Superposed Fluid and Porous Layers: Influence of the Interfacial Modeling. International Journal of Heat and Mass Transfer, 50, 1356-1367. https://doi.org/10.1016/j.ijheatmasstransfer.2006.09.038

  12. 12. Hirata, S.C., Goyeau, B. and Gobin, D. (2009) Stability of Thermosolutal Natural Convection in Superposed Fluid and Porous Layers. Transport in Porous Media, 78, 525-536. https://doi.org/10.1007/s11242-008-9322-9

  13. 13. Hirata, S.C., Goyeau, B., Gobin, D., Chandesris, M. and Jamet, D. (2009) Stability of Natural Convection in Superposed Fluid and Porous Layers: Equivalence of the One- and Two-Domain Approaches. International Journal of Heat and Mass Transfer, 52, 533-536. https://doi.org/10.1016/j.ijheatmasstransfer.2008.07.045

  14. 14. Alloui, Z. and Vasseur, P. (2010) Convection in Superposed Fluid and Porous Layers. Acta Mechanica, 214, 245-260. https://doi.org/10.1007/s00707-010-0284-y

  15. 15. Yovogan, J. and Degan, G. (2013) Effect of Anisotropic Permeability on Convective Heat Transfer through a Porous River Bed Underlying a Fluid Layer. Journal of Engineering Mathematics, 81, 127-140. https://doi.org/10.1007/s10665-012-9605-6

  16. 16. Yovogan, J., Degan, G. and Fagbemi, L. (2018) Effect of Constant Magnetic Field on Convective Heat Transfer through Anisotropic River Beds. Journal of Crystallization Process and Technology, 8, 57-71. https://doi.org/10.4236/jcpt.2018.82004

  17. 17. Beavers, G.S. and Joseph, D.D. (1967) Boundary Conditions at a Naturally Permeable Wall. Journal of Fluid Mechanics, 30, 197-207. https://doi.org/10.1017/S0022112067001375

  18. 18. Rudraiah, N. and Veerabhadraiah, R. (1977) Temperature Distribution in Couette Flow past a Permeable Bed. Proc. Indian Academy of Sciences, 86, 537. https://doi.org/10.1007/BF03046910

  19. 19. Rudraiah, N. and Veerabhadraiah, R. (1978) Effect of Buoyancy on the Free Surface Flow past a Permeable Bed. Waerme-Stoffuebertragung, 11, 265. https://doi.org/10.1007/BF02587790

  20. 20. Cormack, D.E., Leal, L.G. and Inberger, J. (1974) Natural Convection in a Shallow Cavity with Differentially Heated End Walls. Part 1. Asymptotic Theory. Journal of Fluid Mechanics, 65, 209-230. https://doi.org/10.1017/S0022112074001352

  21. 21. Sparrow, F.M., Goldstein, R.J. and Jonsson, V.K. (1964) Thermal Instability in a Horizontal Fluid Layer: Effect of Boundary Conditions and Nonlinear Temperature Profile. Journal of Fluid Mechanics, 18, 513-528. https://doi.org/10.1017/S0022112064000386

  22. 22. Vasseur, P., Wang, C.H. and Sen, M. (1989) The Brinkman Model for Natural Convection in Shallow Porous Cavity with Uniform Heat Flux. Numerical Heat Transfer, 15, 221-242. https://doi.org/10.1080/10407788908944686

  23. 23. Sen, M., Vasseur, P. and Robillard, L. (1987) Multiple Steady States for Unicellular Natural Convection in an Inclined Porous Layer. International Journal of Heat and Mass Transfer, 30, 2097-2113. https://doi.org/10.1016/0017-9310(87)90089-5

  24. 24. Sen, M., Vasseur, P. and Robillard, L. (1988) Parallel Flow Convection in a Tilted Two-Dimensional Porous Layer Heated from Below. Physics of Fluids, 31, 3480-3487. https://doi.org/10.1063/1.866916

  25. 25. Bahloul, A., Vasseur, P. and Robillard, L. (2007) Convection of a Binary Fluid Saturating a Shallow Porous Cavity Subjected to Cross Heat Fluxes. Journal of Fluid Mechanics, 574, 317-342. https://doi.org/10.1017/S0022112006004113

  26. 26. Nield, D.A. (1968) Onset of Thermohaline Convection in a Porous Medium. Water Resources Research, 4, 553-560. https://doi.org/10.1029/WR004i003p00553

  27. 27. Degan, G. and Vasseur, P. (2003) Influence of Anisotropy on Convection in Porous Media with Nonuniform Thermal Gradient. International Journal of Heat and Mass Transfer, 46, 781-789. https://doi.org/10.1016/S0017-9310(02)00352-6