Jtam-A4.dvi


JOURNAL OF THEORETICAL

AND APPLIED MECHANICS

51, 1, pp. 53-61, Warsaw 2013

ON THE EFFECTIVENESS OF VARIATION OF PHYSICAL VARIABLES ON

STEADY FLOW BETWEEN PARALLEL PLATES WITH HEAT TRANSFER

IN A POROUS MEDIUM

Hazem Ali Attia

El-Fayoum University, Faculty of Engineering, El-Fayoum, Egypt

e-mail: ah1113@yahoo.com

Mostafa A.M. Abdeen

Cairo University, Faculty of Engineering, Giza, Egypt

e-mail: mostafa a m abdeen@hotmail.com

The effect of variation of physical variables on a steady flow through a porousmediumwith
heat transfer betweenparallel plates is examined.Theviscosityand the thermal conductivity
are assumed to be temperature dependent. A constant pressure gradient is applied in the
axial direction and the twoplates are kept at two constant but different temperatures, while
the viscous dissipation is considered in the energy equation. A numerical solution for the
governingnon-linearcoupledequationsofmotionandtheenergyequation isdetermined.The
effect of porosity of themedium, the variable viscosity, and the variable thermal conductivity
on both the velocity and temperature distributions is reported.

Key words: variable properties, porous medium, heat transfer, parallel plates, steady state

List of symbols

a,b – viscosity and thermal conductivity parameter
cp – specific heat at constant pressure
Ec,Pr – Eckert and Prandtl number
k – thermal conductivity
K – Darcy permeability
P – pressure gradient
M – porosity parameter
T – temperature of fluid
T1,T2 – temperature of lower and upper plate
u – velocity component in x-direction
v – velocity of fluid
x,y – axial and distance in vertical direction
µ,ρ – viscosity and density of fluid

1. Introduction

The flow between infinite horizontal parallel plates has received great attention due to its ap-
plications in magnetohydrodynamic (MHD) power generators and pumps, etc. Hartmann and
Lazarus (1937) examined the effect of a transverse uniform magnetic field on the flow of a vi-
scous incompressible electrically conducting fluid between two infinite parallel plates. Analytical
solutions for the velocity fields were developed by Tao (1960), Alpher, 1961, Sutton and Sher-
man (1965), Cramer and Pai (1973) under different physical effects. Some exact and numerical



54 H.A. Attia, M.A.M. Abdeen

solutions for the heat transfer problemwere obtained inNigam and Singh (1960). Soundalgekar
et al. (1970), Soundalgekar and Uplekar (1986) studied the effect of Hall currents on the steady
MHDCouette flowwith heat transfer. The temperatures of the two plates were assumed either
to be constant (Soundalgekar et al., 1970) or linearly varying along the plates in the direction
of the flow (Soundalgekar and Uplekar, 1986). Attia (1998) studied the Hall current effects on
the velocity and temperature fields of an unsteady Hartmann flow with uniform suction and
injection.
Most of these studies were based on constant physical properties, however more accurate

prediction for the flow and heat transfer can be obtained considering the dependence of the
physical properties on temperature (Herwig and Wicken, 1986). Klemp et al. (1990) examined
the influence of temperature dependent viscosity on the entrance flow in a channel in the hydro-
dynamic case. Attia and Kotb (1996) handled the steady MHD fully developed flow and heat
transfer between two parallel plates with temperature dependent viscosity which was extended
to the transient state byAttia (1999). The effect of variation in physical properties in theMHD
Couette flow between parallel plates was studied by Joaqúın et al. (2010), Eguia et al. (2011).
In this paper, the steady flow of a viscous incompressible fluid with heat transfer through a

porousmedium is examined.Theflow in the porousmediumdealswith the analysis inwhich the
governing differential equations are based on Darcy’s law which accounts for the drag exerted
by the porous medium (Joseph et al., 1982; Ingham and Pop, 2002; Khaled and Vafai, 2003).
The viscosity and the thermal conductivity are assumed to change with temperature. The two
plates are kept at two constant but different temperatures, and a constant pressure gradient is
applied in the axial direction. The coupled nonlinear system of equations of motion and energy
equation including the viscous dissipation is solved numerically using finite differences to obtain
the velocity and temperature distributions.

