Int. Comm. HeatMass Transfer, Vol. 24, No. 4, pp. 543552, 1997 Copyright © 1997 Elsevier Science Ltd Printed in the USA. All rights reserved 07351933/97 $17.00 + .00
Pergamon
PII S07351933(97)000390
FREE CONVECTION FROM A VERTICAL PLATE IN A POROUS MEDIA SUBJECTED TO A SUDDEN CHANGE IN SURFACE TEMPERATURE
S.D. Harris and D.B. Ingham Department of Applied Mathematical Studies University of Leeds, Leeds, LS2 9JT, England I. Pop Faculty of Mathematics University of Cluj, R3400 Cluj, CP 253, Romania
(Communicated by P.J. Heggs) ABSTRACT The unsteady heat transfer process involved in free convection flow along a vertical surface embedded in a porous medium is investigated. An analytical solution has been obtained for the temperature/velocity field for small times in which the transport effects are confined within an inner layer adjacent to the plate. Then, a numerical solution of the full boundarylayer equation is obtained for the whole transient from the initial unsteady state to the final steady state. Detailed results of the effect of the temperature inputs on the transient process are given. © 1997FAscvicrScienceLtd Introduction Convective heat transfer process is of fundamental importance in a variety of practical applications, ranging from mechanical engineering to geophysics and recent reviews by Nield and Bejan [1] and Nakayama [2] give the extent of the research information on this area. The unsteady convective heat transfer problems encountered in these applications are rather complex and they can be solved either analytically or numerically. Numerical techniques, such as finite differences or boundary elements, are commonly used for complex problems, while analytical methods leading to exact solutions are preferred for their simplicity in engineering applications. In spite of extensive effort to analyse the transient process during cooling or heating of a surface which is embedded in a porous medium, it appears that the literature is lacking of simple solutions which determine the heat transfer characteristics from such surfaces. In this paper, we present a method for determining the heat transfer quantities for 543
544
S.D. Harris, D.B. Ingham and I. Pop
Vol. 24, No. 4
a vertical surface embedded in a fluidsaturated porous medium which is subjected to an impulsive change of the temperature of the plate. Thus, it is assumed that for time ~ < 0 a steady temperature/velocity has been attained in the boundarylayer which occurs due to a uniform temperature T1 of the plate. Then at time ~ = 0 the temperature of the plate is suddenly changed to T2 and maintained at this value for ~ > 0. The analytical and numerical results show that the transient process is strongly affected by the levels of the existing energy inputs. Basic Equations Consider a vertical flat plate embedded in a porous medium which is at a constant temperature T~ and the plate is maintained at a uniform temperature T1. At time 9 = 0 the temperature of the plate is suddenly changed to T2 and is maintained at this constant value for ~ > 0. Using the DarcyBoussinesq and boundarylayer approximations, the conservation equations for the unsteady free convective flow are
Ou
u =
OT
+
Ov
913K //
OT
+
=0
(1)
(T  Too)
OT
=
(2)
O2T
(3)
which have to be solved subject to the boundary conditions
u(x, oo, 9) : O,
T(x, oc, ~) : Too
u(0, y,¢) : v(0, y,¢) = 0,
v(x,O,~)=O,
]
T(0, y, 9) : T~o /
(4)
T(x,O,9)= T2
for "~ > 0 and 0 < x,y < o0. Here u and v denote the velocity components along the x and y directions, with x being measured along the plate starting at the leading edge and y measured normal to it, T is the fluid temperature and the other quantities are defined in the Nomenclature. We shall now proceed to transform Equations (1)(3). For ~ > 0 the nondimensional, reduced streamfunction, f, and the temperature, 0, are defined as
TToo ¢ = Vc~(:c) f ( ~ , r ) , where/kT1 = T1  Too and
0(7,~) 
AT1
(5)
Vol. 24, No. 4
, 
Here
Y
~(~),
C O N V E C T I OFROM N A PLATE IN POROUS MEDIA
(2\½ 6(x)=x(,~a~)
6(x) is the
,
"r
ae
o~ U¢=Rax.
~ ( 6 ( ~ ) ) 2'
•
Ra.
545
gfl K ATxx ~
(6)
boundarylayer thickness at ~" = 0 and ¢ is the streamfunction which is
defined by u = ~ and v =  ~ .
Substitution of the expressions (5) into Equations (1)(3)
yields the following equation for f:
03f + (_1+ 2.r~) 002f ( Of) ~02f ~+ f2~'~ =0
Or?~
(7)
where ~ = 8. Equation (7) has now to be solved for ~" > 0, subject to the boundary conditions
f(o,~) = o, where AT2 =
AT2 ~(o~, Of ~) 0
(o, ~) = h~'
(8)
T2  T¢~. For the steady boundarylayer at ~" = 0 one can write f(~?, 0) = f0(~?),
so that, from Equation (7), f0(~) satisfies the ordinary differential equation
]o" + fof~' =
0
(9)
along with the boundary conditions /o(0) = O,
f~(O) = I,
fo(OO) = 0
(10)
where primes denote differentiation with respect to 7/Small Time Solution, ~ << 1 In this case there exists an inner boundarylayer and in order to study this layer it is convenient to use the new variables
f(~,.) =.½F(~,T),
~ = ~~2~~
(11)
as in [3, 4]. Equation (7) then becomes
io~~+~  l + , N
o2g+
~
N/5~=o
(12)
which has to be solved subject to the boundary conditions F(0,~) = 0,
OF
~(0,r)=
AT2
2~~
(13)
The solution in this growing inner layer is taken to match the outer steady boundarylayer, which at small 7/can be approximated by the series expansion
546
S.D. Harris, D.B. Ingham and I. Pop
1
f0(~) ~ 77 + ~a~
2
1 4 ~a'7 

