^{1}

^{*}

^{2}

^{*}

^{1}

^{*}

^{3}

^{*}

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 .

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 [

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 [

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.

The physical model illustrating the problem under different considerations is shown in

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 .

The equations governing the conservation of mass, momentum, energy and the Boussinesq approximation (see, Alloui et al. [

• 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 [

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 [

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)

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. [

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)

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 ) = C x + 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)

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)

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.

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 [

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. [

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. [

R a c = 320. (55)

which is also a value reported in the past by Alloui et al. [

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. [

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.

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.

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)

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. [

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)

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

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

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

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

horizontal lower boundary and a free upper surface. these results are in agreement with those gotten in the past by Alloui et al. [

The effect of the permeability ratio, K * , and of Rayleigh number Ra on the Nusselt number is presented in

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

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. [

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

observe that, when Da is very small, the effect of the anisotropic permeaility ratio K * is nearly negligible as indicate in

In

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, Ra_{c}, 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 ).

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

The characteristic parameter of the criterion of occurrence of convection Ra_{c} 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.

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.

The authors declare no conflicts of interest regarding the publication of 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