Numerical simulation of free convection within an eccentric annulus filled with micropolar fluid using spectral method

Applied Mathematics and Computation 219 (2013) 5397–5409

Contents lists available at SciVerse ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Numerical simulation of free convection within an eccentric annulus filled with micropolar fluid using spectral method F.M. Mahfouz Mechanical Power Department, Menoufia University, Egypt

a r t i c l e

i n f o

Keywords: Natural convection Fourier Spectral Method Micropolar fluid Eccentric annulus

a b s t r a c t In this paper, the problem of natural convection in a circular eccentric annulus filled with micropolar fluid has been numerically investigated using Fourier Spectral Method. The annulus inner wall is heated and maintained at constant temperature while the outer wall is cooled and kept at constant temperature. The full governing equations of momentum, angular momentum and energy have been solved to give the details of flow and thermal fields. The heat convection process in the annulus is mainly controlled by Rayleigh number Ra, Prandtl number Pr, radius ratio Rr, eccentricity, e and material parameters of micropolar fluid. The material parameters are dimensionless spin gradient viscosity k, dimensionless micro-inertia density B and dimensionless vortex viscosity D. This study considers Ra up to 5  104 and the eccentricity is varied between 0.6 and +0.6 while dimensionless material parameters D of micropolar fluid is considered between 0 and 10. Both Pr, and Rr are fixed at 0.7, 2.6, respectively while B and k are assigned the value of 1. The study considers the effect of controlling parameters on flow and thermal fields with emphasis on the effect of these parameters on mean Nusselt number. The study has shown that for certain controlling parameters the heat transfer in the annulus is minimum at a certain eccentricity. The study also has shown that as parameter D increases the rate of heat transfer through the annulus decreases. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Micropolar fluids are fluids with microstructure. They represent fluids consisting of rigid, randomly oriented or spherical particles suspended in a viscous medium (such as polymeric suspensions, animal blood, liquid crystals), where the deformation of fluid particles is ignored. These fluids are of importance in the chemical, pharmaceutical and food industries. In order to study such type of fluids Eringen [1] introduced the theory of micropolar fluids. In this theory, the microscopic effects arising from local structure and microrotation of the fluid elements were taken in consideration. In this case many classical concepts such as the symmetry of the stress tensor or absence of couple stresses are no longer existed. In practice, the theory of micropolar fluid requires solving an additional transport equation representing the principle of conservation of local angular momentum. To consider heat convection in micropolar fluids, the theory of thermomicropolar fluids was developed by Eringen [2] by extending the theory of micropolar fluids. An extensive reviews of the theory and applications of micropolar fluids can be found in the review articles by Ariman et al. [3,4] and the books by Lukaszewicz [5] and Eringen [6]. Buoyancy driven flow and associated heat convection in an annular enclosure has long been investigated because of its pertinence to many practical engineering applications. These applications include, thermal storage systems, solar collectorreceivers, compressed gas insulating high-voltage electric transmission cables, cooling systems in nuclear reactors, and

E-mail address: [email protected] 0096-3003/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.amc.2012.11.038

5398

F.M. Mahfouz / Applied Mathematics and Computation 219 (2013) 5397–5409

Nomenclature B c d E e fn, Fn Fb g gn, Gn h hn, Hn J k Keq Keq Kv j L M Nu, Nu Nuconc Nucond Pr Rr Ra RaL t T x0 ; y0 x, y R⁄

dimensionless micro-inertia density half distance between the two poles of Bi-polar coordinates c/Ri eccentricity dimensionless eccentricity = E/(Ro  Ri ) Fourier coefficients buoyancy force gravitational acceleration Fourier coefficients local heat transfer coefficient Fourier coefficients Jacobian of coordinate transformation thermal conductivity equivalent conductivity (Nu=Nuconc ) mean equivalent conductivity (Nu=Nuconc ) vortex viscosity micro-inertia density radius difference = (Ro  Ri ) dimensionless microrotation local and mean Nusselt numbers conduction Nusselt number in case of concentric annulus (=2/LnRr) conduction Nusselt number in case of eccentric annulus Prandtl number (¼ m=a) radius ratio ð¼ Ro =Ri Þ Rayleigh number based on 2Ri (=g bð2Ri Þ3 ðDT=am) Rayleigh number based on L(=gbðLÞ3 DT=am) dimensionless time temperature Cartesian coordinates dimensionless Cartesian coordinates the dimensionless distance along the gap between the two cylinders = ðrðhÞ  Ri Þ=ðRðhÞ  Ri Þ

Greek symbols a thermal diffusivity b coefficient of thermal expansion DT temperature difference (=T i  T o ) g; n bi-polar coordinates / dimensionless temperature c spin gradient viscosity m kinematics viscosity k dimensionless spin gradient viscosity q density r microrotation in x0  y0 plane s time dimensional and dimensionless stream functions w0 , w f0 , f dimensional and dimensionless vorticity Subscripts i at the inner tube surface o at the outer tube surface

many others. The previous experimental and numerical investigations on thermal convection within an annulus that formed between two concentric/eccentric cylinders are numerous. Most of the these studies [7–16] focused on natural convection in case of concentric/eccentric annulus filled with Newtonian fluids. Kuehn and Goldstein [7] has provided a thorough literature review and conducted somehow detailed study in case of concentric annulus with fixed radius ratio of 2.6. In their study, they carried out both numerical simulations using finite difference method and experimental study using Mach–Zehnder interferometer. Their heat transfer results were presented in terms of the Nusselt number as well as the equivalent conductivity. The case of isothermal heating of concentric annulus has been further considered by a number of authors [8–10], to mention a few. In case of natural convection within

