Effects of a random porosity model on double diffusive natural convection in a porous medium enclosure

(1)

Pergamon

PII S0735-1933(00)00090-7

EFFECTS OF A RANDOM POROSITY MODEL ON DOUBLE DIFFUSIVE NATURAL CONVECTION IN A POROUS MEDIUM ENCLOSURE

Wu-Shung Fu and Wen-Wang Ke Department of Mechanical Engineering

National Chiao Tung University Hsinchu, Taiwan 30010

R.O.C.

(Communicated by J.P. Hartnett and W.J. Minkowycz) ABSTRACT

A double diffusive natural convection in a rectangular enclosure filled with porous medium is investigated numerically. The distribution of porosity is based upon the random porosity model The Darcy-Brinkman-Forchheimer model is used and the factors of heat flux, mean porosity and standard deviation are taken into consideration. The SIMPLEC method with iterative processes is adopted to solve the governing equations.

The effects of the random porosity model on the distributions of local Nusselt number are remarkable and the variations of the local Nusselt number become disordered. The contribution of latent heat transfer to the total heat transfer of the high Rayleigh number is larger than that of the low Rayleigh number and the variations of the latent heat transfer are not in order. © x000 Elsevier Science Ltd

lalraflltC.ti~

Porous structures are very often applied to enhance heat transfer rate. For analyzing heat transfer in porous medium, a simplified model of constant porsity was first proposed in the past to analyze the phenomena o f flow and thermal fields in the porous. However, the porosity which varied significantly in the near wall region was observed by Roblee et al. [1] and

(2)

120 W.-S. Fu and W.-W. Ke Vol. 27, No. 1

Benenati and Brosilow [2] in their experimental results. Further, Cheng et al. [3] indicated that the distribution of porosity could be simulated as a damped oscillatory function of the distance from the wall. Georgiads et al.[4, 5, 6] studied the unidirectional transport phenomena of flow and heat transfer in porous medium with the stochastic model and indicated that the results obtained by the random porosity were larger than those ottained by the Forchheimer model with the constant porosity model. Saito et al. [7] studied the effects of the porosity and void distributions on the permeability by Direct Simulation Monte Carlo method. The results showed that the characteristics of the porosity distribution were disordered and random. Recently, Fu and Huang [8] developed a random porosity model to investigate the forced convection in the porous medium. The variations of distributions of heat transfer rate on the heated wall were no longer smooth and varied apparently. As for the results of the natural convection accompanying with evaporation in porous medium under the random porosity model are seldom investigated in detail. The aim of this study i~ to use a random porosity model to investigate double diffusive natural convection phenomena in an enclosure filled with porous medium numerically. The distribution of the random porosity is generated by Kinderman-Ramage procedure [9] to follow the normal (Gaussian) distribution criteria. The Darcy-Brinkman-Forchheimer model of the porous medium is taken into consideration. The SIMPLEC numerical method is adopted to solve the governing equations with the iterative computing processes.

aaal.v~

A physical model of a two-dimensional rectangular porous medium enclosure is shown in FIG. 1. The random porosity model proposed by Fu and Huang [8] is adopted, and the distributions of the porosity are disordered and random.

(3)

Vol. 27, No. 1 NATURAL CONVECTION IN A POROUS ENCLOSURE 121

H Aq~

0T=0 ~---=0,u=v=0 ay 'ay

Extreme thin film

~

qt(Y) t~qs(Y) TL(y ) 1 CL(y) tlL=O v~ = 0 Porous Medium rjt = 293K ~R = 30% tl~v=O ~ = 0 ac=0 . . . 0 ~ay 'ay ~i _{r I} L FIG. 1

Physical model and porous medium with the random porosity distribution.

The width and length o f the enclosure are H and L (=H), respectively. Two horizontal walls of the enclosure are adiabatic and the gravity is downward. The left wall is subject to a uniform heat flux qw and low temperature, concentration, and relative humidity conditions of the fight wall are

