Natural convection in a shallow cavity filled with a micropolar fluid

International Journal of Heat and Mass Transfer 53 (2010) 2750–2759

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Natural convection in a shallow cavity filled with a micropolar fluid Z. Alloui *, P. Vasseur Ecole Polytechnique, C.P. 6079, Succ ‘‘Centre Ville”, Montréal, Que., Canada H3C 3A7

a r t i c l e

i n f o

Article history: Received 6 November 2009 Received in revised form 5 February 2010 Accepted 5 February 2010 Available online 17 March 2010 Keywords: Natural convection Micropolar fluid layer

a b s t r a c t This paper reports an analytical and numerical study of natural convection in a shallow rectangular cavity filled with a micropolar fluid. Neumann boundary conditions for temperature are applied to the horizontal walls of the enclosure, while the two vertical ones are assumed insulated. The governing parameters for the problem are the thermal Rayleigh number, Ra, the Prandtl number, Pr, the aspect ratio of the cavity, A and various material parameter of the fluid. For convection in an infinite layer (A  1), analytical solutions for the stream function temperature and angular velocity are obtained using a parallel flow approximation in the core region of the cavity and an integral form of the energy equation. The critical Rayleigh numbers for the onset of supercritical convection are predicted explicitly by the present model. Furthermore, a linear stability analysis is conducted yielding numerically the critical Rayleigh numbers for the onset of motion. Also, results are obtained from the analytical model for finite-amplitude convection for which the flow and heat transfer are presented in terms of the governing parameters of the problem. Numerical solutions of the full governing equations are obtained for a wide range of the governing parameters. A good agreement is observed between the analytical model and the numerical simulations. The influence of the material parameters on the flow and heat transfer is demonstrated to be significant. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Micropolar fluid theories were introduced by Eringen [1–3] in order to describe some physical systems which do not satisfy the Navier–Stokes model. In industry examples of applications of such fluids include the extrusion of polymer fluids, continuous casting glass-fiber, animal bloods, porous media, cooling of a metallic plate in a bath production, paper production, metal extraction, colloidal and suspension solutions, etc. A review of the various applications of micropolar fluid mechanics has been presented by Ariman et al. [4]. Eringen [5] later generalized the micropolar fluid theory to include thermal effects. The first application of this last model seems to be due to Balaram and Sastry [6] who investigated the problem of convective heat transfer of a micropolar fluid in a vertical channel. The horizontal configuration was considered by Maiti [7]. More recently, natural convection of a micropolar fluid, confined in a vertical rectangular enclosure heated differentially from the vertical sides, was considered numerically by Hsu and Chen [8]. Parametric studies of the effects of microstructure on the fluid flow and heat transfer in the enclosure have been conducted by these authors. It was found that an increase in the vortex viscosity parameter gives rise to a decrease in heat transfer rate. On the other hand it was demonstrated that the increase of spin-gradient viscosity increases the heat flux. In another work Hsu et al. [9] * Corresponding author. Tel.: +1 514 340 4711 3619; fax: +1 541 340 5917. E-mail address: [email protected] (Z. Alloui). 0017-9310/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijheatmasstransfer.2010.02.026

investigated the case of an inclined cavity with uniform heat sources applied on the system. Various boundary conditions, for microrotation, were considered and they were found to have a significant influence on the heat transfer. Natural convective flow and heat transfer of micropolar fluids confined in a square cavity, subject to various thermal conditions, was investigated numerically by Aydin and Pop [10,11] and Zadravec et al. [12]. In all these studies it was demonstrated that the microrotation of particles in suspension in general decreases the overall heat transfer. The Rayleigh Benard problem of onset of convective instabilities in a horizontal layer of a micropolar fluid heated from below has also be considered in the past. Rama Rao [13] studied the onset of convection of a heat conducting micropolar fluid layer confined between two horizontal rigid boundaries. The instability of rotating micropolar fluids has been investigated by Sastry and Ramamohan [14] and Qin and Kaloni [15]. The effect of throughflow and magnetic field on the onset of Benard convection in a horizontal layer of micropolar fluid confined between two rigid, isothermal and microrotation free, boundaries has been studied by Murty [16]. The critical Rayleigh number was predicted on the basis of a single-term Galerking technique. The effects of non-uniform temperature profiles on Marangoni convection in micropolar fluids confined between a lower rigid isothermal boundary and an upper free, constant heat flux boundary was investigated by Rudraiah et al. [17]. It was demonstrated by these authors that micropolar fluids heated from below are more stable when compare to the pure viscous fluid situation. Rayleigh Benard convection in a

2751

Z. Alloui, P. Vasseur / International Journal of Heat and Mass Transfer 53 (2010) 2750–2759

Nomenclature A B C g 0 H j k K 0 L N Nu n Pr q0 Ra RaC t T u v

aspect ratio of the cavity, L0 =H0 micro-inertia parameter, H02 =j constant temperature gradient in x direction gravitational acceleration height of fluid layer micro inertia per unit mass thermal conductivity vortex viscosity parameter, j/l width of fluid layer dimensionless angular velocity, N0 H02 =a Nusselt number, Eq. (27) dimensionless micro-gyration parameter, Eq. (5) Prandtl number, m/a constant heat flux per unit area thermal Rayleigh number, gb0T DT 0 H03 =am critical Rayleigh number, Eq. (26) dimensionless time, t 0 a=H02 dimensionless temperature, ðT  T 00 Þ=DT 0 dimensionless velocity in x direction, ðu0 H0 =aÞ dimensionless velocity in y direction, ðv 0 H0 =aÞ