F.M. Mahfouz / Applied Mathematics and Computation 219 (2013) 5397–5409

5399

isothermally heated eccentric annuli, an early study has been conducted by Kuehn and Goldstein [11] who provided somehow detailed experimental results, using the experimental setup they have used in [7]. Since then, several experimental and theoretical works [12–16] have been conducted to explore further the effect of different parameters such as the Rayleigh number, the eccentricity and the Newtonian fluid properties on heat convection within such annuli. The study of both flow and heat convection problems related to micropolar fluids has received a considerable interest. The previous works up to 1973 was compiled by Arimen et al. [3]. Since then the studies related to heat convection in micropolar fluids have focused mainly on unconfined buoyancy driven flows over simple geometry such as flat plates [17–19] and regular surfaces [20–23]. While only a few studies have considered free convection of micropolar fluid flow in confined enclosures. Hsu et al. [24] considered the free convection of micropolar fluids in a rectangular enclosure with heat sources. Chiu et al. [25] investigated the transient free convection of micropolar fluids in concentric spherical annuli while Chen [26] examined the transient natural convection of micropolar fluids between two concentric and vertically eccentric spheres. Aydin and Pop. [27] studied natural convection from a discrete heater in enclosures filled with a micropolar fluid. The main objective of this study is to analyze the free convection in a circular eccentric annulus filled with micropolar fluid. The heat transfer from the surface of a cylindrical cable to the surrounding eccentric enclosure that is filled with a dusty fluid represents an approximate practical example of isothermal heating of an eccentric annulus filled with micropolar fluid. Such a problem, to the best of author’s knowledge, has not yet been investigated. The analysis is based on the solution of full governing equations using the bi-polar coordinate transformation along with Fourier Spectral Method. The buoyancy driven flow is assumed laminar and two dimensional. 2. Problem formulation Fig. 1 shows the physical domain under consideration which is the annulus formed between two long horizontal eccentric cylinders. The outer cylinder has radius Ro , while the inner cylinder has radius Ri , and is placed in such away the vertical line passing through its center is passing through the center of outer cylinder. The eccentricity, E between the two cylinders is taken positive if the center of the inner cylinder above that of the outer cylinder and negative if below it. The annulus contains incompressible, micropolar fluid and is heated through its inner cylinder which is maintained at constant wall temperature, T i , while the outer wall is cooled and kept at constant temperature T o higher than that of the cooling ambient. The resulting buoyancy driven flow in the annulus is assumed to be laminar and two dimensional. The density of the fluid is considered constant except in the buoyancy term where it varies according to Buossensq approximation. 2.1. The governing equations

x

After impulsively heating the annulus through its inner wall, a buoyancy driven flow is initiated. The induced steady state flow is reached after a period of time during which the time dependent governing equations have to be solved. The time

ξο

η= 0

g

Fb Ro

ξ

ξ E

η

η=π

η Ri

y

c

c

Fig. 1. Physical domain and coordinate system.

5400

F.M. Mahfouz / Applied Mathematics and Computation 219 (2013) 5397–5409

dependent conservation equations of mass, linear and angular momentum and energy can be written in Cartesian coordinates in terms of stream function, vorticity, microrotation and temperature as:

@f0 ¼ @s



mþ

Kv

q



r2 f 0 þ

  ½w0 ; f0  1 @F x0 2 ; K r r þ  v ½x0 ; y0  q @y0

ð1Þ

f0 ¼ r2 w;0

ð2Þ

@r c ½w0 ; r K v 0 ¼ r2 r þ 0 0 þ ðf  2rÞ; ½x ; y  qj @ s qj

ð3Þ

@T ½w0 ; T ¼ ar2 T þ 0 0 ; @s ½x ; y 

ð4Þ

½w; s @w @s @w @s   ; ½x; y @x @y @y @x

s  ½f0 ; r; T;

ð5Þ

where s is the time, w0 is the stream function, f0 is the vorticity, r is microrotation in x0  y0 plane, T is the temperature, q is the density, m and a are the momentum and thermal diffusivities respectively. F x0 ¼ qgbðT  T o Þ is the net upward driving force, Fb. Kv, c and j are material parameters. 2.2. Initial and boundary conditions At time t = 0, the inner wall is assumed to be impulsively heated to constant wall temperature while the outer wall is kept at fixed temperature. The micropolar fluid in the annulus is initially at rest and at the same temperature of the outer wall. The hydrodynamics boundary conditions at the inner and outer walls are the no-slip condition and impermeability. The thermal boundary condition at the inner wall is the constant wall temperature (T = Ti) while the temperature of outer wall is kept constant at To. The steady state boundary conditions can be expressed mathematically as follows: -On the inner wall -On the inner wall

w0 ¼

@w0 ¼ 0; @x0

@w0 ¼ 0 and T ¼ T i : @y0

ð5aÞ

Sn is the normal direction to the surface. -On the outer wall

w0 ¼

@w0 ¼ 0; @x0

@w0 ¼ 0 and T ¼ T o : @y0

ð5bÞ