T R = 293K, C a and ~ R = 30%, respectively. In order to enhance heat transfer rate of the left wall, the left wall surface is wetted with water film in this study. The thickness of the water film is assumed to be extremely thin and evaporated immediately on the wall, then the water film neither penetrates into the porous medium region nor flows downward. The water film is evaporated by the heat flux imposed on the left wall, and the condition of relative humidity of the left wall surface is saturated. Consequently, the sensible and latent heat transfers of natural convection occur in the enclosure simultaneously, and the temperature TL(y ) and concentration CL(y) on the left wall are no longer regarded as presumably constant in advance and affect mutually. The working fluid in the enclosure becomes air-moist fluid.

For facilitating the analysis, some following assumptions are made. (1). The porous medium is made

of non-deformable pure copper spherical beads ( k s = 386 W m l o C-I ) which are not chemically reactive

with the fluid. (2). The flow in the enclosure is laminar, steady and two-dimensional. Corresponding to the thermal and concentration boundary conditions selected, the condensation does not occur in the

(4)

enclosure. (3). The effects of Soret and Dufour induced by mass transfer are neglected, and Boussinesq approximation is held. (4). Except the left wall region, the properties of the air-moist fluid medium are constant and based on the temperature of the right wall. The properties of the air-moist fluid on the left wall which are calculated by [10] are based on the temperature of the left wall under a saturated condition. (5) The form of the effective thermal conductivity l%fr are defined as kefr / kf = 4 ln(k s / kf ) - 11 which is proposed by Shonnard and Whitaker [ 11 ] for high thermal conductivity ratio. (6). The permeability K =

g3d2p/[150(1- s)2] and inertia factor F = 1 . 7 5 / ( ~ 5) refer to [12]. Based on the above assumptions

and with the following characteristic scales of H, T R , q w, v R, k ~ , the governing equations adopted in [13], geometry dimensions and boundary conditions are normalized as follows. The porosity s in the following equations is random and obtained from [8]. In the meantime Darcy-Brinkman-Forchheimer model is used in the momentum equations.

Governing equations: ~U ~V - - + - - = 0 ( 1 ) 8X 0Y

ffX-t,-~)

~ - I , - ~ - ) = - ~ - t,~--~T ~--~T)

Da

- ~ a

sU

(2)