2. Formulation of the problem

The fluid is flowing between two infinite parallel plates located at the y =±h planes through
a porousmediumwhere the Darcy model is assumed (Khaled and Vafai, 2003). The two plates
are kept at two constant temperatures T1 for the lower plate and T2 for the upper plate with
T2 >T1, and a constant pressure gradient is applied in the x-direction. The viscosity of the fluid
is assumed to vary with temperature exponentially, while the thermal conductivity is assumed
to depend on temperature linearly.
The flow of the fluid is governed by the Navier-Stokes equation

ρ(v ·∇v)=−∇P +∇· (µ∇v)+ fb (2.1)

where v is the velocity vector, P is the pressure, µ is the viscosity of the fluid, and fb is the
body force per unit volume. Using Eq. (2.1) the component of theNavier-Stokes equation in the
x-direction is given by

−
dP

dx
+
d

dy

(
µ
du

dy

)
−
µ

K
u=0 (2.2)

where K is theDarcy permeability (Khaled andVafai, 2003). Theno-slip condition at the plates
implies that

y=−h :u=0 y=h :u=0 (2.3)

The energy equation for the fluid considering the viscous dissipation is given by Ingham and
Pop (2002)

1

ρcp

[ d
dy

(
k
dT

dy

)
+µ
(du
dy

)2]
=0 (2.4)



On the effectiveness of variation of physical variables... 55

where T is the temperature of the fluid, cp is the specific heat at constant pressure of the
fluid, and k is the thermal conductivity of the fluid. The last term in the left side of Eq. (2.4)
represents the viscous dissipation respectively.

The temperature of the fluidmust satisfy the boundary conditions

T =T1 y=−h and T =T2 y=h (2.5)

The viscosity of the fluid is assumed to vary with temperature in the form µ = µof1(T),
and µo is the viscosity of the fluid at T = T1. By assuming the viscosity to vary exponentially
with temperature, the function f1(T) takes the form f1(T) = exp(−a1(T −T1)), see Klemp
et al. (1990), White (1991). The thermal conductivity of the fluid varies with temperature as
k = kof2(T), where ko is the thermal conductivity of the fluid at T = T1. We assume a linear
dependence for the thermal conductivity with temperature in the form k= ko(1+ b1(T −T1))
(Klemp et al., 1990;White, 1991), where the parameter b1may be positive or negative (Klemp
et al., 1990; White, 1991).

Introducing the following non-dimensional quantities

ŷ=
y

h
P̂ =
Pρh2

µ2o
û=
uρh

µo
T̂ =

T −T1
T2−T1

α=−
dP̂

dx̂

and
f̂1(T)= exp(−a1(T2−T1)T)= exp(−aT), a is the viscosity parameter,

f̂2(T)= 1+ b1(T2−T1)T =1+ bT , b is the thermal conductivity parameter,

M =µoh
2/K is the porosity parameter,

Pr=µocp/ko is the Prandtl number,

Ec=µ2o/[h
2cpρ

2(T2−T1)] is the Eckert number,

τL =(∂û/∂ŷ)ŷ=−1 is the axial skin friction coefficient at the lower plate,

τU =(∂û/∂ŷ)ŷ=1 is the axial skin friction coefficient at the upper plate,

NuL =(∂T/∂ŷ)ŷ=−1 is the Nusselt number at the lower plate,

NuU =(∂T/∂ŷ)ŷ=1 is the Nusselt number at the upper plate.

Equations (2.3) to (2.5) read (the hats are dropped for simplicity)

α+f1(T)
d2u

dy2
+
df1(T)

dy

du

dy
−f1(T)Mu=0 (2.6)

y=−1 u=0 and y=1 u=0 (2.7)

1

Pr
f2(T)

d2T

dy2
+
1

Pr

df2(T)

dy

dT

dy
+Ecf1(T)

(du
dy

)2
=0 (2.8)

T =0 y=−1 and T =1 y=1 (2.9)