The boundary conditions for microrotation at inner and outer walls are taken as r ¼ qf0w where 0 6 q 6 1 .Many authors considered the case of q = 0 and some other considered q = 1/2 while the case of q = 1 has been considered to simulate turbulent flows. In this study the author recommends the value of q = 0 because this value represents the plausible no spin condition at the walls. Moreover, the boundary conditions of r ¼ 0 (or q = 0) at walls has been considered in a good number of previous studies [28–31]. 2.3. Coordinate transformation In order to set the inner and outer walls boundary conditions properly without interpolation, a coordinate system that makes the boundary lines coincide with coordinates lines should be used. In this problem the bipolar coordinates n; g defined as x0 ¼ c sinh n=ðcosh n  cos gÞ, y0 ¼ csing=ðcosh n  cos gÞ is adopted. In bipolar coordinates, after introducing the dimensionless variables x ¼ x0 =Ri , y ¼ y0 =Ri , d ¼ c=Ri , t ¼ sa=R2i , w ¼ w0 =a, f ¼ f0 R2i =a, D = Kv/mq, k ¼ c=mqR2i , B = j=R2i and / ¼ ðT  T o Þ=DT, the governing equations (1)–(4) read

J

  @f ½w; f JPrRa @/ @/ ; ¼ Prð1 þ DÞr2 f0 þ PrDr2 M þ ðcosh ncosg  1Þ  sinh nsing þ @t ½n; g 8d @g @n

Jf ¼ r2 w

ð6Þ ð7Þ

J

@M Prk 2 ½w; M DPr ¼ Jðf þ 2MÞ; r Mþ  @t B ½n; g B

ð8Þ

J

@/ ½w; / ¼ r2 / þ ; @t ½n; g

ð9Þ

F.M. Mahfouz / Applied Mathematics and Computation 219 (2013) 5397–5409

5401

where J ¼ xn yg  xg yn is the determinant of the Jacobian of coordinate transformation matrix, M ¼ rR2i =a is dimensionless microrotation, Pr ¼ m=a is the Prandtl number and Ra ¼ gbð2Ri Þ3 DT=ma is the Rayleigh number. Some of the previous studies have defined the Rayleigh number as RaL ¼ gbL3 DT=ma, L ¼ Ro  Ri . The RaL is considered in case of comparison with the previous studies. In the new coordinates, the steady boundary conditions (5) can now be expressed as: -On the inner tube surface, n ¼ ni

w¼

@w ¼ 0; @n

@w ¼ 0; @g

M ¼ qf and / ¼ 1:

ð10aÞ

-On the outer tube surface, n ¼ no

w¼

@w ¼ 0; @n

@w ¼ 0; @g

M ¼ qf and / ¼ 0;

ð10bÞ