Vol. 24, No. 4
1 2s ZT6 ~ 7/
+ 0(76)
(14)
where a = f~'(0) = 0.62756, see [5]. The substitution of expression (11) into Equation (14) yields, for large values of ~,
F(,L',) ~ 2~ + 2 a ~
l

2 4 3 ~,~ ,~  4 ~,~.,_~ + o(A)
(15)
The behaviour of the inner layer as ~ ~ oo is to be matched with the steady outer solution (15). Therefore, the solution of Equation (12) within the inner layer results as follows:
F ( ( , T ) = Fo(() + ~½FI(~) ÷ "r~F2(~) + v2F3(() ÷ O(~~) The functions F~(~),
(16)
i = 0, 1 , 2 , . . . , can easily be obtained from Equations (12), (13) and
(15). The resulting expression for of is given by ~=22
4a~'½
o{
(1~)
÷
[2;2 1 ~~2
÷

÷
(17)
where erfc is the complementary error function. Under the transformation (11) the resulting velocity/temperature function is given at small times T by
Also, the nondimensional heat flux from the plate can be calculated from the expression
q~(~)=
_
0'7 ~1~0= ~ +

h~J
~"
~

~
~
+ o(J)
(19)
Numerical Solution The governing partial differential Equation (7), or its equivalent form (12), are parabolic and can be integrated numerically using a stepbystep method similar to that described by Merkin [6], provided that the coefficient of ~
or °~f remain positive. This marching
method gives a complete solution for ~ < ~*, where T* is the m a x i m u m value of ~ reached in the numerical scheme, which is less than the exact time v = ~T~" aT~ The matching of the solution at ~  ~* to the asymptotic steady state solution may now be achieved using a variation of the method first described by Dennis [7]. In order to accurately evaluate the nondimensional heat flux from the plate initially we apply the stepbystep scheme to Equation (12) and begin the numerical solution at the
Vol. 24, No. 4
CONVECTION FROM A PLATE IN POROUS MEDIA
547
small time r = r0 using the profile given by expression (17). The evolution of the function @ = ~ is governed by the integrodifferentiai equation
r(
[email protected])
[email protected]~r
[email protected][l{
[email protected] ({ ,r)d{'lj
7~~
(9.0)
which has to be solved subject to the boundary conditions =
'a&,')
= 9.
where the undisturbed state f~ is enforced at a large value of ~, denoted by ~oo. The stepbystep procedure of Merkin [6] is now applied to Equation (9.0) for ~(~, r) in precisely the form described by Harris et S~I,j+
al. [4]. Thus
the finitedifference equation
}  9.S/{,,,4+1 [. S~_I,j+~. it
<,,')=
_
,
 <,,+,)[,,
(9.2)
_o
represents an approximation to Equation (20) at ~ = (i  1)h ~ and r =
1 rj + 51Xr5, where
i1
S,~,j+½ = ¢~,J+~ + @~,J,
ei,i=
i,j+½  ~ t,°~,j+½
'~ ~ /
i,=2
'~*~
i1
g((Ih,./
[email protected])+Y]@ej,
A=
(9.3)
+1
i*=2
h~ = ~
and the boundary conditions require that AT2
(24)
This system of nonlinear algebraic equations is solved iteratively using the NewtonRaphson method and by employing the method proposed by Doolittle to decompose the associated Jacobian matrix into the product of a lowertriangular and an uppertriangular matrix. This iterative process is repeated until the absolute difference between successive approximations reaches a value less than some tolerance q . The initial time increment At0 at time r = r0 is set to some prescribed small value and a time step doubling procedure is adopted. Each time step is covered using first one and then two time increments. If the absolute difference between the two solutions is less than a preassigned tolerance e2 then the time step is doubled. The size of the discretised ~space under consideration, 0 _< ~ _< ~oo, increases with time. If we regard • r/oo to correspond to ~/ = oo then at the nearest time "~ below r  ( ~ ) ~
it
is necessary to transfer to a version of the stepbystep procedure which uses a constant mesh width, provided that this time @ is less than r = ~2AT2
or,
equivalently, R = ~AT1 < 2 (~)~ ~
"
In such cases we retain the current parameter values for the time increment and tolerances
548
S.D. Harris, D.B. Ingham and I. Pop
Vol. 24, No. 4
and apply the stepbystep method to the nondimensional temperature function 0 = satisfying the integrodifferential equation (1
2,tO)
o,rO0 0200772 [" .( + ~00Jo
2rr~_~r ~00( ,,,r))
(25)
and subject to the boundary and initial conditions
O(0,,r)
1
AT2 AT1'
007°°'z) = 0,
O0?,'?) =
2¢(~'~)
(26)
~= ~r 2¢ 2
where the initial profile for 77 > 2~½ remains the undisturbed state. The finitedifference equation
+(h") 2 (S:"+I,j+½ ~'x,y+½)[¼(1 2X)~2:~+½ + )~O:,i] = 0 can then be derived to approximate Equation (25), where S1j+½
aTe, SN+I,j+½ = 0,
h" =~6"°° and the remaining terms are defined analogously to expressions (23). The resulting nonlinear system of equations is then solved as described for Equation (22). The matching of the steady state similarity solution 0(,7, oo) = Rf0(~R~ ' 1)
(28)
w h e r e R = a TzT at large times with that which is valid at z = ,r * is now achieved using an
adaptation of the method of Dennis [7]. The system of equations
Of
0,
020
00
00
(29)
where
of
p(~,'r) = f  2TTN,
q(~,T) = 1  2,r0
(30)
must now be solved subject to the boundary conditions 001,T*)=0"(~7) , f(0,,r) = 0,
l
!
0(rl,,r~)=_Rf;(~/R~), 0(~oo,,r)
0,
1
1
f(r/,~'~)R~f0(T/R~)~
0(0,,r) = R,
,r* < ,r < ,r~
/
(31)
where Too is some large but finite time corresponding to "r = oo. By replacing 7? derivatives by centraldifferences and oe by either a backward or forward difference depending on whether q(r/,,r) > 0 or q(r/,r) < 0, respectively, Equation (29) becomes (l+2hPij)0i+lj+(11
0
(2 + ~~2'q~J']~ 0ij,
: 
(32)
where 0i*J = 0ij+l if qij < 0 and 0~,j = 0i,j~ if qi,j > 0. The iterative solution of the system
(32)throughout the domain isachieved using precisely the procedure described by Harris et aL [41.
Vol. 24, No. 4
CONVECTION FROM A PLATE 1N POROUS MEDIA
549
Results and Conclusions The restriction to finite dimensional ( and 7 spaces was achieved by taking ( ~ = 10 and 7~ "~ 10, respectively, where the precise value of 7~ = 2x/~(oo is dependent on the final t i m e ~ reached in the first stepbystep solution. The effect on the numerical schemes of variations from these values of ( ~ and 7¢o, while keeping h ~ and h" constant, respectively, was investigated and it was concluded that any larger values of (oo and 7oo produced results which were indistinguishable from those presented in the figures. Thus the application of the stepbystep m e t h o d in 7 and r variables is only required when R < 2. The values of the tolerances el and e~ as average errors over the unknown grid points were taken to be el = 10 4 and e2 = 10 ~. The initial time To and first time increment A~0 were assigned the values TO = 5 × 10 .8 and AT0 = 10 6. More restrictive values of these parameters were considered and discovered to produce numerical results which did not show any significant variation. The m a i n source of deviation in the solutions for the fluid t e m p e r a t u r e arises by considering changes in the number of grid spaces N ~ and N ". It was observed that as N ~ and N ~ increased, and consequently h ~ and h" decreased, the initial development of the numerical solution for the nondimensional heat flux from the plate approached that of the small time solution. Figure 1 shows the variation of the profile of 0(7 , T) at various times T calculated using h ~  0.0125, h ' = 0.025, N ~ = 800 and N" = 400 for ratios of surface temperatures, R = 0.5 and 2, where refinements in the spatial mesh produced an insignificant improvement in the accuracy of the solution. The evolution of the surface heat flux
q~,(T) with time r
is illustrated
in Figure 2 for R = 0.5 and 2. The numerical, transient solution is shown to develop closely following the small t i m e solution (19) and is graphically almost identical when ~ < 0.9 and r < 0.2 for R = 0.5 and 2, respectively. The matching m e t h o d solution over the finite region T* < ~< too is achieved by enforcing the steady state solution to apply at Too = 12 and the convergence criterion was set by assigning the value e3 = 5 x 10 .9 for the average absolute error in 0. No significant improvement in accuracy was achieved by either increasing Too or reducing e3. The o p t i m u m value of the relaxation p a r a m e t e r was found to be w = 1.5. Solutions for 0(7 , r ) were found for the grid spacings h ~ 0.025, k ~ 0.05; h ~ 0.025, k ~ 0.025 and, for R = 2, h ~ 0.0125, k ~ 0.025. As the mesh size was reduced, and consequently the computational time increased, the numerical solution was observed to change only marginally over the solution domain. Thus h ~ 0.025 and k ~ 0.05 were used so t h a t n = 400 and m = 221 and 235 for R = 0.5 and 2, respectively.
550
S.D. Harris, D.B. Ingham and I. Pop
Vol. 24, No. 4
The evolution of the heat flux at the surface q~(z) from ~ = T* to ~ = ~oo is displayed in Figure 2 and shown to proceed past local m i n i m u m and m a x i m u m values near ~" = T* to the steady state solution at large times for R = 0.5 and 2, respectively. A single oscillation of q,,(r) was observed for the case R = 2, a feature which is imperceptible in Figure 2 but whose existence in such problems was remarked upon by Harris et al. [4]. Nomenclature
f F h
nondimensional, reduced streamfunction transformed function step length for 0 < T < 2~r2 aT~ step length in the z/direction for ~2LXT2 < T < %0 nondimensional t i m e increment for ~2ZxT2 < T < %0 p e r m e a b i l i t y of the porous m e d i u m number of grid spacings in the ~'direction for a_Tt_T < ~ < Too 2AT2 number of grid spacings in the ~direction for ~2 A T 2 < ~ <~oo number of grid spacings for 0 < ~ < ~T~  2AT2 variable coefficients in the governing equation for 2Z~T2 ~ < r < Too nondimensional heat flux at the plate ratio of the final surface t e m p e r a t u r e to the initial surface t e m p e r a t u r e local Rayleigh number sum of numerical solutions for t e m p e r a t u r e over consecutive time steps t e m p e r a t u r e of the fluid initial surface t e m p e r a t u r e (~ < 0) final surface t e m p e r a t u r e (r > O) characteristic t e m p e r a t u r e velocity components along x  and yaxes, respectively characteristic velocity Cartesian coordinates along the plate and normal to it, respectively 
K m ?%
N p, q q~ R
Ra~ S T
T1 T2 AT ILiad
uc x,y

Greek Symbols
6 £I ~ ~2 ~ £3
(,~ ~,O,gZ 0 v
O G
T
AT
effective t h e r m a l diffusivity of the fluidsaturated porous medium coefficient of thermal expansion boundarylayer thickness tolerances in the numerical schemes similarity variables expressions defined in Equation (23) nondimensional t e m p e r a t u r e function kinematic viscosity nondimensional t e m p e r a t u r e f u n c t i o n ratio of composite material heat capacity to convective fluid heat capacity time nondimensional time exact value of "r where transfer to the stepbystep method in 7/takes place nondimensional time increment for 0 < ~ <  2~T~ ~T2 relaxation p a r a m e t e r
Vol. 24, No. 4
¢
CONVECTION FROM A PLATE 1N POROUS MEDIA
551
strearnfunction
Subscripts 0 value at T = 0
i,j
e v a l u a t e d at t h e i t h and j t h nodal points in the ~/ and Tdirections, respectively e v a l u a t e d at t h e wall ambient condition
W O0
Superscripts *
p o i n t w h e r e t h e stepbystep n u m e r i c a l solution breaks down associated w i t h t h e stepbystep n u m e r i c a l solution in the ~1 and T variables associated w i t h t h e stepbystep n u m e r i c a l solution in the ~ and T variables 1.0
0.8
~ b
",...,.,
.......
/~x,~ / "~ / / ~_..~
~ + .... _ _
t
0.6
Steady state solution at ~ = 0, f~(q) Numerical solution at r = 0.006150 Numerical solution eA 7 = 0.073798 Numerical solution at T = 0.51!{)46 Steady _ state solution at large T Rf;(~?n½)] _ o
, 7 JU\~
o~ 0 . 4 
0.20.0
'
0.0
I
'
0.5
!
'
1.0
2.0~ ~ , ,
[
1.5 / \ ~ ,
[ / [
. . ~ ~ , , , ,
I
1.5
....... o [] ~ ....
'
I
2.0
'
r
'
2.5
I
'
3.0
I
3.5
'
L
4.0
Steady state solution at ~ = 0, ]~(77) Numerical solution at T = 0.003718 Numerical solution at T = 0.045126 Numerical solution at ~ = 0.203846 Steady state solution at large T, RJ~(~R½)
~' 1.0
0.5
0.0
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
4.0
77 FIG. 1 V a r i a t i o n of t h e n o n  d i m e n s i o n a l t e m p e r a t u r e 8(7 , ~) as a function of 77 at various values of ~" and t h e s t e a d y s t a t e solutions at ~ = 0 and T = oo. qH
(a) a =
= 0.5,
(b) R =
= 2.
552
S.D. Harris, D.B. Ingham and I. Pop
Vol. 24, No. 4
4.0[
3.0~ t.
Numerical solutions ] Small time solution, i.e. Equation (19) Steady state solution at large 7, R~f~'(0)
•. .
zo., "T ~ •
J
[....................
I
0.0
0.2
3=
0.4
0.6
0.8
1.0
I
i
1.2
1.4
1.6
1.8
3
i
~"
T = T"
i
4.0~ / •"'
?(
8.O I~ 0.0
Numerical solutions Small time solution, i.e. Equation [19) Steady
, 0.1
0.2
state
solution
I 0.3
at large
v,
z_ tr R~f~(O)
0.4
0.5
FIG. 2 Variation of the nondimensional heat flux from the plate q,~(r) as a function of r, the small time solution and the steady state solutions at ~ = 0 and ~ = oo, where the transition between solution methods occurs at the indicated times. qH
(a) R = ~  o.~,
(b) ~ = g  2
Re{erences 1. D.A. Nield a n d A. Bejan, Convection in Porous Media. Springer, New York (1992). 2. A. N a k a y a m a , PCAided Numerical Heat Transfer and Convective Flow. CRC Press, Tokyo (1995). 3. D.B. I n g h a m , J.H. M e r k i n a n d I. Pop, Int. J. Heat Mass Transfer 25, 1916 (1982). 4. S.D. Harris, D.B. I n g h a m a n d I. Pop, Transport in Porous Media (in press) (1996). 5. T.Y. Na a n d I. Pop, Int. Comm. Heat Mass Transfer 23, 697 (1996). 6. J.H. Merkin, Int. J. Heat Mass Transfer 1 5 , 9 8 9 (1972). 7. S.C.R. Dennis, J. Inst. Math. Applics. 10, 105 (1972).
Received October 14, 1996