Equations (2.6) and (2.8) represent a coupled system of non-linear ordinary differential equ-
ations which are solved numerically under boundary conditions (2.7) and (2.9) using finite dif-
ferences. A linearization technique is first applied to replace the nonlinear terms at a linear
stage, with corrections incorporated in subsequent iterative steps until convergence is reached.
Then the Crank-Nicolson implicit method is applied and an iterative scheme is used to solve
the linearized system of difference equations. The resulting block tri-diagonal system is solved
using the generalized Thomas-algorithm (Ames, 1977). Finite difference equations relating the
variables are obtained by writing the equations at the mid point of the computational cell and
then replacing the different terms by their second order central difference approximations in the
y-direction. The computational domain is divided intomeshes each of dimension ∆y.We define
the variables v= du/dy and H = dT/dy to reduce the second order differential Eqs. (2.6) and



56 H.A. Attia, M.A.M. Abdeen

(2.8) to first order differential equations. The finite difference representations for the resulting
first order differential take the form

α+
f1(T)i+1+f1(T)i

2

vi+1−vi
∆y

+
f1(T)i+1−f1(T)i

∆y

ui+1−ui
2

−M
f1(T)i+1+f1(T)i

2

ui+1−ui
2

=0

1

Pr

f2(T)i+1+f2(T)i
2

Hi+1−Hi
∆y

+
1

Pr

f2(T)i+1−f2(T)i
∆y

Hi+1+Hi
2

+Ec
f1(T)i+1+f1(T)i

2

vi+1+vi
2

vi+1+vi
2

=0

(2.10)

The variables with bars are given initial guesses. The linearized system of difference equations
represents a banded matrix that can be solved by different methods. One of the very powerful
methods to solve this system is the Thomas algorithm which has the advantage of reducing
memory storage and computational time. Equations (2.10) can bewritten in the following form

a1ui+a2ui+1+a3vi+a4vi+1+a5Ti+a6Ti+1 = a7

b1Ti+ b2Ti+1+ b3Hi+ b4Hi+1+ b5ui+ b6ui+1 = b7
(2.11)

where the a’s and b’s are the coefficients of the difference equations. We set the following diffe-
rence equations for evaluating the unknowns u, v, T ,H

ui =uivi+ ũiHi+ ûi Ti =T ivi+ T̃iHi+ T̂i

vi = vi+1vi+1+ ṽi+1Hi+1+ v̂i+1 Hi =Hi+1vi+1+ H̃i+1Hi+1+ Ĥi+1
(2.12)

where ui, ũi, ûi, Ti, T̃i, T̂i, vi, ṽi, v̂i, Hi, H̃i, Ĥi are the Thomas coefficients. The equations
relating the velocity and temperature to their derivatives using Eq. (2.12) are given as

ui+1 =
(
ui+
∆

2

)
vi+ ũiHi+

∆

2
vi+1+ ûi

Ti+1 =T ivi+
(
T̃i+
∆

2

)
Hi+

∆

2
Hi+1+ T̂i

(2.13)

Substituting Eqs. (2.12) and (2.13) into Eqs. (2.11) and solving the resulting equations for the
unknown coefficients in Eqs. (2.12)3,4, we obtain

vi+1 =
ℓ31h2−h31l2
l1h2− l2h1

ṽi+1 =
ℓ32h2−h32l2
l1h2− l2h1

v̂i+1 =
ℓ33h2−h33l2
l1h2− l2h1

Hi+1 =
h31ℓ1−h31h1
l1h2− l2h1

H̃i+1 =
h32ℓ1− ℓ32h1
l1h2− l2h1

Ĥi+1 =
h33ℓ1− ℓ33h1
l1h2− l2h1

where

l1 = a1ui+a2
(
ui+
∆

2

)
+a3+(a5+a6)T i l2 =(a1+a2)ũi+a5T̃i+a6

(
T̃i+
∆

2

)
+a7

l31 =−
(
a2
∆

2
+a4
)

l32 =−
(
a6
∆

2
+a8
)

l33 = a9− (a1+a2)ûi− (a5+a6)T̂i h1 = b5ui+ b6
(
ui+
∆

2

)
+ b7+(b1+ b2)T i

h2 =(b5+ b6)ũi+b1T̃i+ b2
(
T̃i+
∆

2

)
+ b3 h31 =−

(
b6
∆

2
+ b8
)

h32 =−
(
b2
∆

2
+ b4
)