where ni and no are negative constants define the inner and outer tube surfaces 1 1 (ni ¼ sinh ðdÞ, no ¼ sinh ðd=RrÞ and Rr ¼ Ro =Ri is the radius ratio and e ¼ E=ðRo  Ri Þ is the dimensionless eccentricity. 3. Method of solution The method of solution is based on Fourier Spectral Method in which the dependent variables stream function w, vorticity f, microrotation M and temperature / are approximated using Fourier series expansion. The approach is similar to that used by Badr and Dennis [32] and Mahfouz and Badr [33]. Considering the symmetry of both flow and temperature fields around the vertical axis. The dimensionless stream function, vorticity, microrotation and temperature can be expanded as: N X fn ðn; tÞ sinðngÞ;

wðn; g; tÞ ¼

ð11aÞ

n¼1

fðn; g; tÞ ¼

N X g n ðn; tÞ sinðngÞ;

ð11bÞ

n¼1

Mðn; g; tÞ ¼

N X r n ðn; tÞ sinðngÞ;

ð11cÞ

n¼1

/ðn; g; tÞ ¼

N X 1 Ho ðn; tÞ þ Hn ðn; tÞ cosðngÞ; 2 n¼1

ð11dÞ

where N is the number of terms considered in the Fourier series and fn, gn, Ho and Hn are the Fourier coefficients. Substitution of w, f, M and / defined in (11a)–(11d) in Eqs. (6)–(9) and integrating these equations (after multiplying each by 1, sinðngÞ and cosðngÞ respectively) from 0 to 2p results in the following set of differential equations:

@ 2 fn @n2

 n2 fn ¼ ao g n þ

N 1X g ðajmnj  anþm Þ; 2 1 m

ð12Þ

N @g 1X @g m @ 2 gn  n2 g n ðajmnj  anþm Þ ¼ Prð1 þ DÞ ao n þ 2 1 @t @t @n2

ao

N @r n 1 X @r m Prk @ 2 r n  n2 rn þ ðajmnj  anþm Þ ¼ B @t 2 1 @t @n2

ao

N @Ho X @Hm @ 2 Ho am þ Z o ðn; tÞ; þ ¼ @t @t @n2 1

N @Hn 1 @Ho 1X @Hm þ an þ ðanþm  ajnmj Þ ¼ ao 2 @t 2 1 @t @t

! þ Sn ðn; tÞ;

ð13Þ

! þ K n ðn; tÞ;

ð14Þ

ð15Þ @ 2 Hn @n2

! 2

 n Hn

þ Z n ðn; tÞ;

ð16Þ

P where ao ; an are the coefficients that result from expanding the term J in the form ðao þ 1 n¼1 an cos ngÞ. Sn, Kn, Zo and Zn are all easily identifiable functions. Eqs. (12)–(16) define a system of differential equations that should be solved simultaneously at

5402

F.M. Mahfouz / Applied Mathematics and Computation 219 (2013) 5397–5409

every time. The steady boundary conditions for all Fourier functions presented in Eqs. (12)–(16) are deduced from (10) and can be expressed as: -At n = ni

fn ¼ @fn =@n ¼ 0;

rn ¼ 0;

Ho ¼ 2 and Hn ¼ 0:

ð17aÞ

rn ¼ 0;

Ho ¼ 0 and Hn ¼ 0:

ð17bÞ

-At n = no

fn ¼ @fn =@n ¼ 0;

In order to get the vorticity distribution at inner wall an integral condition is developed by integrating Eq. (13) (after multiplying both sides by enn) from n = ni to n = no and using boundary conditions (17). The resulting integral condition reads:

Z ni

no

! N 1X ao g n þ g ðajmnj  anþm Þ enn dn ¼ 0: 2 1 m

ð18Þ

The above integral condition is used to calculate the values of the functions gn on the inner tube surface from which the vorticity distribution along the wall is predicted. The vorticity distribution at outer wall is calculated at (n = no) from Eq. (12). The finite difference method is used to dicretize the ordinary differential equation (13) and partial differential equations (13)–(16). In discretization process all spatial derivatives are approximated using second order accurate central finite difference scheme. The first derivatives at the inner and outer walls are approximated by one side second order accurate forward (at inner wall) or backward (at outer wall) three point finite difference scheme. The discretization of time dependent equations (13)–(16) is handled using Crank–Nicolson scheme. The forward time step is taken very small to achieve the accuracy of the scheme. At every time step the coupled non-linear system of equations has been linearized through an iterative procedure, resulting in tri-diagonal system of equations that has been solved by Tri-Diagonal Matrix Algorithm, TDMA. More details about solution procedure is given by Badr and Dennis [32] and Mahfouz and Badr [33]. 3.1. Heat transfer parameters After obtaining the distribution of vorticity, stream function and temperature the flow and heat transfer characteristics can be easily determined. The vorticity distribution along the walls gives an indication about the velocity gradients as well as shear stresses along the walls. The stream function distribution helps find flow velocity components in n and g directions as V n ¼ J 1=2 @w=@ g; V g ¼ J 1=2 @w=@n. While the temperature distribution helps find the heat transfer results in terms of local and mean Nusselt numbers. The local and mean Nussult numbers are the direct indicator of local heat flux and total rate of heat transfer at the walls. The local Nusselt numbers at the inner and outer walls are defined as

Nui ¼ 2hi Ri =k;

Nuo ¼ 2ho Ro =k;

ð19Þ

where k is the fluid thermal conductivity, and hi and ho are the local heat transfer coefficients at inner and outer walls and are defined as:

hi ¼ qi =DT; ho ¼ qo =DT;

qi ¼ kð@T=@Sn Þi ; qo ¼ kð@T=@Sn Þo ;

Sn is the local direction of heat flux normal to the wall. From the above definitions one can deduce

  @/ Nui ¼ 2 J 1=2 ; @n i

  @/ Nuo ¼ 2Rr J 1=2 : @n o

ð20Þ

RP The mean Nusselt number is defined as Nu ¼ 1p 0 Nudp where P is the perimeter of the circular section. The mean Nusselt number at inner and outer walls can now be expressed as:

  @Ho Nui ¼  ; @n i

  @Ho Nuo ¼  : @n o

ð21Þ

Based on the above definition of mean Nusselt number at inner and outer walls it can be easily concluded that the ratio of mean Nusselt number at the inner and outer walls is equal to the ratio of the total rate of heat transfer from the two walls, that is Nui =Nuo ¼ Q_ i =Q_ o . Based on energy conservation within the annulus, it can be easily concluded that in the steady state condition the heat transfer rate into the annulus through the inner wall must be equal to the heat transfer rate rejected out the annulus through the outer wall. Thus when the mean Nussult number at outer wall approaches that at inner wall (within tolerance of 0.001) the steady state solution is reached and the calculation is terminated. In case of steady state pure conduction it is proved, through the analytical solution of energy equation, that the mean Nusselt number can be expressed as:

Nucond ¼ Nui ¼ Nuo ¼ 2=X io ; 1

where X io ¼ cosh

ððR2o

þ

R2i

2

 E Þ=2Ro Ri Þ

ð22Þ

5403

F.M. Mahfouz / Applied Mathematics and Computation 219 (2013) 5397–5409

Num. isotherms and ( streamlines) Present

Exp. isotherms Kuehn and Goldstein [11]

Num. isotherms and ( streamlines) Chung, et al. [34]

Fig. 2. Comparison of present numerical results with previous experimental and numerical results at RaL = 4.8  104, Rr = 2.6 and e = 0.652.

3.2. Accuracy of numerical solution The accuracy of the method of solution is first verified by comparing the present results with the most relevant results in the literature. Due to the absence of previous experimental and numerical results on the same present problem, the accuracy of the numerical solution is first assessed by comparing the results for the special case of Newtonian fluid (D = 0) with the

12 4 RaL= 4.8 x 10 , Pr=0.7 Rr=2.6, e=0.652

10

Present Exp. Kuehn and Goldstein [11)

8

Num. Chung et al.[34]

Keq 6 Outer wall

4

2

Inner wall

0 0

45

90

135

180

θ Fig. 3. Distribution of local equivalent conductivity along inner and outer walls at RaL = 4.8  104, e = +0.652.

Table 1 Effect of RaL and e on steady Keq at D = 0 and comparison with Badr [13].  (inner wall) Keq

RaL

a

e = 0.6

e = 0.3

e = 0.0

e = 0.3

e = 0.6

0

Badr [13] Present

1.235 1.234 (1.234)a

1.046 1.045 (1.045)a

1.001 1.001 (1.000)a

1.046 1.045 (1.045)a

1.235 1.234 (1.234)a

102

Badr [13] Present

1.239 1.237

1.047 1.047

1.001 1.001

1.047 1.046

1.237 1.236

103

Badr [13] Present

1.492 1.477

1.208 1.205

1.079 1.083

1.121 1.125

1.325 1.327

104

Badr [13] Present

2.315 2.175

2.220 2.061

2.030 1.982

1.995 1.828

1.965 1.825

 Analytical Keq.

5404

F.M. Mahfouz / Applied Mathematics and Computation 219 (2013) 5397–5409

available results in the literature. Fig. 2 shows a qualitative comparison of thermal field in terms of isotherms between the present numerical results and the experimental results of Kuehn and Goldstein [11] in case of RaL = 4.8  104 and e = 0.652. Shown also in the same figure are the numerical results of Chung et al. [34] for the same case. The complete pattern of experimental isotherms is presented in the figure while (due to symmetry) only half pattern of numerical isotherms is presented (and half pattern of the corresponding streamlines). The figure shows good qualitative agreement between the present numerical results and both numerical and experimental results. The corresponding distributions of equivalent conductivity for the aforementioned cases are shown in Fig. 3. The equivalent conductivity is defined as the actual heat flux divided by the heat flux that would occur by pure conduction. The figure shows good agreement with numerical results and reasonable agreement with experimental results. Table 1 shows another quantitative comparison between the present result for the mean steady state equivalent conductivity Keq with the numerical results of Badr [13] at Rr = 2.6, Pr = 0.7. Shown also in the table is the steady conduction values of Keq based on the analytical solution of energy equation. The table shows excellent agreement with the steady conduction values and very good agreement with numerical values of Badr [13]. It can be observed that the difference between present results and that of Badr [13] is very small at low Ra and is increasing with Ra reaching less than about 9% in case of Ra = 104 and e = 0.3. 4. Results and discussion The governing equations along with boundary conditions are solved in order to give the details of both flow and thermal fields. The main controlling parameters beside the classical ones (Ra and Pr) are the eccentricity (e), the radius ratio (Rr) and the dimensionless micropolar material parameters. These are the parameter D which characterizes the vortex viscosity, the parameter k which characterizes the spin gradient viscosity and the parameter B which characterizes the microinertia density. The close inspection of governing equations shows that the parameter D represents direct link between the flow field and microrotation field. The other two parameters k and B appear only in microrotation equation which gives expectation that the effect of these two parameters will not be as significant as the effect of parameter D. For this reason and for the sake of brevity only the influence of the vortex viscosity is considered. The parameter D is selected in the range from 0 to 10 while the other two material parameters k and B are assigned value of 1. These values for material parameters are adopted to satisfy the thermodynamics restrictions given by Eringen [2]. The Rayleigh number is varied up to 105 while the Prandtl number is fixed at 0.7. The eccentricity is considered between 0.6 and +0.6 while the radius ratio is fixed at Rr = 2.6. The numerical solution, as indicated before has marched in time till reaching the steady state solution. Fig. 4 shows the time development of Nu at the inner and outer surfaces at Ra = 5  104, Rr = 2.6 and e = 0.3 and for different values of D = 0, 2, 5 10. The time development of Nu at the inner and outer surfaces for the case of pure conduction regime is also shown in the same figure. The transient conduction solution is obtained by solving the energy equation with convection terms set to zeros. The figure shows that the time variation of Nu at inner and outer surfaces is generally similar to that for annulus filled with Newtonian fluid (D = 0). At the start of heating process, the heat is initially transferred from the inner wall through a very thin conduction layer, causing high heat transfer rate and so high Nui . The initial time stages are characterized by fast growth of the thermal layer and so quick decrease in heat transfer rate. In these stages, the conduction heat transfer dominates and the contribution of convection is almost negligible. This can be inferred from the coinciding of Nui t curve with that of conduction at these time stages. Once the conduction thermal layer thickens enough, the convection starts effectively to take

15 Ra=5 x 104, Pr=0.7 Rr=2.6, e= - 0.3

12

D=0 D=2 D=5

9

D==10

__ Nu

Ra=0 (cond.)

Inner wall

6

3 Ra=0 Outer wall

0 0.0

0.4

0.8

1.2

1.6

2.0

Time Fig. 4. The time development of mean Nusslet number of inner and outer walls of the annulus at Ra = 5  104, e = 0.3 and at different values of D.

F.M. Mahfouz / Applied Mathematics and Computation 219 (2013) 5397–5409

5405

part in heat transfer process till the steady state is slowly approached. Meanwhile, and after some time delay, as the thermal layer approaches the outer wall the heat transfer starts dissipating from the outer wall. As the time proceeds the rate of heat dissipated from the outer wall increases till reaching the steady state. At steady state condition the rate of heat transfer pumped into the annulus through the inner wall is equal to that dissipated from the outer wall, thus satisfying the steady energy balance. When such steady state is reached, the mean Nusselt numbers at the inner and outer walls, Nui and Nuo , will have the same value. One more observation related to heat transfer from the inner surface is that for a relatively smaller value of D, the transition from conduction dominated mode to convection dominated mode is accompanied with an overshot in heat transfer rate and so an overshot in Nui . The overshoot phenomenon is very obvious for D = 0 and D = 2 and less obvious for D = 5 and less further for D = 10. It can be observed from Fig. 4 that the mean steady heat transfer rates in the case of micropolar fluids (D = 2, 5, 10) is lower than that for Newtonian fluid (D = 0). This decrease is attributed to the increase of the effect of vortex viscosity which makes the fluid flow more viscous and so weakens the convection currents. In the initial stages where the conduction mode of heat transfer is dominant there is no effect for vortex viscosity and the heat transfer rates for micropolar fluid and Newtonian fluid are identical. As the convection mode starts developing the vortex viscosity of micropolar fluid enhances the flow viscosity and so weakens the convection currents which in turn decreases the heat transfer rate. The larger the value of vortex viscosity D the larger the flow viscosity and the lower the value of steady state Nu. Decreasing of heat transfer rates as the vortex viscosity D increases, is reported in the works of Mahfouz [21] and Hsu et al. [24]. Fig. 5a, b show both steady flow field, in terms of streamlines, and steady thermal field, in terms of isotherms, as well as steady vorticity and microrotation fields in case of Ra = 5  104, e = 0.3, k ¼ 1 and B = 1 and at the abovementioned four values of D = 0, 2, 5, 10. Since these fields are symmetrical about the vertical axis, only one half of each field is considered. Fig. 5a shows the patterns of streamlines (left) and isotherms (right) while Fig. 5b shows the patterns of equi-vorticity (right) and equi-microrotation (left). The streamlines pattern shows that the flow is formed of two main circulating eddies, one on the right side and another on the left side. Each eddy comprised of flow circulating upward along the inner wall and downward along the outer wall. The figure shows that as D increases the maximum value of stream function, wmax (written below each plot) decreases. The decrease of wmax gives an indication about the decrease of convection effects and the development of thicker both momentum and thermal layers. The isotherms pattern plotted in the same figure shows the temperature inversion within the annulus and the formation of the thermal plume at the top of the inner wall. At D = 0 the temperature inversion is steep and mushroom shaped thermal plume is clear while at D equal to 2, 5 and 10 the mushroom shape diminishes as D increases, indicating the weakness of convection currents and so the decrease of heat transfer rates. The equi-vorticitty and equi-microtation plot in Fig. 5b shows three regions of flow mean rotation (expressed by vorticity) and similar three regions for elemental rotation (in terms of microrotion). The three regions are separated by zero vorticity lines and zero microrotation lines (dashed lines in both sides). It can be observed that the flow adjacent to the walls (in the same side right or left) has the same sign of rotation and the same sign of microrotation while that in the core area has apposite rotation (and microrotation). The size of each region is influenced by the parameter D. As D increases the wall regions, especially for microrotation get thicker while core region gets smaller, reflecting the more viscous flow and the thicker boundary layer and the less effective convection. The steady state local Nusselt number distributions along inner and outer walls for the above cases is shown in Fig. 6. Since the thermal field is symmetrical about the vertical axis, only one half of Nu distribution is shown. It can be seen that the minimum Nu at inner wall takes place at h ¼ 0 while the maximum value takes place at h ¼ 180. At h ¼ 0 flow separates from the top point of inner wall where the temperature gradient is minimum due to prior heating of fluid during its upward motion. While at h ¼ 180 the cooled fluid that separated from the bottommost part of the outer wall (h ¼ 180) stagnates at the lower part of inner wall, creating relatively high temperature gradient and high heat flux and so Nu. The figure shows as D increases Nu decreases at all values of h along the inner wall which explains the smaller Nui as D increases as shown in Fig. 4. In case of outer wall the Nu assumes its maximum value at the h ¼ 0 where the already heated upward flow stagnates at h ¼ 0 resulting high temperature gradient and so high Nu. From topmost point (h ¼ 0) to bottommost point (h ¼ 180), the Nu decreases for all values of D, reaching a minimum at almost h ¼ 130 and then gradually increases to assume local maximum at h ¼ 180. Also it can be inferred from the distribution that as D increases the Nuo decreases. This is because as D increases the convection currents get weaker. The weakness of convection currents can be observed in Fig. 7 which shows the velocity distribution along the line of symmetry above the inner wall (h ¼ 0) and below the inner wall (h ¼ 180). The figure clearly shows that as D increases the levels of velocity along the line of symmetry decreases and so the convection effects. The effect of parameter D on steady state Nu at different Ra (including Ra = 0) number and eccentricity, e is presented in Table 2. It can be seen that as the parameter D increases at any fixed value of Ra, and (e) the Nu decreases. At low Ra this decrease is very small while at high Ra this decrease becomes significant reaching about 42% at Ra = 5  104, e = 0.3 as D increases from D = 0 to D = 10. The decrease of Nu, as explained before, is due to the increase of flow viscosity which weakens the convection effects. The effect of Ra can be also seen in the table where at any fixed value of (e) and D as Ra increases the Nu increases. This is quite expected since increasing of Ra means increasing of convection currents intensity and so increasing the heat transfer rate. The effect of eccentricity, e on Nu can be also inferred from Table 2. It is obvious that Nu assumes higher values at high eccentricity (positive or negative) and smaller values otherwise. When the two cylinders get very close to each other the rate of heat conduction increases and so increases the total rate of heat transfer and so Nu. This leads to a conclusion that there

5406

F.M. Mahfouz / Applied Mathematics and Computation 219 (2013) 5397–5409

Newtonian fluid M=0

ψ max = 24 ,

ζ max = 850

D=2

ψ max = 16 ,

M max = 18 ζ max = 310

D=5

ψ max = 11 ,

D=0

D=10

ψ max = 8 , a

M max = 16 ζ max = 180

M max = 13 ζ max = 100

b

Fig. 5. Patterns of (a) Isotherms (right) and streamlines (left) and (b) equi-vorticity (right) and equi-Microroatation (left) in case of Ra = 5  104, Pr = 0.7 and Rr = 2.6 and at D = 0, 2, 5 and 10.

F.M. Mahfouz / Applied Mathematics and Computation 219 (2013) 5397–5409

5407

20

Ra=5 x 104, Pr=0.7 Rr=2.6, e= - 0.3 16

D=0 D=2 D=5

12 Outer wall

D=10

Nu 8

4

Inner wall

0 0

30

60

90

120

150

180

θ Fig. 6. Local Nu distribution along inner and outer walls of the annulus at Ra = 5  104, e = 0.3 and at different values of D.

100

Ra=5 x 104, Pr=0.7 Rr=2.6, e= - 0.3 80

D=0 D=2

60

D=5 D=10

vξ 40 20

θ=0

0

θ=180 -20 0.0

0.2

0.4

0.6

0.8

1.0

R* Fig. 7. Velocity distribution along vertical centerline of the annulus at Ra = 5  104, e = 0.3 and at different values of D.

will be a value of the eccentricity at which the heat transfer in the annulus is minimum. In case of pure conduction (Ra = 0) the minimum occurs at zero eccentricity (i.e. concentric annulus) while in case of considering natural convection the table shows that the minimum is shifted towards the positive eccentricity. This is because the heat transfer rates for negative (e) is greater than those at same but positive (e). The reason behind such behavior may be roughly explained by considering the heat convection between two parallel horizontal plates under fixed temperature difference with one plate heated and the other plate cooled. It has been known that if the heated plate is facing up (and cooled plate is facing down) the convection contribution will be considerable and the heat transfer from the plate will be higher than that if the heated plate is facing down (and cooled plate is facing up). In this latter case the temperature gradient is not adequate for active convection and the heat transfer is mostly due to conduction. In case of vertically eccentric cylinders the upper part of the heated inner cylinder will be facing up and upper part of cold outer cylinder will be facing down while the lower part of the heated inner wall will be facing down and lower part of the cold outer wall will be facing up. In this configuration the heat transfer from the upper part of inner wall will be more as compared with that from the lower part. Therefore, the total heat transfer from the inner wall will be controlled by the effectiveness of heat convection in the gap between the upper part of the inner wall and the outer wall. In case of negative eccentricity the wide upper gap between the two walls helps enhance heat convection as a result of intensive motion of convection currents. While for same but positive eccentricity the narrow upper gap between the

5408

F.M. Mahfouz / Applied Mathematics and Computation 219 (2013) 5397–5409

Table 2 Effect of Ra, Rr, e, and D on steady mean Nusselt number at Rr = 2.6, B = 1 and k ¼ 1. Ra

0 5  102

5  103

5  104

D

– D=0 D=2 D=5 D = 10 D=0 D=2 D=5 D = 10 D=0 D=2 D=5 D = 10

 Nu e = 0.6

e = 0.3

e = 0.0

e = 0.3

e = 0.6

2.582 2.639 2.588 2.583 2.582 3.566 3.045 2.791 2.655 5.545 4.481 3.974 3.619

2.186 2.208 2.186 2.183 2.182 3.083 2.486 2.284 2.218 5.441 4.241 3.605 3.155

2.092 2.105 2.092 2.092 2.091 2.795 2.237 2.132 2.107 5.291 4.121 3.435 2.900

2.186 2.205 2.185 2.183 2.182 2.719 2.329 2.241 2.208 4.989 3.785 3.145 2.758

2.581 2.612 2.587 2.582 2.581 3.079 2.758 2.662 2.619 4.592 3.765 3.372 3.107

two walls and the inadequate temperature gradient below the inner wall hinder the convection currents and result in small convection contribution. 5. Conclusions The problem of natural convection in a circular eccentric annulus filled with micropolar fluid has been numerically investigated. The annulus inner wall is heated and maintained at constant temperature while the outer wall is cooled and kept at constant temperature. The heat convection process in the annulus depends mainly on Rayleigh number Ra, Prandtl number Pr, radius ratio Rr, eccentricity, e and material parameters of micropolar fluid. The material parameters are dimensionless spin gradient viscosity k, dimensionless micro-inertia density B and dimensionless vortex viscosity D. The study considers Ra up to 5  104 and the eccentricity is varied between 0.6 and +0.6 while dimensionless material parameters D of micropolar fluid is considered between 0 and 10. Both Pr, and Rr are fixed at 7, 2.6, respectively while B and k are assigned the value of 1. The study has confirmed that for certain controlling parameters as Ra increases the heat transfer rate (in terms of mean Nu) within the annulus increases. The study has also shown that for certain controlling parameters the heat transfer rate in the annulus is minimum at a certain eccentricity. The study has further concluded that the heat transfer is higher for negative eccentricity than that for same but positive eccentricity. Moreover, the study has shown that the rate of heat transfer through the annulus decreases as the dimensionless vortex vorticity, D increases. It could be also concluded that as D increases the level of velocity field decreases and so the convection effects. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]

