On double diffusion in a Brinkman heat generating porous layer

INT. COMM. HEAT M A S S T R A N S F E R Vol.

12, pp. 149-158,

1985

~P~rganDnPress

0735-1933/85 $3.00 + .00 Ltd. Printed in the United States

ON DOUBLE DIFFUSION IN A BRINKMAN HEAT GENERATING POROUS LAYER

Bill Selimos and Dimos Poulikakos University of Illinois at Chicago Department of Mechanical Engineering P.O. Box 4348 Chicago, Illinois 60680

( C ~ i ~ n i c a t e d by J.P. Hartnett and W.J. Minkowycz)

ABSTRACT This paper relies on linear M%~bility analysis to define the parametric domain in which c0nve~ti~e motions exist in a horizontal unstably stratified p o r o ~ layer. The paper focuses on the combined effect of three distinct factors on the onset of convection: the presence of a second diffusing component (in addition to heat), the presence of internal heat generation in the porous medium and the departure from the Darcy model of flow in porous medittm by using Brinkman's modification of Darcy's law. Introduction The classical "Benard" problem in porous media is the determination of the onset of convective instabilities

in a fluid saturated porous layer

heated from below and cooled from above so that a linear temperature profile is established in the pure conduction state. stability analysis)

The first theoretical

(linear

solutions of the above problem were obtained by Horton

and Rogers EI~ and were motivated by geophysical

considerations.

During

recent years, a number of studies aimed to extent the results in E1,23 and to establish the effect of a number of phenomena on the critical Rayleigh number describing the onset of convection.

The majority of these investi-

gations are concerned with the following three phenomena as they relate to the onset of convection in a porous layer:

the presence of internal heat

generation in the porous matrix

[3-53, the presence of a second diffusing

component

[6-83, and the use of flow models for the

(in addition to heat)

porous layer valid for a wide range of pore diameters,

to improve the pre-

dictions of the commonly used Darcy model, which is valid as long as the size of the pore is relatively small compared to the dimensions of the layer [ 9 - I ~ 149

150

B. Selimos and D. Poulikakos

Vol.

12, No. 2

The main objective of the present study is to address the topic of the combined effect of the above three distinct factors on the onset of convection.

Eventhough the effect of internal heat generation or the effect of

double diffusion or the effect of Brinkman's modification

alone on the onset

of convection has been investigated [3-113, to the authors'

knowledge the

combined effect of all three phenomena has not been determined yet.

The

present paper relies on the formalism of linear stability analysis to study the impact of the above three mechanisms, marginal stability boundaries overstability

acting simultaneously,

for stationary convection.

on the

The phenomenon of

(which is a possibility due to the presence of a second dif-

fusing component)

has not been investigated Mathematical

The configuration

in the course of this study.

Formulation

of interest consists of a heat generating horizontal

porous layer bounded from above and below by two impermeable boundaries. These boundaries

are kept at constant temperature and species concentration.

According to Brinkman's modification

of Darcy's flow model El0,113,

equations describing the conservation of mass, momentum,

the

energy and species

at each point in the porous layer are: V • u = 0 VP =

(i) -

(2)

~T ~ ~ O ~ + u " VT = d V2T + ( O C p ) f l u i d

(3)

~~S

[

+

~

•

(T-To) + 7

?S

=

V2s

D

(4)

s

The various symbols in the above equations are defined in the nomenclature. The effective viscosity in the porous bed, ~, may not, in general, be e q u a l to the fluid viscosity.

Various models EI2~ have been postulated

to determine the effective viscosity.

in order '

Most of these models are elaborate

and it has been found that simply letting ~ = ~ reproduces well experimental data E133.

It is worth notinglthat

side of equation

the last two terms in the right hand

(2) imply that when the permeability K is large, the re-

sistance to flow becomes effectively controlled by the ordinary viscou~ resistance

(~V2~.

in a fluid layer.

In this case the convection phenomenon is similar to that According to the Boussinesq approximation

the density was

assumed constant everywhere except in the body force term in equation where it was approximated by 0 = Qo [I - ~
(2)

VOI.

12, NO. 2

BRINKMANH~tTG~%~TINGPOROUSIAYER

151

Since the scope of the p r e s e n t study is limited to define the m a r g i n a l stability b o u n d a r i e s

for stationary convection,

the left hand sides of eqs. stability b o u n d a r i e s state equations.

(3) and

the time d e p e n d e d terms in

(4) are omitted and the m a r g i n a l

are o b t a i n e d via linear stability analysis of the steady

The t e m p e r a t u r e

wall are assumed constant.

and the species c o n c e n t r a t i o n

The notation adopted in the linear stability

analysis of the next section assumes that if dT = T L- T U > the respective

on either

temperature or c o n c e n t r a t i o n

0or dS=

S L- S U < 0

gradient is destabilizing.

Linear Stability A n a l y s i s A c c o r d i n g to standard linear stability p r o c e d u r e s we assume that all the depended variables of the p r o b l e m can be d e c o m p o s e d a perturbation

component

(T = T + @, S = S + c).

p r o b l e m are easy to o b t a i n E14 3 . The p e r t u r b a t i o n

into a mean

(static) and

The solutions to the static

Hence, they are emitted here for brevity.

q u a n t i t i e s can be w r i t t e n as

,c. ,,,)

(w.,o.c.). (w. ,z.,.

[,x.+

The resulting g o v e r n i n g equations read [14 3 -Rn 2 Q, + ~ e ~2RsC. = ,

-

_

-

[~(~-z.) - 1] W,-(D2- ~ .

(7,

Note that the above equations have been m a d e d i m e n s i o n l e s s

according to the

following definitions: x.,y.,z. =

@ @. = ~ ,

(x,y,z)/H,

The n o n d i m e n s i o n a l i z a t i o n

c w C, = ~-, w. = ~ H

(9)

of the p r o b l e m y i e l d e d five d i m e n s i o n l e s s para-

meters R, Rs, B, Le and Da.

All these p a r a m e t e r s

are d e f i n e d in the no-

menclature.

In addition, n 2 2 + n2 is the h o r i z o n t a l wave number and d = ~x y operator D = d~z, is shorthand notation for the d e r i v a t i v e w i t h respect to z,. In the p r e s e n t study we investigate b o t h the case where the u p p e r surface is free and the case where the upper surface is rigid.

Mathematically,

the b o u n d a r y conditions are expressed by the following equations: z, = 0

:

W. = DW, = Q, = C, = 0

rigid lower surface

z, = 1

:

G, = C, = 0 and W, = D2W. = 0 for free upper b o u n d a r y or W. = DW, = 0 for rigid upper b o u n d a r y

(I0)

152

B. Selimos and D. Poulikakos

Vol. 12, No. 2

To obtain the critical Rayleigh number, at which convection first sets in the porous layer, we employ the Galerkin steps are reported in E I4~.

method E 15-173.

The method yields a system of n algebraic

equations for the unknown constants A I F m I (n ,m) BR] | o ~A m ~2 [ R - R s ---J2 L- (m~)2 + ~2 Da------[ + (m~)2 + ~ 2 + CmISn)+

D I(n)+ m 3

+ BmI(n)6 + C I ( n 7)+m

The successive

E I(n~ m# j

(n) + Bm I 1

P - (n ,mY + BRN 2 l±s ~ -2 ] [ (m~)2+ D2+Da[(m~) 2 + )~ 1

DmI: n) EmI:n)] + [ (rag)2 + D2] Io(n'm) I

It is worth noting that by applying the boundary conditions

(ii)

0

(I0) for each of

the two cases of interest, the constants Bm, C m, D m, Em are evaluated and are reported in the Appendix.

The Appendix also contains the definitions

of the integrals, denoted by I.± in (ii). The system of n equations represented by expression (ii) for n = i, 2, 3, ..., yields an nxm matrix of coefficients for the unknown constants A . In the present study the approximation n = m = 6 was found m to be sufficient to provide accurate results. In order that A be nonm trivial it is required that the determinant of the matrix of coefficients for Am vanishes.

This condition yields Rcr and Ncr for prescribed values of

R , B and Da E14~s Results and Discussion The main results regarding marginal stability for rigid-free boundaries are reported in Fig. i.

Figure la shows the impact of the heat generation

B on the critical Rayleigh number, R

. The value of the Darcy number, cr Da = 0.001 is representative of the Darcian porous medium limit. The curves labeled R

= 0 were found to be in excellent agreement with the results s reported in ref. E 4 3 which dealt with the problem of the onset of convection in a heat generating porous layer in the absence of concentration gradients (the disagreement was of the order of a fraction of one percent and it was due to the fact that the curves in ref. 4 were for smaller Darcy number, Da = 0).

Two families of curves are reported:

one for destabilizing AT.

one for stabilizing AT and

According to our notation, when AT is stabilizing,

the sign of the values of R negative.

read off the graph should be changed to cr With reference to the curves corresponding to destabilizing

values of AT, increasing the internal heat generation has the effect of de-

Vol. 12, NO. 2

~4

. . . . . . . .

[

. . . . . . . .

BRINKMANHEATG~qERATINGPOROUSIAYER

i

. . . . . . . .

r

. . . . . . . .

i

153

i

. . . . . .

6T~O

!'*i

AT~,O

.-,

.......

;e

...................

;

......

10 3 R©r

n~

(a)

(b)

~

~

,

I

, ' ' <.

I

,

........

I

........

l

........

I

........

"\

)1 t D ~ . O( I1 ~ s G :

•

li (d,,t,~zJ,~

~4

........

I

10

........

'iOs

~

,

, , .....

~

. . . . . .

|(c

10~

IOs

......

~o~

10-4

~

10-1

i0-I

I

I0

Oo

(C)

(d) FIG.

a)

........

1

C r i t i c a l Rayleigh number vs. heat g e n e r a t i o n function for rigid-free boundaries, Da = 0.001 b) Da = 0.i c) Da = i00 d) The t r a n s i t i o n from D a r c y p o r o u s m e d i a to classical fluids for rigid-free b o u n d a r i e s I

........

I

........

I

.......

i

......

'

~TCO

........

r

o'1

0.~o(I).

i

.......

1~ s

[ 1~ I

........

i

........

1

L

~o

j_,,

...... lo

m

~od ~,,j

~0s

lO'1

,

I

, ......

~0- I

~o-1

I~r

Ca)

Cb) FIG.

a)

1 D*

2

C r i t i c a l R a y l e i g h number vs. heat g e n e r a t i o n function for rigid-rigid boundaries, Da = 0.001 b) The t r a n s i t i o n from Darcy p o r o u s m e d i a to c l a s s i c a l fluids for r i g i d - r i g i d boundaries

10

RI. ~

154

B. S e l i m o s a n d D. P o u l i k a k o s

creasing

the critical

of the curves

Rayleigh

indicates

value of R after w h i c h

number.

Vol.

The sharp drop in the right

that for small values instabilities

appear

of B there

is to either expand (R = -20) or reduce s s regime. When the effect of AT is stabilizing, is n e e d e d

concentration heat

source

gradients

necessary

the i m p o s e d B where tion.

to d e s t a b i l i z e increase

for thermal

AT becomes

the internal

heat g e n e r a t i o n

This fact is shown

and negative generation

to reduce

Figures stability

the impact

Da = 0.i

(Fig.

and, next to Da = i00

matrix

has a negligible

stronger

la to Fig. internal

stabilities

that the Darcy

is satisfied

in d e s c r i b i n g

in Figs.

la-lc,

the right between

that,

heat g e n e r a t i o n

[113.

on the marginal

fluids

as the p o r o u s m a t r i x becomes to bring

This p h e n o m e n o n

about makes

does not satisfy

to flow.

Another

Da shifts

porous

limit).

less dense, in-

sense p h y s i c a l l y

the no-slip

con-

fluids

by ~ in eq.

interesting

the stability

the

Passing

convective

In the limit of classical

is that increasing

to

of the case where

since the term m u l t i p l i e d

the resistance

Increas-

(highly permeable)

(classical

is required

media

however,

(2) takes

observation

boundaries

to

of increasing R ). Figure id bridges the gap cr flows and flows of c l a s s i c a l fluids for B = 0 and i00.

It is safe to say that for Da = 0(i) on the onset criterion

of c o n v e c t i o n

or greater

is negligible.

the p o r o u s

Increasing

matrix

B renders

effect

the above

more conservative.

The last two Figures

2a,b illustrate

to the case of rigid-rigid

discussion well.

of positive

heat

(the d i r e c t i o n

porous

relevant

cr of destabiliza-

is increased

representative

flow model

on the lower b o u n d a r y

this c o n d i t i o n over

ic)

for the same Rcr and R s.

if we realize dition

of the Darcy number

of

and large

= R + 4 ~ 2 was verified. cr s on R s cr"

on the flow

ic reveals

R

1 (the curves

of a n o n - D a r c i a n

(Fig.

effect

The nature

of internal

of the Darcy number

ib) r e p r e s e n t a t i v e

matrix

from Fig.

the p r o c e s s

of the

[63 R of R

The value

strength

respectively.

dominates

ib and ic show the effect

boundaries.

or d e s t a b i l i z i n g

in the limit of small

In the absence

(B = 0), the c r i t e r i o n

ing B tends

Stabilizing

in the left side of Fig.

AT's come together).

The effect

= 20) the convection s stronger internal heat

the needed

convection,

less important

side

(R

the layer.

or decrease

2

is a distinct

in the system.

of R

generation

12, No.

for the rigid-free

Hence,

investigation.

for brevity, Detailed

boundaries.

boundary

of the investigation Most of the q u a l i t a t i v e

conditions

we will p r e s e n t

results

results

holds

in this case as

here only the highlights

are reported

in ref.

E143.

of the

It was observed

Vol. 12, NO. 2

BRINEM%NH~iT~POBOt~

155

that in the small Da limit the results for the two types of boundary conditions converge to the same values. becomes felt

A s the presence of the solid walls

(increasing Da) the results for the two sets of boundary

conditions deviate

(Figures id and 2b).

Of interest is the fact that the

curves for stabilizing and destabilizing AT "come together" at considerably higher R

values in the case of rigid,rigid boundaries (compare for example, cr Fig. la to Fig. 2a). The results of the present study for rigid-rigid boundaries showed excellent agreement with the findings of ref. [18~ in the classical fluids limit with heat generation and ref. R

= 0, B = 0.

[193

in the limit

Finally, the presence of a second solid boundary

(top wall)

s

yielded an "all around" increase in R

cr

.

Nomenclature A

constants in eq.

m

(ii)

dimensionless volumetric heat generation

B

constants in equation

B m

(ii)

c

perturbation concentration

c(z)

function in eq.

C

constants in equation

m

c

(5) (ii)

fluid specific heat at constant pressure P

D

mass diffusivity of the fluid saturated porous medium

s Da D E

Darcy number, K/H 2

m m

g

constants in equation

(11)

constants in equation

(11)

acceleration of gravity

H

layer depth

I. 1 K

(n,m)

ith intergral in equation

(11)

permeability of the porous matrix

Le

Lewis number, ~/'D s

P

pressure

Q

volumetric heat generation, QH2/(AAT)

R

thermal Rayleigh number, K~8~TH

S

concentration Rayleigh number, K ~ S H 9D species concentration s

AS

species concentration difference between upper and lower surface

R s

temperature AT

temperature time

difference between upper and lower surface

156

B. Selimos and D. Poulikakos

Vol. 12, No.

velocity vector, ~(u,v,w) U

x- component of velocity vector

V

y- component of velocity vector

W

z- component of velocity vector

W(z)

function in eq.

(5)

X

horizontal cartesian coordinate

Y

horizontal cartesian coordinate

Z

vertical cartesian coordinate

Greek symbols thermal diffusivity of the fluid saturated porous medium

B

thermal coefficient of volume expansion

Y

concentration coefficient of volume expansion

n

horizontal wave number

8

perturbation temperature

@(z)

function in eq. (5)

l

thermal conductivity of porous/fluid matrix fluid viscosity effective viscosity of porous bed fluid density

P

heat capacity ratio

[~ (0Cp)f + ( i - ~ ) ( 0 C p ) s ] ~ D C p ) f

porosity of porous matrix Subscripts L

lower surface

m

incremental integer of sine series solution

n

number of equations

u

upper surface

symbols dimensionless quantity static condition References i.

Horton, C.W. and Rogers, F.T., Jr., Journal of Applied Physics, 16, 367, (1945).

2.

Lapwood, E.R., Proc. Camb. Phil. Soc., 44, 508, (1948).

Vol. 12, No. 2

~ H E A T ~ P O R ~ J S I A Y E R

157

3.

Hwang, I., Ph.D Thesis, University of Minnesota, 1971.

4.

Gasser, R.D. and Kazimi, M.S., Journal of Heat Transfer, 98, 49, (1976).

5. Schulenberg, T. and Muller, U., Intl. Journal of Heat and Mass Transfer, 27, 677, (1984). 6. Nield, D.A., Water Resources Research, 4, 553, 7. Rubin, H., Water Resources Research, 2, 211,

(1968)

(1973)

8. Taunton, J.W. and Lightfoot, E.N., Physics of Fluids, 15, 748, (1972). 9. Katto, Y. and Masuoka, T., Intl. Journal of Heat and Mass Transfer, iO, 297, (1967). i0. Brinkman, H.C. Applied Science Research, ~, 27, (1947). ii. Cheng, P., Advances in Heat Transfer, 14, i, (1979). 12. Lundgren, T.S., Journal of Fluid Mechanics, 51, 273,

(1972).

13. G. Neale and W. Nader, Canadian Journal of Chemical Engineering, 52, 475, (1974). 14.

Selimos, B. On Double Diffusion in a Brlnkman Heat Generating Porous Layer, M.S. Thesis, University of Illinois at Chicago, Chicago, (1984).

15. Chandrasekhar, S., Hydrodynamic and Hydroma@netic Stability, Oxford University Press, (1961). 16. Kantorovich, L.V. and Krylov, V.I., Approximate Methods of Hi~her Analysis, John Wiley, N.Y., pp. 292-293, (1964). 17. Chung, T.J., Finite Element Analysis in Fluid Dynamics, McGraw-Hill, N.Y., (1978). 18. Sparrow, E., Goldstein, R. and Johnson, V., J. of Fluid Mechanics, 18, 513, (1964) 19. Walker, K. and Homsy, J.M., J. Heat Transfer, 99, 338,

(1977).

Appendix The coefficients B

- E , m = 1,2, ... in equation (ll) are as follws: m

m

a)

Bm = ~2~ Q 1

2

R i g i d - F r e e H o r i z o n t a l Boundaries

+ ~ ~ 2Q 1

2

b)

2

13Q ,

~3Q 1

-

2

2

1

1

Cm = Bm ~ ,

Dm = Bm ~ ,

1

2

Ri@id-Ri@id Horizontal Boundaries

Em

=

-D

m

158

B. Selimos and D. Poulikakos

Bm = Am [%IQ 3

Cm = A m

Vol. 12, No. 2

I 2Q % + ~ 2 + ~ 2 (-l)m ~cosh~ I - coshl2>]/Qs,

[X2Q3 _ % IQ ~ + % i _ % 1 (_[)m (cosh% 1 - c o s h % 2 ) ] / Q 5

Dm = A m

6

2

2

m

m

where A

=

m

-~ [[ - , ~2 [ ~ (m~) + + Da (m)2

= sinh%

Q3

1

sinh%

2

,

Q

~

+ ~

= cosh'

2] 2

1

QI

cosh%

2

= sinh%

,

Q

2

coshl

Q

i

2

= sin%

I

coshl

= -2% , + % % Q -<121 +%;)Q3 I 2 I z 4

5

Q6 = %2Q2 - % IQ I The integrals

in eq.

(ii) are defined as follows:

.i I(n)i =Io sin(n~z,) F i ( z ~ z ,

F

0

F F

= sin(m~z.), = cosh(%2z,),

7

= z,sinh(%2z,),

F

I

F

i = 0,i,

= sinh(% z,), i

F

.°.,9

2

= z, s i n ( m ~ z , ) ,

5 F

~

where

= sinh(% z.), 2 F

= z,cosh(% z,), i

6 F

9

F

3

= cosh(%iz,)

= z,sinh(llz,~ = z,cosh (~ z,) 2

Acknowledgments Financial

support provided by the University

for this work is greatly appreciated.

of Illinois Research Board

2