h33 = b9− (b5+ b6)ûi− (b1+ b2)T̂i



On the effectiveness of variation of physical variables... 57

ThenEqs. (2.13) can be used to determine the rest of the unknown coefficients in Eqs. (2.12)1,2.
Computation of the coefficients starts from the lower plate (y = −1) with known boundary
conditions for both the velocity and temperature problems and proceeds up to the upper plate
(y=1). Consequently, the computation of the unknown functions is done from the upper plate
up till the lower one.

In this study, computations have been carried out for α = 5, Pr = 1 and Ec = 0.2. Grid-
independence studies show that the computational domain −1 < y < 1 can be divided into
intervals with the step size ∆y=0.005. Smaller step sizes do not show any significant change in
the results. Convergence of the scheme is assumedwhen all of the unknowns u, v, T and H for
the last two approximations differ fromunityby less than 10−6 for all values of y in −1<y< 1.
Less than 7 approximations are required to satisfy these convergence criteria for all ranges of
the parameters studied here.

3. Results and discussion

Figure 1 shows the velocity profile for various values of the parameters a and M and for b=0.
It is clear that an increase in the parameter M decreases the velocity u for all values of a as a
result of the increasing damping force on u. On the other hand, an increase in the parameter a
increases u for all M due to the decrease in the viscosity. The influence of the parameter a on u
is more apparent for smaller porosity parameters. The parameter a has a pronounced effect on
the symmetry of the velocity profile about the plane y=0.

Fig. 1. Variation of u with y for various values of a; (a) b=0,M =0, (b) b=0,M =1,
(c) b=0,M =2

Figure 2 shows the temperature profile for various values of the parameters a and M and
for b = 0. An increase in the parameter M decreases T as a result of decreasing u and, in
turn, decreases the viscous dissipation. On the other hand, increasing the parameter a cause an
increase in T for all values of M due to its effect in increasing u and consequently the viscous
dissipation. It is also clear also that the effect of a on T is more clear for smaller M.

Figure 3 presents the temperature profile for different values of the parameters b and M
and for a=0. The figure depicts that the effect of b on T depends on t, and that increasing b
increases T for small times, but decreases T with time progression. This occurs because, for
small time, the centre of the channel acquires heat by conduction from the upper hot plate, but
after large time, when u is large, the viscous dissipation is large at the centre and, consequently
centre looses heat by conduction. It is also clear that the effect of b on T is more pronounced
for smaller M. The parameter b has no significant effect on u in spite of the coupling between
the momentum and energy equations.



58 H.A. Attia, M.A.M. Abdeen

Fig. 2. Variation of T with y for various values of a; (a) b=0,M =0, (b) b=0,M =1,
(c) b=0,M =2

Fig. 3. Variation of T with y for various values of b; (a) a=0,M =0, (b) a=0,M =1,
(c) a=0,M =2

Table 1presents thevariation of temperatureat the centre of the channelwith theparameters
a and b for M = 1. In Table 1, T increases with increasing a for all values of b. On the
other hand, for smaller values of a, increasing b increases T , while for higher values of a,
increasing b decreases it. This is because a decrease in a decreases the velocity and its gradient
which decreases dissipation and makes the centre gain heat by conduction. Higher values of a
increases dissipation, and the centre looses heat by conduction which results in a decrease in T
with the increment of b.

Table 1. Variation of the steady state temperature T at y = 0 for various values of a and b
(M =1)

T a=−0.5 a=−0.1 a=0.0 a=0.1 a=0.5

b=−0.5 0.6315 0.6885 0.7086 0.7322 0.8822

b=−0.1 0.6475 0.6902 0.7046 0.7210 0.8171

b=0.0 0.6514 0.6916 0.7049 0.7202 0.8082

b=0.1 0.6549 0.6927 0.7052 0.7195 0.8007

b=0.5 0.6638 0.6945 0.7045 0.7157 0.7777

Tables 2 and 3 show the influence of the parameters a and b on the skin friction coefficients
at both walls τL and τU, respectively, for M = 1. Increasing a increases τL and the magni-
tude of τU for all b. The effect of b on τL and τU depends on a. For small a, increasing b
increases τL but decreases the magnitude of τU. On the other hand, for higher values of a,