A.C. Eringen, Theory of micropolar fluids, Journal of Mathematics and Mechanics 16 (1966) 1–18. A.C. Eringen, Theory of thermomicrofluids, J. Math. Anal. Appl. 38 (1972) 480–496. T. Ariman, M.A. Turk, N.D. Sylvester, Microcontinuum fluid mechanics, a review, Int. J. Eng. Sci. 11 (1973) 905–930. T. Ariman, M.A. Turk, N.D. Sylvester, Application of micro continuum fluid mechanics, International Journal of Enginering Science 12 (1974) 273–293. G. Lukaszewicz, Micropolar Fluids: Theory and Application, Birkhäuser, Boston, 1999. A.C. Eringen, Microcontinuum Field Theories II: Fluent Media, Springer, New York, 2001. T.H. Kuehn, R.J. Goldstein, An experimental and theoretical study of natural convection in the annulus between horizontal concentric cylinders, J. Fluid Mech. 74 (1976) 695–719. B. Farouk, S.I. Güçeri, Laminar and turbulent natural convection in the annulus between horizontal concentric cylinders, ASME J. Heat Transfer 104 (1982) 631–636. Y.T. Tsui, B. Trembaly, On transient natural convection heat transfer in the annulus between concentric, horizontal cylinders with isothermal surfaces, Int. J. Heat and Mass Transfer 27 (1983) 103–111. R. Kumar, Study of natural convection in horizontal annuli, Int. J. Heat Mass Transfer. 31 (1988) 1137–1148. N.H. Kuehn, R.J. Goldstein, An experimental study of natural convection heat transfer in concentric and eccentric horizontal cylindrical annuli, ASME J. Heat Transfer 100 (1978) 635–640. U. Projahn, H. Rieger, H. Beer, Numerical analysis of laminar natural convection between concentric and eccentric cylinders, Numer. Heat Transfer 4 (1981) 131–146. H.M. Badr, Study of laminar free convection between two eccentric horizontal tubes, Trans. Can. Soc. Mech. Eng. 7 (1983) 191–198. J. Prusa, L.S. Yao, Natural convection heat transfer between eccentric horizontal cylinders, ASME J. Heat Transfer 105 (1983) 108–116. D.W. Pepper, R.E. Cooper, Numerical solution of natural convection in eccentric annului, AIAA J. 21 (1983) 1331–1337. G. Guj, F. Stella, Natural convection in horizontal eccentric annuli: numerical study, Numer. Heat Transfer 27 (1995) 89–105. I.A. Hassanien, Mixed convection in micropolar boundary layer flow over a horizontal semi-infinite plate, ASME J. Fluids Eng. 118 (1996) 833–838. A.A. Mohammedien, R.S.R. Gorla, I.A. Hassanien, Mixed convection in an axisymmetric stagnation flow of micropolar fluid on a vertical cylinder, Acta Mecanica 114 (1996) 139–149. T.H. Hsu, S.P. How, Natural convection of micropolar fluids over a uniformly heated horizontal plate, Acta Mech. 155 (2002) 191–202. S. Bhattacharyya, I. Pop, Free convection from cylinders of elliptic cross section in micropolar fluids, Int. J. Eng. Sci. 34 (1996) 1301–1310. F.M. Mahfouz, Natural convection from an elliptic tube with major axis horizontal and placed in a micropolar fluid, Int. J. Heat Mass Transfer 47 (2004) 1413–1422.