micropolar ferromagnetic fluid has been investigated analytically by Abraham [18] for a layer with free–free, isothermal, spin-vanishing magnetic boundaries. The influence of the micropolar parameters on convection in the ferromagnetic case is similar to its role in the non-magnetic fluid case. More recently, the effect of a non-uniform basic temperature gradient on the onset of Marangoni convection in a horizontal micropolar fluid layer was considered by Melviana et al. [19]. It was assumed that the layer is bounded below by a rigid plate and above by a non-deformable free surface subjected to a constant heat flux. At these boundaries the microrotation is assumed to be vanished. The influence of various parameters on the onset of convection is discussed. Also, a linear stability analysis was performed by Idris et al. [20] to study the effect of non-uniform basic temperature gradients on the onset of Benard–Marangoni convection in a micropolar fluid. The influence of various parameters on the onset of convection has been analysed by these authors. It was found that the presence of micronsized suspended particles delays the onset of convection. In the present paper, we consider natural convection in a horizontal layer of a micropolar fluid with the horizontal boundaries heated and cooled by constant heat fluxes (Neumann boundary conditions). The paper is organized as follows. In the next sections, the formulation of the problem and numerical method are presented. An approximate analytical solution is then proposed. A linear stability analysis is conducted to predict numerically the critical Rayleigh numbers for the onset of motion from the rest state. The last section contains some concluding remarks.

dimensionless coordinate axis, ðx0 =H0 Þ dimensionless coordinate axis, ðy0 =H0 Þ

x y

Greek symbols a thermal diffusivity thermal expansion coefficient b0T l dynamic viscosity m kinematic viscosity of fluid, l/q k material parameter, c/(lj) q density of fluid W dimensionless stream function, W0 =a c spin-gradient viscosity j vortex viscosity Subscripts 0 reference state c refers to critical conditions Superscript refers to dimensional variable

0

The governing equations which describe the system behaviour are conservation of momentum, microrotation and energy which 0 are given below in terms of the stream function W as (see for instance [12])

@ðr2 W0 Þ 1 j @T 0 þ JðW0 ; r2 W0 Þ ¼ ðl þ jÞr2 ðr2 W0 Þ þ r2 N0  gb0T 0 @t 0 q0 q0 @x ð1Þ @N0 c 2 0 j þ JðW0 ; N0 Þ ¼ r N  ð2N 0 þ r2 W0 Þ q0 j @t0 q0 j

ð2Þ

@T 0 þ JðW0 ; T 0 Þ ¼ ar2 T 0 @t0

ð3Þ

where, as usual, in order to satisfy the continuity equation, the 0 stream function W is defined such that u0 ¼ @ W0 =@y0 , 0 0 0 v ¼ @ W =@x and Jðf ; gÞ ¼ @f =@y@g=@x  @f =@x@g=@y. The appropriate boundary conditions applied on the walls of the layer are

x0 ¼ L0 =2;

u0 ¼ v 0 ¼ 0;

y0 ¼ H0 =2;

u0 ¼ v 0 ¼ 0;

ð4Þ ð5Þ

where 0 6 n 6 1 is a constant.

y ′, v′

2. Mathematical formulation The configuration considered in this study is a horizontal shallow cavity, of thickness H0 and width L0 filled with a micropolar fluid (see Fig. 1). The origin of the coordinate system is located at the centre of the cavity with x0 and y0 being the horizontal and vertical coordinates, respectively. The cavity is insulated on the sides and heated from the bottom by a uniform heat flux q0 . The micropolar fluid is assumed to satisfy the Boussinesq approximation. The density variation with temperature and concentration is described by the state equation q ¼ q0 ½1  b0T ðT 0  T 00 Þ where q0 is the fluid mixture density at temperature T 0 ¼ T 00 and b0T is the thermal expansion coefficients, respectively.

@v 0 @T 0 ; ¼0 0 @x @x0 0 @u @T 0 q0 N0 ¼ n 0 ; ¼ 0 @y @y k

N0 ¼ n

q′

H′

Micropolar fluid

x′, u′

q′ L′ Fig. 1. Schematic diagram of the problem domain and coordinate system.

2752

Z. Alloui, P. Vasseur / International Journal of Heat and Mass Transfer 53 (2010) 2750–2759

The case n = 0 represents concentrate particle flows in which the microelements close to the wall are unable to rotate [21]. The case n = 0.5 represents weak concentration and corresponds to the vanishing of antisymmetric part of the stress tensor [22]. Finally, it was postulated by Peddieson [23] that the case n = 1 is applicable to the modeling of turbulent boundary layer flows. The governing equations are non-dimensionalized by scaling 0 length by H , stream function by the thermal diffusivity a, microro0 0 tation by a/H 2 and time by H 2/a. Also, we introduce the reduced 0 0 0 temperature T ¼ ðT 0  T 0O Þ=DT 0 where DT = q H /k. The dimensionless equations governing the present problem then read

@ r2 W @T þ JðW; r2 WÞ ¼ Pr ð1 þ KÞr2 ðr2 WÞ þ Pr K r2 N  Pr Ra @x @t ð6Þ @N þ JðW; NÞ ¼ Pr kr2 N  Pr BKð2N þ r2 WÞ @t @T þ JðW; TÞ ¼ r2 T @t

ð7Þ ð8Þ

