Microsoft Word - 3-65 page 41-54.doc IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 41 NATURAL CONVECTION HEAT TRANSFER IN CONCENTRIC HORIZONTAL ANNULI CONTAINING A SATURATED POROUS MEDIA AHMED F. ALFAHAID*1, R.Y. SAKR1 AND M. I. AHMED2 1- Department of Mechanical Technology, Riyadh College of Technology, Saudi Arabia. P.O. Box 91471 Riyadh 11633, Saudi Arabia *alfahaid@gotevot.edu.sa 2- Department of Mechanical Engineering, International Islamic University Malaysia. Abstract: Natural convection in horizontal annular porous media has become a subject receiving increasing attention due to its practical importance in the problem of insulators, such as ducting system in high temperature gas-cooled reactors, heating systems, thermal energy storage systems, under ground cable systems, etc. This paper presents a numerical study for steady state thermal convection in a fully saturated porous media bounded by two horizontal concentric cylinders, the cylinders are impermeable to fluid motion and maintained at different, uniform temperatures. The solution scheme is based on two-dimensional model, which is governed by Darcy-Oberbeck-Boussinesq equations. The finite element method using Galerkin technique is developed and employed to solve the present problem. A numerical simulation is carried out to examine the parametric effects of Rayleigh number and radius ratio on the role played by natural convection heat transfer in the porous annuli. The numerical results obtained from the present model were compared with the available published results and good agreement is observed. The average Nusselt number at the heating surface of the inner cylinder is correlated to Rayleigh number and radius ratio. Keywords: Natural convection, numerical investigation, saturated porous media, finite element method, concentric horizontal annuli. 1. INTRODUCTION Natural convection in horizontal porous annuli has a wide variety of practical applications such as the insulation of aircraft cabin or horizontal pipes, cryogenics, the storage of thermal energy, and the underground cable systems. Lui et al. [1] studied experimentally the natural convection in a horizontal annulus for a fluid layer and showed a multicellular regime for a relatively small radius ratio, (R = 1.15). Later, Bishop and Carley [2] provided photographs showing oscillatory flow regime. Grigull and Hauf [3] used a Mach-Zehnder interferometer and visualized different convective regimes including one where three-dimensional effects were present in the upper part of the layer. Using purturbation method, Mack and Bishop [4] solved the steady two-dimensional equations and although the analysis is valid only for small Rayleigh numbers, the results revealed the IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 42 existence of secondary flows in the upper and lower parts of the layer for very small Prandtle numbers. Kuehn and Goldstein [5] have conducted an experimental study, which included flow visualization and heat transfer measurements. The analogous problem of the thermal convection in eccentric annulus containing viscous fluid (non-porous medium) has attracted some attention in recent literature. Yao [6] developed a purturbation solution for slightly eccentric cylinders, using two-parameter expression in terms of eccentricity and Rayleigh number. For this purpose, Yao [6] used a special coordinates transformation in which the inner circle was transformed into pole. More recently, Prusa and Yao [7] constructed a finite difference method simulation using Yao’s [6] transformed coordinates. Projahn, et al. [8] and by Cho et al. [9] provided numerical simulations of viscous fluid in an eccentric annulus. The former applied non- orthogonal body fitted curvilinear coordinates, while the later used bi-cylindrical coordinate. The experimental work of Fant et al. [10] showed that at fairly high Rayleigh numbers, thermal instability for air appears as steady counter rotating cells near the top of the annulus. Perhaps their most interesting result is that the same flow exhibits hysteresis behavior for small gap widths. As Prandtle number tends to zero, unsteady hydrodynamic instability was demonstrated at high Grashof numbers in the middle of the annulus. Cheddadi et al. [11] solved the same equations using the artificial compressibility method to obtain the pressure field. The tangential velocity component was measured using Laser Doppler anemometry in air-filled annular space. Experimental and numerical values agree well; hysteresis behavior was also reported. Few results for porous layer are present. Caltagirone [12] visualized the thermal field using the Christiansen effect and observed a fluctuating three-dimensional regime in the upper part of the layer even though the lower part remained strictly two-dimensional. Both a perturbation method and finite difference technique were used to solve two-dimensional Boussinesq equations. Burns and Tien [13] examined the variations of the overall heat transfer coefficients with the external heat transfer coefficient and radius ratio by steady- state two-dimensional analyses with the finite difference method and purturbation method. It was indicated that a maximum value of overall heat transfer coefficient existed depending upon the radius ratio. Using finite difference method, Echigo et al. (14) also obtained two-dimensional steady state numerical results taking into account the radiation effect. Bau [15], in recent works, took eccentricity into consideration and demonstrated that heat transfer in the annulus could be optimized by a proper choice of eccentricity. As mentioned above although several works have been done on the problem, most of them were restricted to one flow pattern, the unicellular one, except for the case of narrow annulus. However, the experimental measurements of Caltagirone [12] show, the flow pattern is unlike and not as simple as unicellular one and mainly because of this, the overall heat transfer rates predicted numerically disagree in a margin with those experimental data. Rao el al. [16] investigated analytically with the Galerkin method steady and transient analyses of natural convection in horizontal porous annulus heated from the inner surface. They obtained three families of convergent solutions appearing one after the other with increasing modified Rayleigh number corresponding to different initial conditions. Their predictions were limited to modified Rayleigh number up to 300. They also determined numerically the bifurcation point, which coincide very well with that from the IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 43 experimental observations of Caltagirone [12]. Mota and Saatdjian [17] solved two- dimensional Boussinesq equations using finite difference scheme with an ADI method and successive under relaxation to a very fine grid. They showed that for very small radius ratio and on increasing the Rayleigh number, the steady state regime changes from two to four to six to eight cells without exhibiting a hystresis loop. For radius ratio above 1.7 approximately, closed hystersis loops between ranges containing 2 or 4 cells are obtained. Mota and Saatdjian [18] used an accurate finite difference code to solve two-dimensional Darcy-Boussinesq equations for eccentric, horizontal annulus filled with a saturated porous medium. They showed that reducing the radius ratio or increasing the eccentricity has the same impact on the geometry in the top part of the layer where the convective effects are more pronounced. 2. MATHEMATICAL FORMULATION The problem considered here is a porous layer bounded between two horizontal concentric cylinders of radii Ri and Ro as shown in Fig. (1). The two cylinder walls are assumed to be impermeable. The surfaces of the two cylinders are assumed to be maintained at a constant temperatures Ti and To respectively with Ti To. The governing equations for steady state natural convection with Oberbeck-Boussinesq approximation, Darcy flow are given as follows: 0 u v x y ¶ ¶ + = ¶ ¶ (1) = 0 f p u x K m¶ + ¶ (2a) g = 0 f f p v y K m r ¶ + - ¶ (2b) 2 2 2 2 e f f T T k T T u v x y c x yr æ ö¶ ¶ ¶ ¶ ÷ç+ = + ÷ç ÷÷ç¶ ¶ ¶ ¶è ø (3) r 1 ( )f f rT Tr r bé ù= - -ë û (4) By taking the curl of Eq. (2) and using Eq. (4), we obtain the following equations: 2 r = - g f f K T x ¶ y r b m ¶ Ñ (5) 2 f f e c T T - 0 y x y T k x r ¶y ¶ ¶y ¶ ¶ ¶ ¶ ¶ é ù ê úÑ - = ê úë û (6) where, IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 44 = , = u v y x ¶y ¶y ¶ ¶ - (7) Non-dimensionalizing the variables as defined below; o i i i T -T x y , X , Y = T R RoT q = = - = ; where = , e f f k c y a a r Y The governing equations reduce to the following: 0 = YX XY - 2                 (8) X Ra = *2    (9) where, Ra* = Ra . Da, Ra* is the modified Rayleigh number and is given by: feioiffr f * k/R)TK(T c gRa  Ra; is Rayleigh number and given by: fe 3 ioiffr f k/R)T(T c gRa  and Da; is Darcy number and given by 2 iR/K Da  2.1 Boundary Conditions The problem is assumed to be symmetric about the vertical axis, and as a result, only one half of the flow domain will be considered in this analysis. The boundary conditions are handled as follows:- a- Plane of symmetry: 0 X , 0 =     (10a) b- Inner cylinder surface 1.0 , 0 =  (10b) c- Outer cylinder surface 0 , 0 =  (10c) 2.2 Heat Transfer The local Nusselt number at the inner and outer cylinders surfaces can be calculated from the following equations respectively: IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 45 nk Rh Nu 1YXe ii i 22           (11a) nk Rh Nu RYXe oo o 22           (11b) where n: represent the direction normal to the cylinders surfaces. The steady state average Nusselt numbers at the inner or outer cylinders surfaces are assumed to be equal ( oNu = iNu = Nu) and given by:      0 ii d )(Nu 1 Nu (12a)      0 oo d )(Nu 1 Nu (12b) 3. NUMERICAL SOLUTION The solutions of Eqs. (8) and (9) subject to the boundary conditions specified by Eq. (10) is obtained numerically by using the Galerkin based finite element method [19, 20]. The objective of the finite element is to reduce the system of governing equations into a discretized set of algebric equations. The procedure begins with the division of the continuum region of interest into a number of simply shaped regions called elements. The grid used in the present calculation is illustrated in Fig. (1). The element type which used here is linear triangular element. Fig. 1: Physical geometry and sample grid of the present problem. x y r Ro  IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 46 The approximate expressions of temperature and stream function in an element are given by polynomials in terms of the nodal values and interploation functions. The interpolation functions are derived from the assumption of linear variation of temperature and stream function through the element and are given by the following equation:    3 1m mm e N = (13a)    3 1m mm e N = (13b) where Nm is the usual interpolation function and is defined by: Nm= 1 2A (am + bm X + cmY) (14) where A is the element area and 23321 YXYXa  321 YYb  (15) 231 XXc  The other components are given by cyclic permutation of the subscripts in the order 1,2 and 3. If the approximation given by Eq. (13) is substituted in the governing Eqs.(8-9), and the global errors are minimized using the above interpolation functions Nm as weighting functions, and after performing the weighted inegration over the domain G and the application of Green’s theorem, The present model can be written in the equivalent forms: [K1] {} = {F1} (16a) [K2] {} = {F2} (16b) where; e 1 TTe 1 1 dG ][ [N] ][ [N] dG ) ][ . ][][ . ][ (][               E e XY TT G E e Y N X N Y N Y N X N X N K e         eTT E 1e 1 dY N ]N[ X N ]N[}F{ e               e TT G E 1e 2 dG )Y ]N[ . Y ]N[ X ]N[ . X ]N[ (]K[ e                         E e G eeTT E e ee dG X N Rad Y N N X N NF 1 T 1 2 ][ [N] ][][}{      IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 47 and, E = total number of elements, G= bounded domain,  = domain boundary, Y , X YX       Equations (8) and (9) result in two sets of linear equations which have been solved by Gauss elemination method. The resulting two sets of equations have been solved iteratively through a computer code written here in FORTRAN langauge. The iterative procedure was terminated when the following relative convergence criterion was satisfied:      1i i1i 10-4 (17) where i denote the iteration number performed. 4. MODEL VALIDATION The code was validated by solving the convection problem of two concentric horizontal cylinders for which solutions are available. The obtained results compared with the available published data. Table 1 shows the average Nusselt Number for different previous researchers. A good agreement is found between the present work and the other researchers. Table 1: Comparison of the average Nusselt number for natural convection for convective flow between two concentric cylinders with radius ratio of 2. Caltagirone [12] Rao et al. [16] Bau [21] Facas [22] Facas & Farouk [23] Present Code Grid size 49x49 10x10 30x44 50x50 25x25 10x18 Ra=50 1.328 1.341 1.335 1.342 1.362 1.317 Ra=100 1.829 1.861 1.844 1.835 1.902 1.865 IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 48 5. RESULTS AND DISCUSSION A parametric study of five different radius ratios namely 1.2, 1.4, 1.6, 1.8 and 2 is carried out. In Fig. 2, the variation of the average Nusselt Number (Nu) with the Rayleigh Number (Ra) for a radius ratio R of 1.2 is illustrated. It is noticed from this figure that, for a very small Ra the dominant mode is the pure conduction. As the Ra increases the Nu decreases until Ra reaches 400. Above this value, as Ra increases the Nu increases too. Also, the figure shows a break point at Ra about 800, which may indicate different flow characteristics, as stated by [16-17]. The variation of Nu with Ra for R=1.4 is illustrated in Fig. 3. A slight decrease of Nu is noticed at very low Ra<100, then as Ra increases the Nu increases. A break point at about Ra=140 is also noticed. The variation of the Nu with Ra for R=1.8 is depicted in Fig. 4. It is noticed from the figure that as Ra increases the Nu increases. Also, two breaking points are noticed at Ra=200 and 400 respectively. This may be due to change in flow characteristics at these conditions. Figure 5 illustrates the variation of Nu with Ra for R=2. The variation of Nu with R for different Ra is shown in Fig. 6. It can be concluded that for small Ra (i.e. Ra ≤ 10), where the conduction is the dominant mode of heat transfer, the increase of the insulation thickness (R) leads to decrease of Nu (i.e the decrease of heat loss). As Ra increases (where the natural convection becomes the dominant mode) the useful insulation radius ratio (R) becomes smaller. The heat transfer and flow characteristics in the porous layer for R = 2 and 1.6 for Ra = 120 are illustrated in Figs. 7 and 8 respectively. The natural convection is dominant in the top portion of the figure while the conduction is dominant at the bottom. 500 500 1000 1000 1500 1500 Ra 4.7 4.7 4.8 4.8 4.9 4.9 5 5 5.1 5.1 5.2 5.2 5.3 5.3 5.4 5.4 5.5 5.5 5.6 5.6 5.7 5.7 5.8 5.8 N u 500 500 1000 1000 1500 1500 Ra 2.5 2.5 3 3 3.5 3.5 4 4 4.5 4.5 5 5 N u Fig. 2: Variation of Nu and Ra for R = 1.2. Fig. 3: Variation of Nu and Ra for R = 1.4. IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 49 200 200 400 400 600 600 800 800 Ra 1.5 1.5 2 2 2.5 2.5 3 3 3.5 3.5 4 4 N u 100 100 200 200 300 300 400 400 Ra 1.2 1.2 1.4 1.4 1.6 1.6 1.8 1.8 2 2 2.2 2.2 2.4 2.4 2.6 2.6 2.8 2.8 N u Fig. 4: Variation of Nu and Ra for R = 1.8. Fig. 5: Variation of Nu and Ra for R = 2. 1.2 1.2 1.4 1.4 1.6 1.6 1.8 1.8 2 2 Radius Ratio R 1.5 1.5 2 2 2.5 2.5 3 3 3.5 3.5 4 4 4.5 4.5 N u Ra = 10 Ra = 80 Ra = 160 Ra = 400 Fig. 6: Variation of Nu with the Radius Ratios for different Ra. The heat transfer and flow characteristics in the porous layer for R=1.2 and Ra=120 are illustrated in Fig. 9. Due to the small thickness of the porous layer, the conduction mode is dominant mode of heat transfer in the whole layer. Figure 10 shows the comparison between the predicted and correlated average Nusselt number. IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 50 Fig. 7: Temperature filled () and stream function () contours for R = 2 and Ra =120. T 0.9375 0.875 0.8125 0.75 0.6875 0.625 0.5625 0.5 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 -37.9131 -32. 858 -10.1102 -2.52754 Fig. 8: Temperature filled () and stream function () contours for R = 1.6 and Ra =120. T 0.9375 0.875 0.8125 0.75 0.6875 0.625 0.5625 0.5 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 51 T 0.9375 0.875 0.8125 0.75 0.6875 0.625 0.5625 0.5 0.4375 0.375 0.3125 0.25 0.1875 0.125 0.0625 -2.77668 - 1.85 112 -0.185112 Fig. 9: Temperature filled () and stream function () contours for R = 1.2 and Ra =120. 1 1 2 2 3 3 Ra0.2018 R-1.47 1 1 1.5 1.5 2 2 2.5 2.5 3 3 3.5 3.5 4 4 4.5 4.5 5 5 5.5 5.5 6 6 A ve ra ge N u Nu (Predicted) Nu (Correlated) Fig. 10: The comparison between the predicted and correlated average Nusselt number Nu. 6. CONCLUSIONS The numerical solution of natural convection in a porous medium bounded by two horizontal isothermal cylinders is obtained. The solution of five different radius ratios namely 1.2, 1,4, 1,6, 1.8 and 2.0 are presented. Also, the present code is validated through the comparison of its results with those available in the literature. IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 52 For a small insulation thickness, the dominant heat transfer mode is the conduction regardless of Ra. The useful maximum insulation thickness decrease as Ra increases. From the data obtained, the average Nusselt number is correlated with Rayleigh number and radius ratio as:-Nu = 1.82 Ra0.2018 . R-1.4713 with 6.7 % sum of square errors. REFERENCES: [1] C.Y. Liu, W. K. Muller, and F. Landis, “Natural Convection Heat Transfer in Long Horizontal Cylindrical Anuuli”, Int. Devl. Heat Transfer, Vol. 5, pp. 976-984, 1961. [2] E. H. Bishop, and C.T. Carley, “Photographic Studies of Natural Convection Between Concentric Cylin ders”, Proc. of Heat Transfer & Fluid Mechanics Institute, Stanford Univ., pp. 63-78, 1966. [3] U. Grigull, and W. Hauf, “Natural Convection in Horizontal Annuli,” Proc. 3 rd Int. Heat Transfer Conf., Vol. 2, 1966, pp. 182-195. [4] L.R. Mack, and E.H. Bishop, “Natural Convection between Concentric Cylinders for Low Rayleigh Numbers,” Quart. J. Mech. Appl. Math., Vol. 21, pp. 223-241, 1968. [5] T.H Kuehn, and R.J. Goldstien, “An Experimental and Theoritical Study of Natural Convection in the Annulus between Concentric Cylinders,” J. Fluid Mechanics, Vol. 74, pp. 695-719, 1976. [6] L.S. Yao, “Analysis of Heat Transfer in Slightly Eccentric Annuli,” ASME J. of Heat Transfer, Vol. 102, 1980, pp. 279-284. [7] J. Prusa, and L.S. Yao, “Natural Convection Heat Transfer between Eccentric Cylinders,” ASME J. of Heat Transfer, Vol. 105, 1983, pp. 108-116. [8] U. Projahn, H. Rieger, and H. Bear “Numerical Convection between Concentric and Eccentric Cylinders,” Numerical Heat Transfer, Vol. 4, pp. 131-146. [9] C.H. Cho, K.S. Chung, and K.H. Park, “Numerical Simulation of Natural Convection in Concentric and Eccentric Cylindrical Annuli,” ASME J. of Heat Transfer, Vol. 104, 1982, pp. 624-630. [10] D.B. Fant, A. Rothmayer, and J. Prusa, 1989, “Natural Convective Flow Instability Between Horizontal Concentric Cylinders”, Numerical Methods in Laminar and turbulent Flow, Pineridge Press, Swansea, U.K., Vol. 6, Part 2, pp. 1047-1065. [11] A. Cheddadi, J.P. Caltagirone, and A. Motjtabi, 1989, “An Experimental and Numerical Study of Natural Convection in Horizontal Cylindrical Annuli”, Numerical Methods in Laminar and turbulent Flow, Pineridge Press, Swansea, U.K., Vol. 6, Part 2, pp. 1047-1065 [12] J.P. Caltagirone, 1976, “Thermo-Convective Instability in Porous Medium Bounded By Two Concentric Horizontal Cylinders,” J. of Fluid Mechanics, Vol. 76, pp. 337-362. [13] P.J. Burns, and C.L. Tien, “Natural Convection in Porous Media Bounded by Concentric Spheres and Horizontal Cylinders,” Int. J. of Heat Mass Transfer, Vol. 22, pp. 929-939. [14] R. Echigo, S. Hasegawa, S. Tottori, H. Shimonure, and Y. Okamoto, 1978, “An Analysis on the Radiation and Free Convective Heat Transfer in Horizontal Annulus with Permeable Insulators,” Proc. 6 the International Heat Transfer Conf., Vol. 3, pp. 385-390. [15] H.H. Bau, 1984, “Low Rayleigh Number Thermal Convection in Saturated Porous Medium Bounded by Two Horizontal Eccentric Cylinders,” ASME J. of Heat Transfer, Vol. 106, 1984, pp. 166-175. IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 53 [16] Y.F. Rao, K. Fukuda, and , S. Hasegawa “Steady and Transient Analysis of Natural Convection in a Horizontal Porous Annulus with Galerkin Method,” ASME J. of Heat Transfer, Vol. 109, 1987, pp. 919-927. [17] J.P.B. Mota, and F. Sattdjian, “Natural Convection in Porous Horizontal Cylindrical Annulus,” ASME J. of Heat Transfer, Vol. 116, 1994, pp. 621-626. [18] J.P.B. Mota, and F. Sattdjian, “On the Reduction of Natural Convection Heat Transfer in Horizontal Eccentric Annuli Containing Saturated Porous Media,” Int. J. Num. Methods Heat Fluid Flow, 7 (4), pp. 401-416, 1997. [19] S. S. Rao, The Finite Element Methods in Engineering, Pergamon Press, 1982. [20] D.W. Pepper, and J.C. Heinrich, The Finite Element Method- Basic Concepts and Applications, Hemisphere Publishing Corporation, 1992. [21] H.H. Bau, “Thermal Convection in a Horizontal, Eccentric Annulus Containing a Saturated Porous Medium- An Extended Purturbation Expansion”, Int. J. Heat Mass Transfer, Vol. 27, pp. 2277-2287, 1984. [22] G.N. Facas, “Natural Convection From a Buried Pipe with External Baffles”, Numerical Heat Transfer, Part A, 27:595-609, 1995. [23] G.N Facas, and B. Farouk, “Transient and Steady State Natural Convection in a Porous Medium between Two Concentric Cylinders”, ASME J. Heat Transfer, Vol. 105, pp. 660- 663, 1983. NOMENCLATURE c Specific heat of fluid at constant pressure Da Darcy number; 2iR/K Da  g Gravity acceleration G Bounded domain h Heat transfer coefficient ke Effective thermal conductivity of the porous medium K Permeability Nu Nusselt number n Normal direction R Radius ratio, Ro/Ri Ri Inner cylinder radius Ro Outer cylinder radius Ra Rayleigh number; fe 3 ioiffr f k/R)T(T c gRa  T Temperature u Velocity component in x – direction v Velocity component in y – direction x,y Cartesian coordinates X,Y Dimensionless Cartesian coordinates Greeks  Thermal diffusivity of the porous medium  Thermal expansion coefficient of the fluid  Circumferential angle  Domain boundary IIUM Engineering Journal, Vol. 6, No. 1, 2005 A. F. Alfahaid et al. 54  Dimensionless temperature; (T-To)/(Ti-To)  Dynamic viscosity  Density of fluid  Kinematic viscosity  Stream function  Dimensionless stream function; / Subscripts e effective f fluid i. inner o outer r reference x,y Cartesian components Superscripts e element level i iteration number T transpose  average