F.M. Mahfouz / Applied Mathematics and Computation 219 (2013) 5397–5409

5409

[22] C.Y. Cheng, Free convection heat and mass transfer from a horizontal cylinder of elliptic cross section in micropolar fluids, Int. Commun. Heat Mass Transfer 33 (2006) 311–318. [23] C.Y. Cheng, Natural convection heat and mass transfer from a sphere in micropolar fluids with constant wall temperature and concentration, Int. Commun. Heat Mass Transfer 35 (2008) 750–755. [24] T. Hsu, P. Hsu, S. Tsai, Natural convection of micropolar fluids in an enclosure with heat sources, Int. J. Heat Mass Transfer 40 (1997) 4239–4249. [25] C.P. Chiu, J.Y. Shich, W.R. Chen, Transient natural convection of micropolar fluids in concentric spherical annuli, Acta Mech. 132 (1999) 75–92. [26] W.R. Chen, Transient natural convection of micropolar fluids between concentric and vertically eccentric spheres, Int. J. Heat Mass Transfer 48 (2005) 1936–1951. [27] O. Aydin, I. Pop, Natural convection from a discrete heater in enclosures filled with a micropolar fluid, Int. J. Eng. Sci. 43 (2005) 1409–1418. [28] N.T. Eldabe, E.F. Elshehawey, M.E. Elsayed Elbarbary, Nasser S. Elgazery, Appl. Math. Comput. 160 (2005) 437–450. [29] H.S. Takhar, R.S. Agarwal, B. Bhargava, S. Jain, Mixed convection flow of a micropolar fluid over a stretching sheet, Heat Mass Transfer 34 (1998) 213– 219. [30] Ching-Yang Cheng, Natural convection boundary layer flow of a micropolar fluid over a vertical permeable cone with variable wall temperature, Int. Commun. Heat Mass Transfer 38 (2011) 429–433. [31] Sarifuddin a, Santabrata Chakravarty b, Prashanta Kumar, Heat transfer to micropolar fluid flowing through an irregular arterial constriction, Int. J. Heat Mass Transfer 56 (2013) 538–551. [32] H.M. Badr, S.C.R. Dennis, Time-dependent viscous flow past an impulsively started rotating and translating circular cylinder, J. Fluid Mech. 158 (1985) 447–488. [33] F.M. Mahfouz, H.M. Badr, Flow structure in the wake of a rotationally oscillating cylinder, ASME J. Fluids Eng. 122 (2000) 290–301. [34] J.D. Chung, C.-J. Kim, H. Yoo, J.S. Lee, Numerical investigation on the bifurcative natural convection in a horizontal concentric annulus, Numer. Heat Transfer Part A 36 (1999) 291–307.