see for instance [25], the present solution is rather independent of this parameter provided that this latter is of order one of greater. Furthermore, this point is confirmed by the analytical solution discussed in the following section which, in its range of validity, is independent of Pr. Typical numerical results are presented in Fig. 2a–c for RaT = 5000, Pr = 7, A = 6, B = 1, k ¼ 0:1, n = 0 and various values of K. On the graphs, streamlines, isotherms and microrotation patterns are presented from top to bottom. The effects of varying the intensity of the vortex viscosity parameter K from 0 to 5, on the strength of convection (Wmax) and the resulting heat (Nu) are observed to be considerable. This point will be discussed in detail in the following sections. However, Fig. 2a–c clearly illustrate the fact that for a shallow cavity (A  1) the flow and angular velocity in the core of the enclosure is essentially parallel while the temperature is linearly stratified along the x direction. The numerically determined profiles of horizontal velocity u, temperature T and microrotation velocity, N at the center of the convective cell, are compared in Figs. 3–6 with their analytical counterparts derived below. The agreement between the two solutions is seen to be excellent.

The corresponding boundary conditions are

@W @v @T ¼0 ¼ 0; N ¼ n ; @x @x @x @W @u @T y ¼ 1=2 W ¼ ¼ 1 ¼ 0; N ¼ n ; @y @y @y

x ¼ A=2;

W¼

ð9Þ ð10Þ

The heat transfer rates can be expressed in terms of the Nusselt number and can be obtained from the following expression:

Nu ¼

q0 1 ¼ kDT 0 =H0 DT

ð11Þ

where DT = T(0,  1/2)  T(0, 1/2) is the temperature difference, evaluated at x = 0. This follows from the fact that DT is independent of x, such that they are arbitrarily evaluated at the center of the cavity. From the above equations it is observed that the present problem is governed by seven parameters, namely the thermal Rayleigh number Ra ¼ gb0T DT 0 H03 =am, Prandtl number Pr = m/a, aspect ratio of 0 0 the layer A = L /H , materials parameter k ¼ c=ðljÞ; vortex viscosity 0 parameter K = j/l, micro-inertia parameter B = H 2/j and microgyration parameter n.

4. Analytical solution In the limit of a shallow cavity, the governing equations (6)–(8) can be considerably simplified under the parallel flow assumption.

Ψ

Τ

(a)

Ν

Ψ

3. Numerical solution The solution of the governing equations and boundary conditions, Eqs. (6)–(10), is obtained using a control volume approach and SIMPLER algorithm [24]. A finite difference procedure with variable grid size is considered for better consideration of boundary conditions. The power-law scheme is used to evaluate the flow, heat and mass fluxes across each of the control volume boundaries. A second order back-wards finite difference scheme is employed to discretize the temporal terms appearing in the governing equations. A line-by-line tridiagonal matrix algorithm with relaxation is used in conjunction with iterations to solve the non-linear discretized equations. We consider that convergence is reached when

Τ

(b)

Ν

Ψ

PP i

zþ1 z j jWi;j  Wi;j j 6 106 P P zþ1 i j jWi;j j

ð12Þ

where the superscripts z and (z + 1) indicate the value of the zth and (z + 1)th iterations, respectively, i and j indices denote grid location in the (x, y) plane. Numerical tests, using various mesh sizes, were done for the same conditions in order to determine the best compromise between accuracy of the results and computer time. Thus, most of the calculations presented in this paper were performed using a 60  180 grid. All the numerical results presented in this study are limited to water-based solutions, i.e. Pr = 7. As discussed by many authors,

Τ

(c)

Ν

Fig. 2. Contour lines of stream function (top), temperature (middle) and microrotation (bottom) predicted by the numerical solution of the full governing equations for RaT = 5000, Pr = 7, A = 6, B = 1, k ¼ 0:1, n = 0 and (a) K = 0, Wmax = 4.00, Nu = 2.51, (b) K = 1, Wmax = 2.70, Nu = 2.03, (c) K = 5, Wmax = 1.06, Nu = 1.25.

2753

Z. Alloui, P. Vasseur / International Journal of Heat and Mass Transfer 53 (2010) 2750–2759

0.5

0.5

(c)

5

8

(a)

Analytical

Analytical

Numerical

Numerical

0.5 K=0

y

0

y

0

K=0 0.5 8

5

-0.5 -0.01

-0.005

0

0.005

-0.5

0.01

-0.0005

0

0.5

(b)

0.0005

θ*

u* 4

(d) Analytical

3.5

Numerical

3

K = 0,

5

8

0

y

Nu

0.5

2.5

0.005

K 0 1 5

2

*

N(y = 0)

1.5 10

-2

10

0

10

0

0.002

0.004

0.006

Numerical

2

K

-0.5 -0.002

Analytical

1

0.008

103

104

105

Ra

*

N

Fig. 3. Distribution of (a) horizontal dimensionless velocity component, u*, (b) dimensionless microrotation, N*, and (c) dimensionless temperature, h*, profiles at the vertical mid-plan (x = 0) of the cavity and (d) Nu versus Ra for B = 1, k ¼ 0:1, n = 0 and K = 0, 1, 5.

The procedure to obtain the analytical solution presented in this section has already been described in details (see for instance [26]) Thus, in the limit of, A  1, it is assumed that W(x, y)  W(y), N(x, y)  N(y) and T(x, y)  Cx + h(y), where C is unknown constant temperature gradient in x-direction. Using the above approximations, it is found that Eqs. (6)–(8) reduce to the following systems of ordinary differential equations: 4

ð1 þ KÞ 2

k

d N dy

2

2

dy

dy

4

þK

"

¼

ð18Þ ð19Þ

where

W ¼ XC 1 ½coshðXyÞ  coshðX0 Þ þ C 2 ð16y4  1Þ=192 þ C 3 ð4y2  1Þ=8 N  ¼ ð1 þ KÞC 1 X3 coshðXyÞ=K þ C 2 ½y2 þ k=ðKBÞ=2 þ C 3 =2 h ¼ C 1 sinhðXyÞ þ C 2 y5 =60 þ C 3 y3 =6  C 4 y