U 0 C V ~ + v c~ ( V~ _¢gP+["~92V

02V'~

1 sV

FIOI

(3)

1 ¢

a 0)

U O0 + V 6~0 = P r (~--~- + ~ - ) c~X 0Y (4) u O W + v C3w

- - / ~

=

f"O2W

_

l

O2W'~ 0 x o v -- ScC 0x~ + 0 y ~ ) (5)

Where' X = ~ ' Y = y H ' U =--,UHvR V = v H , P - p H 2 v R - - - o r 2 , 0 - kes (T - T R ) q w H , W - (eL (y) _ CR ) ( C - CR) ,

_ g~tqw H4 g[3c(C - C R)H 3 _ vlt v R _ keff K

,Pr - - - , Sc = - - , ct

Da : I01

=,/~

+ v2

Ca- t

Grc

' V2 ff'R DR -- ~ - P ' ' ' kerr 1 M a [St = ~ , fie = (-:'7--= - 1) (6) rvl~

(5)

Boundary conditions:

~0 aW = 0 (7)

On OP (top wall) U = 0,V = 0, ~ = 0, 0~f

On PQ (right wall) UR = 0,VR = 0, Or( = 0, W R = 0 (8) 0 aW

On QR (bottom wall) U = 0,V = 0, ~-~ = , ~ - = 0 (9) O n R O (left w a l l ) U L =0, V L =0, 0 L =0L(Y), W L = W L ( Y ) (10)

The other important values of evaporation mass flow rate riaL(y ) , sensible heat flux qs (Y) and latent heat flux q e (Y) are calculated by the following equations, respectively.

pLDL 0 C = thE(Y)= (l_CL(y))0X[x=0 pLDL (CE(XI,Y)-CL(Y)) (1 -CL(Y)) _ldx 1 2 (11) q w = q s ( y ) + q e ( y ) (12) q s ( Y ) = - k e f f - ~ x = 0 = - k e f f (TE(Xl'Y)-TL(y))

~dx~

2 (13) PLDLhfg ~ x = 0 = pLDLhfg ( C E ( X l ' y ) - C L ( Y ) ) (14) qe(Y) = rhL(Y)hfg = (I_CL(Y)) (1-CL(Y)) l d x 1

2

Based on the experimental results of Fu and Huang [8], in this study the mean porosity ~ is conveniently regarded as 0.5 and the standard deviation tse is 0.05 which is 10% of the mean porosity.

h

The mean diameter dp is selected as a characteristic bead diameter and equal to 1.0 mm. For the necessity of computing procedures, the theoretic form of the porosity distribution of the random porosity model is generated by the Kinderman-Ramage procedure in [8] and the distribution of the general random variable e is regarded as the porosity distribution of the random porosity model in this study.

~Ialae, lJ.catm~aml

The SIMPLEC algorithm [14] with TDMA solver [15] is used to solve the governing equations (1)- (4). Under-relaxation factors are 0.1-4).5 for all the variables of velocity, temperature, and concentration in the cases, respectively. The staggered mesh is used, the finer meshes are set near the solid wall regions

(6)

and the meshes are expanded outward from the boundary wall with a scale ratio 1.03. The conservation residues of the governing equations and the relative errors of all variables are used to examine the convergence criteria defined as follows :

2 /2

[Residue of ¢ equation c v 1 -<10-4' ~ = U, V, 0 and W. (15)

L

Cn+' ¢°V

- / i q ~ , . i i _< 1 0 - 3 , ~ = U , V , P , 0, W a n d q . (16)

The accuray of the similar numerical method is validated in Fu et al. [16, 17, 18]. The comparison between the results obtained from the present methos and those from Goyeau et al. [19] are indicated in Table 1. for the conditions of double diffusive natural convection in a porous cavity with aspect ratio A=I and Darcy number Da=10 -7 . The deviations between these two results are less than 40/60.

TABLE 1

The comparison of the results of [20] under the constant porosity model with the present study.

Ra* Le Nu (Goyeau et al.) 50 1 1.98 10 100 100 1 10 100 1.98 1.9984 1.98 1.9984 3.11 3.1558 3.11 3.1558 3.11 3.1558 Nu Error (present (%) study) 1.9984 0.92 0.92 0.92 1.45 1.45 1.45 Sh (Goyeau et al.) 1.98 8.79 8.9514 27.97 28.754 3.11 3.1558 13.25 13.660 41.53 42.358 Sh Error (present (%)

study)

1.9984 0.92 1.8 2.73 1.45 3.00 1.95

Results and Discussion

The distributions of porosity are generated by the random porosity model [8] with a given mean

porosity g and a standard deviationo~. However, it is difficult to solve all of the patterns. Therefore, in

the present study only ten patterns with each mean porosity (~ =0.5, 0.7) are presented to investigate the effects of the random porosity model on the flow, thermal and concentration fields. The parameters of two

(7)

mean porosities, three standard deviations and eight random variables tabulated in TABLE 2 are adopted. The results of the grid tests for computation are listed in TABLE 3 and the grids of 74× 74 are selected for the following calculation processes with Ax I = Ay I = 0.0075 and A x i + l / A x i = l . 0 3 which are symmetric to central lines of the computing domain. In this study, the range of Rayleigh numbers Ra

[=(gl3tqwH4)/(kegVR~tR)] from 1.839x106 to 2 . 9 6 3 x 1 0 6 , the Prandtl number Pr is 0.02583 and Schmidt number Sc is 0.5425.

TABLE 2

The main porosity ~, standard deviation o and random variable ~ adopted for the cases.

Case ~ ~c (%) ~ Case 8 o'~ (%)

(Ra= 3.678 x 106 ) RUN1 0.5 0.7 10 ~ RUN6 0.5 0.7 10 ~6 RUN2 0.5 0.7 15 ~l RUN7 0.5 0.7 10 ~7 RUN3 0.5 0.7 5 ~l RUN8 0.5 0.7 10 ~8 RUN4 0.5 0.7 I0 ~4 RUN9 0.5 0.7 I0 ~9 RUN5 0.5 0.7 10 ~5 RUN10 0.5 0.7 10 ~10

Grid tests forRa = 3.678 Ax1, AYl lRatio(AXi+l__

₎

/xX i 0.005 0.005 0.0075 0.075 1.03 1.05 1.03 1.05 Grids " N---u ( N x x N y ) 94 × 94 1.863868 0.0075 74 × 74 1.864367 0.01 74 × 74 1.864232 0.15 60 × 60 1.864193 0.15 TABLE 3 x 106,8=0.5, o'r,=lO% and Axl,AYl

Op--o.oo5.

Ratio( Ax i+l Ax i 1.1 1.03 1.03 1.05 Grids " N--u " ) ( N x x N y ) 42 x 42 1.865727 62 × 62 1.862837 46x 46 1.863325 40 x 40 1.864103

There are two selected cases (RUNs 2 and 3) shown in FIGs. 2(a)-(d). The global porosity distribution maps and the near wall local porosity e v distributions along the Y direction at X=0.00375 are illustrated. In the global porosity distribution maps, the total area are divided into several main porosity regions with different grey scale, and the darker scale represents the smaller porosity. The

(8)

porosity distributions are not in order, then sparse and dense structures are generated at random. The

variations of porosity in most region are from ~ - ~c to ~ + c r which are consistent with the results of

normal distribution. As the standard deviation is smaller, the variation of the distribution of porosity is smaller, too. 08 ~ 070 j ~ i 0.6 0,60 ~ , v 0.50 y 0.5 ~ t 04 I °'~° ~

~ ~ 1

0.30 ! 0.2 =0,20 0 . 0 I ~ i 0.0 0.3 0.5 0.7 O0 0.2 0 4 0.6 0 8 1.0

x (a) RUN2 ~;Y

I

,

1

~057 0.53 Y 0.49 y 0 . 5 0.45 041 0.0 ~ - - 0.0 0.3 0.5 0.7 0.0 0.2 0 4 0,6 0.8 I 0 x (c) RUN3 ~ Y (b) RUN2 (d) RUN3 FIG. 2

The global porosity distribution maps (a) and (c) and the local porosity e y distributions (b) and (d) at X=0.00375 along the Y direction under R a = 3.678 x 106 and ~ = 0.5.

The distributions of streamlines, isotherms and isosolutal lines are illustrated in FIG. 3. For the same reason, as the standard deviation is larger, the distribution or deformation of the above different lines becomes apparent. The isothermal lines are almost vertical as approaching the right wall, which means the heat conduction mode to be dominant. Since the distributions of porosity are random, the streamlines are no longer smooth.

(9)

Vol. 27, No. 1 NATURAL CONVECTION 1N A POROUS ENCLOSURE 127

In TABLEs 4 and 5, the local Nusselt number Nu(Y) at different locations of Y-axis for ~=0.5

and ~ --0.7 are presented, respectively.

TABLE 4

The local Nusselt number at different locations of Y-axis under Ra = 3.678 x l 06 and =0.5. ( RUN11 : constant poros[ty mod d )

I I

NaVY ~

I,]RUNI RUN2 RUN3 RUN4 RUN5 RUN6 RUN7

RUNSIRUN91

RUN10 RUNll

Y=0.00375 i 1.942 1.950 1.937 1.943 1.941 1.942 1.944 1.9431 1.944i 1 . 9 4 1 1.943 Y=0.2083 1.922 1.928 1.919 1.923 1.922 1.923 1.923 1.9221 1 . 9 2 4 1 . 9 2 2 1.923 Y=0.4034 1.885 1.$88 1.883 1.886 1.884 1.885 1.885 1.8841 1.8861 1 . 8 8 5 1.885 Y=0.5966 i 1.843 1.841 1.845 1.843 1.844 1.843 1.843 1.8441 1.843 i 1 . 8 4 4 1.843 Y=0.8050 1.807:1.802 1.81(1 1.80,* 1.806 1.806 1.'806 1.8071 1 . 8 0 5 1 . 8 0 7 11805 Y--0.9963 1.792 ~ 1.786 1.795 1.785 1.790 1.7911 1.790 1.791] 1 . 7 8 9 1.79(] 1.789 TABLE 5

The local Nusselt number at different locations of Y-axis under Ra = 3.678 x 106 and e =0.7.

N~Y) ;RUN1 Y=0.00375 2.426 Y=0.2083 , 2.306 Y=0.4034 2.114 Y=0.5966 1.903 Y=0.8050 1.735 Y=0.9963 1.665

RUN2 IRUN3 IRUN4 [RUN5 IRUN6 IRUN7 IRUN8 IRUN9 I~ ~II] IP4II~J~ kl I P I [ I ~ | Ill'J(1/;11~Iti I~'~ IP J I I ~ ] IP4g[ll I IPJII{I !,

RUN10 1RUNll 2.4231 2.363

2.301 2.25i

2.1041 2.074

1.9161

1.891 1.7411 1.722

116531 1.65(1

The local Nusselt number is defined as follows.

Nu(Y)= ht(y)H = N u s ( Y ) + N u t ( Y )

ke~ qw

,where ht(y ) - (17)

_ qs 1 + qe 1 ( T L ( y ) - T g )

qw 0L(Y) qw 0L(Y)

Due to the random distribution of porosity, even at the same locations the values of the local Nusselt number of each run are slightly different. The average local Nusselt number of the RUN1 and RUN4-10

(10)

numbers is about 3.11% which is smaller than that ( t ~ =10%) of porosity distributions of the RUN1 and RUN4-10. The effect of the standard deviation of porosity distribution on the standard deviation of the average local Nusselt numbers is not very apparent.

Streamlines Isothermal lines Isosolutal lines

-' ~- -;,,,/!, , . . . . , , ~ ~ - ',I, . . . . ' ' i " ' l', i,~i ; ~ ' ";',~ll , ,' ,, I : a )

i i i , !

~b)

FIG. 3

The distributions of streamlines, isothermal lines and constant concentration lines. (a) RUN2 (b) RUN3 for Ra = 3.678 x 106 and e =0.7.

In FIGs. 4(a) and (b), there are the distributions of local sensible and latent Nusselt numbers, N u s ( Y ) and N u e ( Y ) , respectively. Similar to the results shown in TABLE 5, the local sensible and latent Nusselt numbers are different for each case. The cold fluids reach the left (hot) wall at the lower region first, then the strong convection enhances the latent heat transfer and t h e N u e ( Y ) is larger t h a n N u s ( Y ) in the lower region. Oppositely, in the upper region, the value o f t h e N u s ( Y ) is larger than

that of the N u e ( Y ) - In general, the larger the g is, the stronger the evaporation is. As a result the local Nusselt numbers on the left wall becomes larger. Besides, the larger the standard deviation ~ (RUN2), the variation of the porosity distribution becomes more apparent, then the fluctuations of the local Nusselt numbers of N u e ( Y ) and N u s ( Y ) along the left wall increase as the standard deviation c~ increases.

(11)

Vol. 27, No. 1 N A T U R A L CONVECTION IN A POROUS ENCLOSURE 129 1 . 4 , . . . ~

1.o f / _ -

. ~ , . . . .

. .

0 6 ' , " 0.4 ~ ~ 0.0 0.2 0.4 0.6 0.8 1.0 V ( a ) N u e ( Y ) 2.5 ~ r i r ~ 1.5 gt~4 1.0 0 • 5 0.0 0.2 0.4 0.6 0.8 Y 1.0

(b)

FIG. 4

The dislributions o f ( a ) local sensible Nusselt number N u s ( Y ) , (b) local latent Nusselt number N u t ( Y ) o f ~ =0.7 for the selected cases (RUNs 1-4 ) on the left wall for R a = 3 . 6 7 8 x 106 .

R~ 1.0

f l ; s

r r

0.8 0.2 F" ... " . . . . . . . _ _ _ . . . : 0.6 ! . . . i 0.4 0.4 0.6 R s 02I J ~ • ~ R~ 0.8 0 0 ~ , ~ ~ _ ~ _ ' 0 1 0 ~ 0.2-- 0.4"- 0.6 0.8 1.6 .0 1.0 0.8 0.6 Rt 0.4 0.2 0.0 ~ , i , r , i , 0 . 0 0.0 0.2 0.4 0.6 Rs 0.8 1.0 0.2 0.4 0.6 0.8 1.0 Y v (a) (b) FIG. 5

The distributions o f R e and R s on the left wall for different Ra situations (a)for ~ =0.5 and a ~ = 1 0 % in RUN1. (b)for ~ =0.7 and 6~ =10% in RUN1.

. . . R a = 9 . 1 9 5 × 1 0 6 , - . . - R a = 7 . 3 5 6 x 1 0 6 , _ . - . R a = 5 . 5 1 7 × 1 0 6 , - - - Ra = 3.678 × 1 0 6 , - - Ra = 1.839 × 106 and e e e Ra = 3.678× 106(~)

Shown in FIG. 5, the values o f R e and R s indicate the ratios o f latent heat flux qe and sensible heat flux q s to the total heat flux q w - The value o f the ratio R means the contribution o f the above individual heat flux to the total heat flux. The dot symbol indicates the results o f the constant porosity model. A s the value o f Ra is smaller, the driving force o f fluid flow becomes small which results in

(12)

130 W.-S. Fu and W.-W. Ke Vol. 27, No. I

weakening of the convection and the evaporation simultaneously. The value of R e is then smaller under the smaller Ra case shown in FIG. 5(a) . Oppositely, as the value of Ra is larger, based on the reason mentioned above, the evaporation becomes dominant and the value of R e becomes larger. Besides, as the ~ is larger, which means the variation of porosity to be more drastic, the fluctuations of the distribution of R e and R s become more apparent as shown in FIG. 5(b).

A double diffusive natural convection in a square enclosure filled with porous medium under evaporation situation is investigated with the random porosity model numerically. The effects of Rayleigh

number Ra, mean porosity e and deviation ~e on heat and mass transfer are examined and can be

summarized as follows:

(1). The standard deviation is larger, the variation of the distribution of porosity is more chaotic which causes the flow pattern to be distorted apparently.

(2). The effect of the random porosity model on the distributions of local Nusselt number is remarkable. (3). The latent heat flux plays an important role in heat transfer mechanism under high Rayleigh number.

As the mean porosity is larger, the contribution of latent heat flux becomes more important.

~lame, n d a t l l ~

C mass or concentration fraction Cp specific heat of fluid

dp mean porous bead diameter D binary diffusion coefficient Da darcy number

Dp dimensionless bead diameter, dp / H F inertial factor; Forchheimer factor g gravitational acceleration

Gr¢ Grashof number for mass diffusion Gr t Grashofnumber for thermaldiffusion h c concentration transfer coefficient along the

(13)

hfg vaporization enthalpy

h t thermal heat transfer coefficient along the vertical wall

H dimensional width

keff effective thermal conductivity k f thermal conductivity of the fluid k s thermal conductivity of solid phase K permeability

L dimensional length rfi evaporation mass flow rate M molecular weight

Nu local Nusselt number along the vertical wall N u mean Nusselt number

p dimensional pressure P dimensionless pressure Pr Prandtl number, v R/ct R q heat flux Ra Rayleighnumber, (gl3tqwH4)/(k~frvRctR) Ra ° Darcy-Rayleigh number, (gl3tqwH2K)/(k~ VR~t R )

R e ratio of latent heat flux in the left wall to total heat flux in the left wall, q ~ / q w R s ratio of sensible heat flux in the left wall to

total heat flux in the left wal q s / q w Sc Schmidt number, v ~ / D R

T temperature u,v dimensional velocity U,V dimensionless velocity

x, y dimensional Cartesian coordinate X, Y dimensionless Cartesian coordinate

GREEK SYMBOLS a thermal diffusivity

[3 t coefficient o f volumetric expansion with temperature

13c coefficient of volumetric expansion with mass fi:action

porosity, e = ~ + ~ x o~

e v near wall local porosity at X=0.00375 computational variable

~t viscosity

v kinematic viscosity random variable

0 dimensionless temperature

0* modified dimension temperature, 0/kee r p fluid density

eye standard deviation relative humidity

u/ dimensionless stream function,

u = 0 v / 0 Y

SUBSCRIPTS a air

C.V. control volume eft effective value E enclosure f fluid phase i index latent heat L left wall p porous medium R fight wall s solid phase S sensible heat v vapor W solid wall SUPERSCRIPTS n the nth iteration index - - mean value

--~ velocity vector OTHER

(14)

1. L.H.S. Roblee, R. M. Baird and J. W. Tiern, American Institute o f Chemical Engineers Journal, 4, 460 (1958).

R. F. Benenati and C. B. Brosilow, A. I. Chem. Eng. J., 8, 359 (1962).

P. Cheng, A. Chowdhury and C. T. Hsu, Convective Heat Mass Transfer in Porous Media, ed. S. Kaka~ et al., 625, (1991).

J. G. Georgiads and I. Catton, The Physics o f Fluids, 30, 1017 (1987).

I. Catton, J. G. Georgiads and E Adnani, Proceedings o f 1988 National Heat Transfer Conference, 1, ASME HTD-96, Houston, 767 (1988).

J. G. Georgiads, Convective HeatMass Transfer in Porous Media, 499 (1991).

A. Saito, S. Okawa, T. Suzuki and H. Maeda, Proceedings o f ASME/JSME Thermal Engineering Joint Conference, 3, 297 (1995).

W. S. Fu and H. C. Huang, Int. J. Heat andMass Transfer, 42, 13 (1999). J. Jr. Kennedy and J. E. Gentle, Statistical Computing, New York (1980).

T. Fujii, Y. Kato and K. Mihara,, Sei San Ka Gaku Kenkyu Jo, 66, Kyu Shu Dai Gaku, Kyu Shu, Japan (1977).

11. D.R. Shonnard and S. Whitaker, Int. J. Heat andMass Transfer, 23, 503 (1989). 12. K. Vafai, Journal o f Fluid Mechanics, 147, 233 (1984).

13. S.W. Hsiao, E Cheng and C. K. Chen, Int. J. Heat andMass Transfer, 35, 3407 (1992). 14. J. E Van Doormaal and G. D. Raithby, Numerical Heat Transfer, 7, 147 (1984). 15. S.V. Patankar, Numerical Heat Transfer and Fluid Flows, Hemisphere (1980). 16. W. S. Fu and H. C. Huang, Int. ,~ Heat andMass Transfer, 40, 2261 (1997).

17. W.S. Fu and H. C. Huang and W. Y. Liou, Int. J. Heat andMass Transfer, 39, 2165 (1996). 18. W. S. Fu and W. W. Ke, J. Chinese Institute of Engineers,

22,

325 (1999).

19. B. Goyeau, J. P. Songbe andD. Gobin, Int. J. Heat and Mass Transfer, 39, 1363 (1996).

2. 3. 4, 5. 6. 7. 8. 9. 10. Received October 28, 1999