On the effectiveness of variation of physical variables... 59

increasing b decreases τL but increases the magnitude of τU. Tables 4 and 5 show the effect of
the parameters a and b on the Nusselt numbers at both walls NuL and NuU, respectively, for
M =1. Increasing a increases NuL and themagnitude of NuU. However, although an increase
in b decreases the magnitude of NuU, its effect on NuL depends on the value of b. Therefore,
increasing b decreases NuL, but increasing b increases NuL more.

Table 2. Variation of the steady state skin friction coefficient at the lower plate τL for various
values of a and b (M =1)

τL a=−0.5 a=−0.1 a=0.0 a=0.1 a=0.5

b=−0.5 3.4975 3.7445 3.8098 3.8764 4.1570

b=−0.1 3.5081 3.7483 3.8097 3.8717 4.1226

b=0.0 3.5096 3.7488 3.8098 3.8710 4.1184

b=0.1 3.5111 3.7493 3.8097 3.8705 4.1143

b=0.5 3.5164 3.7509 3.8097 3.8688 4.1031

Table 3.Variation of the steady state skin friction coefficient at the lower plate τU for various
values of a and b (M =1)

τU a=−0.5 a=−0.1 a=0.0 a=0.1 a=0.5

b=−0.5 −2.4874 −3.4992 −3.8097 −4.1478 −5.8263

b=−0.1 −2.4739 −3.4942 −3.8097 −4.1545 −5.8958

b=0.0 −2.4718 −3.4936 −3.8097 −4.1553 −5.9022

b=0.1 −2.4699 −3.4929 −3.8097 −4.1559 −5.9076

b=0.5 −2.4645 −3.4913 −3.8097 −4.1579 −5.9197

Table 4.Variation of the steady stateNusselt number at the lower plate NuL for various values
of a and b (M =1)

NuL a=−0.5 a=−0.1 a=0.0 a=0.1 a=0.5

b=−0.5 1.1654 1.3366 1.3910 1.4521 1.7902

b=−0.1 1.1334 1.3056 1.3593 1.4185 1.7358

b=0.0 1.1347 1.3073 1.3609 1.4199 1.7342

b=0.1 1.1380 1.3110 1.3646 1.4234 1.7350

b=0.5 1.1611 1.3354 1.3888 1.4474 1.7511

Table 5.Variation of the steady stateNusselt number at the upperplate NuU for various values
of a and b (M =1)

NuU a=−0.5 a=−0.1 a=0.0 a=0.1 a=0.5

b=−0.5 −0.3986 −0.7259 −0.8288 −0.9434 −1.5755

b=−0.1 −0.1596 −0.3464 −0.4046 −0.4689 −0.8189

b=0.0 −0.1396 −0.3084 −0.3609 −0.4189 −0.7321

b=0.1 −0.1255 −0.2796 −0.3274 −0.3802 −0.6632

b=0.5 −0.1006 −0.2148 −0.2499 −0.2886 −0.4919

4. Conlusions

The steadyflow throughaporousmediumbetween twoparallel plates is examined.Theviscosity
and thermal conductivity of the fluid are assumed to vary with temperature. The effects of the
porosity parameter M, the viscosity parameter a, and the thermal conductivity parameter b



60 H.A. Attia, M.A.M. Abdeen

on the velocity and temperature fields are investigated. It is found that the parameters a and b
have a more apparent effect on the velocity and temperature distributions for smaller values
of M. On the other hand, the parameter b has no significant effect on the velocity u, however,
it has a marked effect on the temperature, and its effect greatly depends on the parameters M
and a. It is of interest to find that the variation of the Nusselt number at the lower plate with
the thermal conductivity parameter b depends on the values of b.

References

1. Alpher R.A., 1961, Heat transfer in magnetohydrodynamic flow between parallel plates, Inter-
national Journal of Heat and Mass Transfer, 3, 108

2. AmesW.F., 1977,Numerical Solutions of Partial Differential Equations, 2ndEd.,AcademicPress,
NewYork

3. Attia H.A., 1998, Hall current effects on the velocity and temperature fields of an unsteady
Hartmann flow,Can. J. Phys., 76, 739-746