2

d N

¼ RaT C

2

dy

2

¼ BK 2N þ

2

d h

d W

hðyÞ ¼ Ra C 2 h  y uðyÞ ¼ Ra Cu

d W 2

dy

ð13Þ

#

u ¼ C 1 X2 sinhðXyÞ þ C 2 y3 =3 þ C 3 y

X ¼ fBKð2 þ KÞ=½kð1 þ KÞg1=2 ; X0 ¼ X=2 ð14Þ

ð20Þ and

dW C dy

ð15Þ

C 1 ¼ 2½BKð1  2nÞ þ 6keX0 =½12Bð2 þ KÞX2 ða1 eX þ a2 Þ; C 2 ¼ 1=ð2 þ KÞ

The solution of Eqs. (13)–(15), satisfying the boundary conditions (10), is

C 3 ¼ 2X2 C 1 sinhðX0 Þ  C 2 =12; C 4 ¼ XC 1 coshðX0 Þ þ ðC 2 =24 þ C 3 Þ=8

WðyÞ ¼ Ra C W

a1 ¼ nKðX þ 2Þ  XðK þ 1Þ  K; a2 ¼ nKðX  2Þ  XðK þ 1Þ þ K

NðyÞ ¼ Ra CN

ð16Þ 

ð17Þ

ð21Þ

2754

Z. Alloui, P. Vasseur / International Journal of Heat and Mass Transfer 53 (2010) 2750–2759

0.5

(a)

y

0.5

(c) Analytical

Analytical

Numerical

Numerical

0

0

y

n=0

0.5

1

1

0.5 n=0

-0.5

-0.004

-0.002

0

0.002

-0.5 -0.0004

0.004

-0.0002

0.5

(b)

0

0.0002

0.0004

θ*

u* 4

(d)

3.5

Analytical Numerical

3

y

0

Nu

2.5

n 0 0.5 1

2

1 0.5 -0.5 -0.05

-0.04

-0.03

-0.02

Ν

-0.01

Analytical Numerical

1.5

n=0 1 103

0

104

*

105

106

Ra

Fig. 4. Distribution of (a) horizontal dimensionless velocity component, u*, (b) dimensionless microrotation, N*, and (c) dimensionless temperature, h*, profiles at the vertical mid-plan (x = 0) of the cavity and (d) Nu versus Ra for K = 5, B = 1, k ¼ 0:1 and n = 0, 0.5, 1.

The parallel flow approximation is only applicable in the core of the layers. Flow in the end regions is much more complicated and cannot be approximated in such a simple manner. For this reason, the thermal boundary condition in the x-direction, Eq. (9), cannot be reproduced exactly with this approximation. We can however, impose an equivalent energy flux condition [27] in that direction

C¼

Z

1=2

uðyÞhðyÞdy

ð22Þ

1=2

1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1  A1 RaÞ=A2 Ra

A2 ¼ sinhðX0 Þ coshðX0 Þe1 þ coshðX0 Þe2 þ sinhðX0 Þe3 þ e4

e2 ¼ ½ð3840 þ 160X2 þ X4 ÞC 2

þ 40ð48 þ X2 ÞX2 C 3  960C 4 X4 C 1 =ð240X3 Þ e3 ¼ ½ð768 þ 96X2 þ X4 ÞC 2  24ð16 þ X2 ÞX2 C 3 þ 192C 4 X4 C 1 =ð24X4 Þ

ð25Þ

e4 ¼ C 21 =ð2X3 Þ þ C 22 =414720 þ C 23 =480  C 3 C 4 =12 Eq. (23) leads us to the critical value for existence of convection in the layer. The supercritical Rayleigh number Rac, for the onset of motion from the rest state, can be obtained from the following expression:

ð23Þ Rac ¼

where

A1 ¼ ½sinhðX0 Þ  X coshðX0 ÞC 1  C 2 =240  C 3 =12

e1 ¼ XC 21 ;

þ 13C 2 C 3 =80640  C 2 C 4 =240