4. Attia H.A., 1999, Transient MHD flow and heat transfer between two parallel plates with tem-
perature dependent viscosity,Mechanics Research Communications, 26, 1, 115-121

5. Attia H.A., Kotb N.A., 1996, MHD flow between two parallel plates with heat transfer, Acta
Mechanica, 117, 215-220

6. CramerK., Pai S., 1973,MagnetofluidDynamics for Engineers andApplied Physicists,McGraw-
Hill

7. Eguia P., Zueco J., Granada E., Patiõ D., 2011, NSM solution for unsteady MHD Couette
flow of a dusty conducting fluid with variable viscosity and electric conductivity, Applied Mathe-
matical Modelling, 35, 303-316

8. Hartmann J., Lazarus F., 1937, Kgl. Danske Videnskab. Selskab,Mat.-Fys. Medd., 15, 6/7

9. Herwig H., Wicken G., 1986, The effect of variable properties on laminar boundary layer flow,
Warme-und Stoffubertragung, 20, 47-57

10. Ingham D.B., Pop I., 2002,Transport Phenomena in Porous Media, Pergamon, Oxford

11. Joaqúın Z., Pablo E., Enrique G., José L.M., Osman A.B., 2010, An electrical network for
the numerical solution of transientmhd couette flow of a dusty fluid: Effects of variable properties
and hall current, Int. Comm. Heat and Mass Trans., 37, 1432-1439

12. Joseph D.D., Nield D.A., Papanicolaou G., 1982, Nonlinear equation governing flow in a
saturated porousmedia,Water Resources Research, 18, 4, 1049-1052

13. Khaled A.R.A., Vafai K., 2003, The role of porousmedia inmodeling flow and heat transfer in
biological tissues, Int. J. Heat Mass Transf., 46, 4989-5003

14. KlempK.,HerwigH., SelmannM., 1990,Entranceflow inchannelwith temperaturedependent
viscosity including viscous dissipation effects, Proceedings of the Third International Congress of
Fluid Mechanics, Cairo, Egypt, 3, 1257-1266

15. Nigam S.D., Singh S.N., Heat transfer by laminar flow between parallel plates under the action
of transversemagnetic field,Quart. J. Mech. Appl. Math., 13, 85

16. SoundalgekarV.M.,UplekarA.G., 1986,Hall effects inMHDCouetteflowwithheat transfer,
IEEE Transactions on Plasma Science,PS-14, 5

17. SoundalgekarV.M., Vighnesam N.V., TakharH.S., 1979, Hall and ion-slip effects inMHD
Couette flowwith heat transfer, IEEE Transactions on Plasma Sciences,PS-7, 3

18. Sutton G.W., Sherman A., 1965,Engineering Magnetohydrodynamics, McGraw-Hill

19. Tao L.N., 1960, Magnetohydrodynamic effects on the formation of Couette flow, Journal of Ae-
rospace Sci., 27, 334

20. White M.F., 1991,Viscous Fluid Flow, McGraw-Hill



On the effectiveness of variation of physical variables... 61

O efektywności uzmienniania parametrów fizycznych stacjonarnego przepływu przez

ośrodek porowaty z uwzględnieniem przesyłu ciepła pomiędzy równoległymi płytami

Streszczenie

W pracy zbadano wpływ zmiany parametrów fizycznych na ustalony przepływ przez porowaty ośro-
dek przy jednoczesnym przesyle ciepła pomiędzy równoległymi płytami ograniczającym obszar ośrodka.
Założono, że lepkość i współczynnik przewodnictwa cieplnego są zależne od temperatury. Do utrzyma-
nia stacjonarnego przepływu wprowadzono stały gradient ciśnienia w kierunku osiowym, a dwie płyty
utrzymano w stałych, lecz różniących się temperaturach. Efekt dyssypacji wiskotycznej uwzględniono
w równaniu energii. Zaprezentowano numeryczne rozwiązanie nieliniowych, sprzężonych równań ruchu
oraz równania energii. Omówiono wpływ porowatości ośrodka, zmiennej lepkości i współczynnika prze-
wodnictwa cieplnego na rozkład prędkości i temperatury w obszarze przepływu.

Manuscript received December 3, 2011; accepted for print March 7, 2012