at any x. Substituting Eqs. (18, (19) into Eq. (22) and integrating yields, after some straightforward but laborious algebra, an expression of the form

C¼

and

ð24Þ

1 A1

ð26Þ

corresponding to C = 0 in Eq. (23). The local heat transfer rate is obtained according, to Eqs. (11) and (18), by

2755

Z. Alloui, P. Vasseur / International Journal of Heat and Mass Transfer 53 (2010) 2750–2759

(a)

0.5

(c)

0.5

0.1

0.01

Analytical

Analytical

Numerical

Numerical

B=1

y

y

0

0

B=1 0.1 0.01 -0.5

-0.5 -0.001

0

-0.0001

0.001

(b)

0

0.0001

θ*

u* 0.5

(d)

3.5

Analytical Numerical

3

2.5

y

0

0.01

B=1

0.1

Nu 2

B 0.01 0.1 1

Analytical Numerical

1.5

-0.5 -0.002

0

*

0.002

0.004

N

1

104

Ra

105

Fig. 5. Distribution of (a) horizontal dimensionless velocity component, u*, (b) dimensionless microrotation, N*, and (c) dimensionless temperature, h*, profiles at the vertical mid-plan (x = 0) of the cavity and (d) Nu versus Ra for K = 5, k ¼ 0:1, n = 0 and B = 0.01, 0.1, 1.

Nu ¼

1 1  A3 Ra C 2

ð27Þ

where A3 ¼ 2C 1 sinhðX0 Þ þ C 2 =960 þ C 3 =24  C 4 . The above analytical model is valid asymptotically in the limit of both A  1 and Pr  1. It follows that the resulting problem is now governed by only five parameters, namely the thermal Rayleigh number Ra, the materials parameters K, B, k and the microrotation boundary condition n. The effects of these parameters on the flows and heat transfer rates will be now discussed. Fig. 3a–c shows the influence of the vortex viscosity parameter K on the distribution of horizontal velocity u*, angular velocity N*, and profiles temperature h*, at the vertical mid-plan (x = 0) of the cavity for B = 1, k ¼ 0:1 and n = 0. The Nusselt number versus the Rayleigh number is depicted in Fig. 3d, for the same conditions. The curves depicted in these graphs are the prediction of the parallel flow approximation. The numerical solution of the full governing equations, depicted by dots, is seen to be in good agreement with the analytical model. In general K is a function of both the shape and the concentration of the microelements. Thus, for a given shape of the microelements, K is directly related to the

concentration of the micropolar solution, the concentration of microconstituents increasing with the increase of K. From Fig. 3a it is observed that, upon increasing parameter K, the intensity of the convective velocity u* is reduced as compared to the Newtonian fluid situation (K = 0). The presence of micropolar elements is responsible for the retardation of the convective flow. This follows from the fact that an increase of the value of the vortex viscosity K results in a larger resistance to the fluid motion. In fact, according to the present theory, it is found that u* ? 0 as the vortex viscosity parameter K ? 1. Fig. 3b illustrates the influence of parameter K on the angular velocity N*. Naturally, N* = 0 for K = 0, since no rotation can occur in the absence of micropolar elements (pure fluid situation). Upon increasing K the intensity of N*, at the center of the layer (y = 0), first increase. However, for K above unity, the value of N* (y = 0) starts to decrease. The fact that a critical value of the viscosity ratio exists on N* has already been reported in the past by Kim and Kim [28], while studying plane Couette micropolar flows. For this flow configuration it was found that the critical value of N* occurs at about K = 0.3. Furthermore, Fig. 3b indicates the formation of two layers near the solid boundaries, with negative values of N*, sandwiching a core layer in which

2756

Z. Alloui, P. Vasseur / International Journal of Heat and Mass Transfer 53 (2010) 2750–2759

(a)

0.5

(c) 10 1

0.5

Analytical

Analytical

Numerical

Numerical

λ = 0.1

y

0

y

0

λ = 0.1 1 10 -0.5 -0.002

0

u (b)

-0.5

0.002 *

0.5

(d)

-0.0001

0

0.0001

θ

*

3.5

Analytical Numerical

3

2.5

y

10

0

1

λ = 0.1

Nu

λ 0.1 1 10

2

Analytical Numerical

1.5

-0.5 -0.002

0

0.002

0.004

N*

1

104

105

Ra

Fig. 6. Distribution of (a) horizontal dimensionless velocity component, u*, (b) dimensionless microrotation, N*, and (c) dimensionless temperature, h*, profiles at the vertical mid-plan (x = 0) of the cavity and (d) Nu versus Ra for K = 5, B = 1, n = 0 and k ¼ 0:1; 1; 10.

N* is positive. The influence of K on the convective heat transfer is illustrated in Fig. 3c and d. Since a large portion of the energy of the system is used in producing gyrational velocity of the micropolar fluid the convective heat transfer is reduced, upon increasing K, as indicated by the temperature distributions of Fig. 3c.Thus, for K ? 1 the temperature profile is a vertical line corresponding to a pure conduction situation. Fig. 3d illustrates the influence of the vortex viscosity K on the variation of the Nusselt number, Nu, versus the Rayleigh number, Ra. For a given value of K it is seen that below a critical Rayleigh which depends upon the governing parameters, Eq. (26), the fluid is at rest and the heat transfer purely conductive (Nu = 1). Above the critical Rayleigh number, the Nusselt number is seen to increases first monotonically with Ra to reach asymptotically a constant value which depends upon K. The fact that for large Rayleigh numbers the Nusselt number reaches a constant is due to the particular thermal boundaries conditions considered here, namely Neumann conditions. For this situation it is well known that, upon increasing Ra, the temperatures at each point of the system decrease continually. However, at sufficiently large values of Ra, a regime is reach where the temper-

atures difference between any two given points within the cavity, and thus the Nusselt number, Eq. (11), remains constant. As discussed above, since an increase of K inhibits the convective flow it follows that the resulting heat transfer is reduced as it can be observed from the curves presented in Fig. 3d. A similar trend was reported in the past by many authors, and in particular by Zadravec et al. [12], for the case of a cavity differentially heated from the sides. In the following paragraphs the effects of parameters n, B and k on the convective heat transfer will be discussed. Since the values of K considered in literature covers the range 0 6 K 6 10 a mean value K = 5 was considered in this study. The influence of micro-gyration parameter n on the velocity u*, temperature h* and microrotation N* profiles for K = 5, B = 1 and k ¼ 0:1 is depicted in Fig. 4a–c. The corresponding Nusselt number Nu versus the Rayleigh number Ra is presented in Fig. 4d. Fig. 4a and b indicates that the intensity of the convective flow u* and that of the angular velocity N*, are minimum for n = 0. This particular value of n represents the case where the concentration of the microelements is sufficiently large so that the particles, close to the walls, are unable to rotate. Naturally, as expected, the correspond-

2757

Z. Alloui, P. Vasseur / International Journal of Heat and Mass Transfer 53 (2010) 2750–2759

ing heat transfer is also minimum for this particular boundary condition, as it can be seen from Fig. 4c and d. Upon increasing the value of n, the concentration of the solution becomes weaker such that the rotation of the particles near the walls are free to rotate. Thus, as n is augmented the microrotation term is augmented which induces enhancement of flow velocities. The resulting heat transfer is promoted as illustrated in Fig. 4c and d. Fig. 5a and b exemplifies the effect of micro-inertia parameter B on the velocity, and microrotation profiles for K = 5, k ¼ 0:1 and n = 0. The results indicate that, upon increasing parameter B, the intensity of the convective flow is slightly enhanced, as illustrated by Fig. 5a. The effect of B on the angular velocity N* is illustrated in Fig. 5b. Naturally, for B = 0, i.e. very large values of the micro inertia j, the angular velocity is nil. On the other hand, the increasing of B favours greatly the rotation of the microelements within the micropolar fluid. A similar trend was observed by Hsu et al. [9] while studying natural convection of micropolar fluids in an enclosure with heat sources. According to Danesh et al. [29] this is due to favourable non-Newtonian effects. From Fig. 5c it is found that the distortion of the temperature profiles is slightly increased as the value of B is made larger. This follows from the fact that, as discussed above, the flow intensity is promoted upon increasing B. The influence of parameter B on the heat transfer depends on the Rayleigh number as it can be seen from Fig. 5d. Thus, when Ra is relatively small it is found that the Nusselt number, Nu, increases as the value of B is made larger, in agreement with the results of Hsu et al. [9]. However, the present analytical and numerical results indicate that this trend is reversed when Ra is above a given value that depends upon the controlling parameters. The influence of the material parameter k, on the velocity u* and microrotation N* profiles for K = 5, B = 1 and n = 0, is depicted in Fig. 6a and b, respectively. k is a viscosity ratio, namely the spingradient viscosity over the dynamic viscosity. Thus, decreasing the dynamic viscosity, i.e., increasing k, results in an increase in the friction between the fluid layers. It follows that increasing the spin-gradient viscosity decreases the velocity of the convective flow, Fig. 6a. Consequently, the distortion of the vertical temperature distributions is reduced, upon decreasing k, as exemplified in Fig. 6c. On the other hand, Fig. 6b shows a decrease of the micrororotation with an increase of k. This follows from the fact that increasing k, i.e. decreasing the dynamic viscosity, reduces the friction between the fluid layer such that the rotation N* is inhibited. These results are in agreement with those of Hsu and Hsu [9] who reported that the value of the maximum rotation decreases sharply for k < 1, continues gradually for k > 1 and the approaches zero for large values of k ðk ! 1Þ. On the other hand, it is observed from Fig. 6d, that the effect of k on the heat transfer does not show consistent variations and depends upon the value of the Rayleigh number. Thus, when the convective motion is relatively weak (Ra < 4  104) Fig. 6d indicates that an increase of the spin-gradient viscosity parameter decreases the heat transfer rate. A similar trend was reported by Hsu et al. [9], while studying natural convection of micropolar fluids in a rectangular enclosure heated from the bottom. However, for larger values of the Rayleigh number (Ra > 4  104), it is observed from Fig. 6d that this trend is reversed. As reported by Hsu and Chen [8], the increase in k increases the couple stress of the fluid such that the convective heat transfer is higher at large spin-gradient viscosity. 5. Linear stability analysis The stability to small perturbations from the quiescent state (wC = NC = 0 and TC = 1  y) of the physical situation described by (6)–(8) is examined now. In order to do so, it is convenient to reb ¼ N  N C and ^ ¼ww , N write the governing equations using w C

b ¼ T  T C . As usual, the perturbed solution is assumed to have T the following functional form:

9 ptþikx ^ x; yÞ ¼ wðyÞe ~ wðt; > = g^ ðt; x; yÞ ¼ g~ ðyÞeptþikx > ^hðt; x; yÞ ¼ ~hðyÞeptþikx ;

ð28Þ

where p = r + ix is the complex amplification rate of the perturbation, k is the real wave number and x the frequency. Introducing (28) into (6)–(8) and neglecting second higher-order non-linear terms yields the following linear system: 2 2 ~ 2 ~ Prð1 þ KÞðD2  k ÞðD2  k Þw þ Pr KðD2  k Þg 2 ~ 2 ~  ik Pr Ra h ¼ pðD  k Þw

ð29Þ

2

~ ¼ pg ~  Pr BK½2g ~ þ ðD2  k2 Þw ~ Pr kðD  k Þg 2 ~ 2 ~ ¼ ph~ ðD  k Þh  ikw 2

ð30Þ ð31Þ

The corresponding boundary conditions are

y ¼ 1=2;

~ ¼ Dw ~ ¼ 0; w

~ g~ ¼ nD2 w; D~h ¼ 0

ð32Þ

where D = d/dy. The perturbed state equations (29)–(31) with the boundary conditions (32) may be written in a compact matrix form as

L1 ðkÞY ¼ pL2 ðkÞY

ð33Þ

~ g ~; ~ where Y ¼ ½w; h is a three-component vector of the perturbation and L1(k) and L2(k) are two linear differential operators that depend on the control parameters Ra, Pr, K, B, k and n. The set of Eq. (33) is solved using a finite differences scheme. The system is discretized using a fourth-order scheme in the domain between y = 1/2 and y = 1/2. For M computational points, the resulting discrete system has 3M eigenvalues that can be found using a standard IMSL subroutine such as DGVCCG. The precision of the value of the critical Rayleigh number predicted by the present numerical procedure depends on the number of points M taken in the y direction. Numerical tests, using various values of M were done for the same conditions in order to determine the best compromise between accuracy of the results and computing time. The validation of the present stability analysis was made for the case of a horizontal pure fluid layer heated from the bottom by a constant heat flux. For this situation, investigated in the past by Sparrow et al. [30], the onset of motion occurs at zero wave number, at a critical Rayleigh number Rac = 720. Typical results are presented in Table 1. Based on these results, the value M = 40 was adopted for this study. Fig. 7 a–c presents series of streamline (top), isotherm (middle) and microrotation (bottom) patterns at the onset of motion, as predicted by the linear stability analysis, for Ac = 6, K = 10, B = 1 and k ¼ 0:1 and various values of the micro-gyration parameter n. For n = 0, Fig. 7a, the critical Rayleigh number occurs at Rac = 5965, while it occurs at Rac = 4609 for n = 0.5, Fig. 7b, and Rac = 2705, for n = 1, Fig. 7c. A bird eye view of Fig. 7a–c indicates that the effect of n on the streamline and isotherm patterns is relatively weak. The flow structures consist in a single centro-symmetrical cell while the corresponding temperature fields indicate a constant horizontal stratification. However, the microrotation fields are found to be affected by parameter n. As depicted from the graphs,

Table 1 Comparison between the present numerical solution of the linear stability analysis and the results of Sparrow et al. [30].

Rac

Sparrow et al. [30]

Present M ¼ 10

Present M ¼ 20

Present M ¼ 30

Present M ¼ 40

Present M ¼ 50

720

729.66

720.63

720.11

720.04

720.02

2758

Z. Alloui, P. Vasseur / International Journal of Heat and Mass Transfer 53 (2010) 2750–2759

the microrotation patterns consist in superposition two shallow layers adjacent to the horizontal boundaries of the cavity, sandwiching a thicker layer in the core region. From the numerically determined microrotation vertical velocity profiles (not presented here), it is found that the thickness of the tree layers is not affected by the value of n. It is interesting to observe that the structure of the microrotation pattern predicted by the linear stability theory, Fig. 7, is quantitatively in agreement with the results predicted by the parallel flow theory and the numerical solution, Figs. 4d–6d. Fig. 8a–c is a plot of the critical Rayleigh number Rac versus the vortex viscosity parameter K for various values of the material parameters, B, k and n. In these graphs, the solid lines are the prediction of the parallel flow approximation which is observed to be in excellent agreement with the numerical results obtained from the linear stability theory, depicted by solid squares. It can be demonstrated that in the limit K ? 0, Eq. (26) yields the following relationship:

(a)

6000

n 0 0.5 1

4000

Analytical Linear stability

Rac 2000

720

120 2 20 Rac ¼ 720ð1 þ KÞ  BK ð1 þ 7nÞ þ 2 B2 K 3 ð1 þ 9nÞ þ OðK 7=2 Þ 7k 7k ð34Þ such that, in the case of a pure fluid (K = 0), the classical value Rac = 720 is recovered.

10

-2

10

-1

10

0

10

1

10

0

10

1

10

0

10

1

K

(b)

6000

B 0.01 0.1 1

4000

Analytical Linear stability

ψˆ

Rac θˆ

2000

(a)

ηˆ 720 10

ψˆ θˆ

-2

10

-1

K

(c) (b)

6000

λ 0.1 1 10

4000

Analytical Linear stability

ηˆ

Rac 2000

ψˆ θˆ

(c) 720

ηˆ

10

-2

10

-1

K Fig. 7. Streamline (top), isotherm (middle) and microrotation (bottom) patterns at the onset of convection for K = 10, B = 1, k ¼ 0:1, Ac = 6; (a) n = 0, Rac = 5965; (b) n = 0.5, Rac = 4609; (c) n = 1, Rac = 2705.

Fig. 8. Variation of critical Rayleigh number with parameter K for (a) B = 1, k ¼ 0:1 and n = 0, 0.5, 1, (b) k ¼ 0:01, n = 0 and B = 0.01, 0.1, 1, (c) B = 1, n = 0 and k ¼ 0:1; 1; 10.

Z. Alloui, P. Vasseur / International Journal of Heat and Mass Transfer 53 (2010) 2750–2759

In general, Fig. 8a–c indicates that an increase in K is responsible of delaying the onset of motion, i.e. the critical Rayleigh number Rac, and this independently of the values of the other material parameters. As already discussed this is due to the increase of the total viscosity upon increasing K. The influence of the dimensionless micro-gyration parameter n is depicted in Fig. 8a for B = 1 and k ¼ 0:1. The results clearly indicate that, for a given value of K, the critical Rayleigh number Rac decreases as the value of n is made larger. Since a decrease of n implies an increasing concentration of the microelements a greater part of the energy of the system is required in developing gyrational velocities of the fluid and, as a result, the onset of convection is delayed. Fig. 8b illustrates the influence of the micro-inertia parameter B for n = 0 and k ¼ 0:01. The results indicate that, for a given K, a higher critical Rayleigh number is required as the value of B is made smaller, i.e. that the micro-inertia density is increased. Similarly, Fig. 8c indicates that the increase in k, which cause an enhancement of the couple stress of the fluid, is to stabilize the system such that a higher critical Rayleigh number is required for the onset of motion. However, as depicted by Fig. 8c, the influence of this parameter is relatively small. 6. Conclusions In this paper an analytical and numerical study of natural convection in a shallow rectangular cavity filled with a micropolar fluid is carried out. Analytical solutions for finite-amplitude convection are derived on the basis of the parallel flow approximation for a shallow enclosure. A finite volume method is used to solve numerically the present problem. The governing parameters of the problem are the thermal Rayleigh number, Ra, the materials parameters K, B, k and the microrotation boundary condition n. The main conclusions of the present analysis are as follows: (1) An approximate analytical model has been derived to predict the flow, microrotation and heat transfer for the problem considered here. The model is found to be in excellent agreement with a numerical solution of the full governing equations. (2) The critical Rayleigh number for the onset of motion is predicted explicitly by the parallel flow theory, in terms of the governing parameters of the problem. The results are in agreement with the numerical results obtained on the basis of a linear stability analysis. (3) The effects of the characteristic parameters of micropolar fluids on the flows and heat transfer in a shallow cavity are found to be significant. However, depending on the strength of convection, i.e. Ra, the influence of B and k on Nu can be quite different.

2759

References [1] A.C. Eringen, Nonlinear theory of simple micro-elastic solids, Int. J. Eng. Sci. 2 (1964) 189–203. [2] A.C. Eringen, Theory of micropolar fluids, J. Math Mech. 16 (1966) 1–18. [3] A.C. Eringen, Non-local Polar Field Theory, Academic Press, New York, 1976. [4] T. Ariman, M.A. Turk, N.D. Sylvester, Microcontinuum fluid mechanics – a review, Int. J. Eng. Sci. 11 (1973) 905–930. [5] A.C. Eringen, Theory of thermomicrofluids, J. Math. Anal. Appl. 38 (1972) 480–496. [6] M. Balaram, V.U.K. Sastry, Micropolar free convection flow, Int. J. Heat Mass Transfer 16 (1973) 437–441. [7] G. Maiti, Convective heat transfer in micropolar fluid flow through a horizontal parallel plate channel, ZAMM 55 (1975) 105–111. [8] T.H. Hsu, C.K. Chen, Natural convection of micropolar fluids in rectangular enclosures, I, J. Eng. Sci. 34 (1996) 407–415. [9] T.H. Hsu, P.T. Hsu, S.Y. Tsai, Natural convection of micropolar fluids in an enclosure with heat sources, Int. J. Heat Mass Transfer 40 (1997) 4239–4249. [10] O. Aydin, I. Pop, Natural convection from a discrete heater in enclosures with micropolar fluid, Int. J. Eng. Sci. 43 (2005) 1409–1418. [11] O. Aydin, I. Pop, Natural convection in differentially heated enclosure filled with a micropolar fluid, Int. J. Eng. Sci. 46 (2007) 963–969. [12] M. Zadravec, M. Hribersek, L. Skerget, Natural convection of micropolar fluid in enclosure with boundary element method, Eng. Anal. Boundary Elem. 33 (2009) 485–492. [13] K.V. Rama Rao, Onset of instability in a heat conducting micropolar fluid layer, Acta Mech. 32 (1974) 79–84. [14] V.U.K. Sastry, V. Ramamohan, Numerical study of thermal instability of a rotating micropolar fluid layer, Int. J. Eng. Sci. 21 (1983) 449–461. [15] Y. Qin, P.N. Kaloni, A thermal instability problem in a rotating micropolar fluid, Int. J. Eng. Sci. 30 (1992) 1117–1126. [16] Y.N. Murty, Effect of throughflow and magnetic field on Bénard convection in micropolar fluid, Acta Mech. 150 (2001) 11–21. [17] N. Rudraih, V. Ramachandramurty, O.P. Chandna, Effects of magnetic field and non-uniform temperature gradient on Marangoni convection, Int. J. Heat Mass Transfer 28 (1985) 1621–1624. [18] A. Abraham, Rayleigh–Benard convection in a micropolar ferromagnetic fluid, Int. J. Eng. Sci. 40 (2002) 449–460. [19] J.F. Melviana, M.A. Norihan, N.S. Mohd, M.N. Roslinda, Effect of non-uniform temperature gradient on Marangoni convection in a Micropolar fluid, Eur. J. Sci. Res. 4 (2009) 612–620. [20] R. Idris, H. Othman, I. Hashim, On effect of non-uniform basic temperature gradient on Bénard–Marangoni convection in micropolar fluid, Int. Commun. Heat Mass Transfer 36 (2009) 255–258. [21] S.K. Jena, M.N. Mathur, Similarity solutions for laminar free convection flow of a thermomicropolar fluid past a nonisothermal vertical flate, Int. J. Eng. Sci. 19 (1981) 1431–1439. [22] G. Ahmadi, Self-similar solutions of incompressible micropolar boundary layer flow over a semi-infinite plate, Int. J. Eng. Sci. 14 (1976) 639–646. [23] J. Peddieson, An application of the micropolar fluid model to the calculation of a turbulent shear flow, Int. J. Eng. Sci. 10 (1972) 23–32. [24] S. Patankar, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC, 1980. [25] O. Trevisan, A. Bejan, Natural convection with combined heat and mass transfer, Int. J. Heat Mass Transfer 28 (1985) 1597–1611. [26] P. Vasseur, C.H. Wang, M. Sen, The Brinkman model for natural convection in a shallow porous cavity with uniform heat flux, Numer. Heat Transfer 15 (1989) 221–242. [27] A. Bejan, The boundary layer regime in a porous layer with uniform heat flux from side, Int. Heat Mass Transfer 26 (1983) 1339–1346. [28] Y.J. Kim, T.A. Kim, A study on the plane Couette flow using micropolar fluid theory, K.S.M.E. Int. J. 18 (3) (2004) 491–498. [29] R.A. Damseh, T.A. Al-Azab, B.A. Shannak, M. Al Husein, Unsteady natural convection heat transfer of micropolar fluid over a vertical surface with constant heat flux, Turkish J. Eng. Env. Sci. 31 (2007) 225–233. [30] E.M. Sparrow, R.J. Goldstein, V.K. Jonsson, Thermal instability in a horizontal fluid layer: effect of boundary conditions and nonlinear temperature profile, J. Fluid Mech. 18 (1964) 13–528.