11-3-2011.pdf Al-Qadisiya Journal For Engineering Sciences Vol. 4 No. 3 Year 2011 NATURAL CONVECTION IN A WAVY POROUS ENCLOSURE HEATED BY AN INTERNAL CIRCULAR CYLINDER Dr. Muneer A Ismael Mechanical Engineering Department-Engineering College – University of Basra Email: muneerismael@yahoo.com P.O. Box 2067, Ashar Post Office, Basra, Iraq ABSTRACT: Natural convection fluid flow and heat transfer of fluid-saturated porous media heated by an internal circular cylinder inside a wavy enclosure is investigated numerically. The 2D enclosure is composed of two isothermal vertical wavy walls and two adiabatic horizontal flat walls. Darcy assumption and Boussinesq approximation were relied on in this steady, incompressible study. The governing equations were solved using Galerkin finite element method implemented in FlexPDE software package. The performance of enclosure was evaluated by three non-dimensional parameters namely, the Darcy-modified Rayleigh number Ram (100-1000), the waviness ratio (0- 0.35), and the position of the inner heated cylinder (0.45-1.05). The results were presented by visualization of the streamline and isothermal contours and by the local and average Nusselt numbers. It was found that the lower the position of the inner cylinder ( =0.45) is the largest the values of Nusselt number while the influence of the wall waviness ratio is found to be very small. KEYWORDS: Natural convection, porous media, enclosure, wavy wall, Darcy assumption. . --. : . . (Darcy assumption)(Boussinesq approximation) . (FlexPDE) .: - Ram (100-1000) ،(0-0.35)(0.45-1.05) . (streamlines)(isotherms) Al-Qadisiya Journal For Engineering Sciences Vol. 4 No. 3 Year 2011 .( =0.45) . NOMENCLATURES a Wavy wall amplitude (m) X,Y Dimensionless coordinates A Aspect ratio Greek D Diameter of inner cylinder (m) Thermal diffusivity of fluid (m2/s) Da Darcy number Thermal expansion of fluid (K-1) g Gravitational acceleration (m/s2) Waviness ratio h Heat transfer coefficient (W/m2K) Kinematic viscosity (m2/s) H Enclosure height (m) dynamic viscosity of fluid (N/m2S) K Permeability of porous media (m2) Stream function (m2/s) k Thermal conductivity (W/m.K) Dimensionless stream function L Length (m) Density of fluid (kg/m3) Nu Nusselt number o Density of fluid at temperature To (kg/m3) P Pressure (Pa) Dimensionless temperature Ram Darcy-modified Rayleigh number Dimensionless position of inner cylinder S Distance along L (m) Subscripts T Temperature (K) av Average U Velocity vector (m/s) c Cold U Velocity in x-direction (m/s) h Hot V Velocity in y-direction (m/s) i Inner W Average width of enclosure (m) L Left x,y 2D Coordinates R Right 1. INTRODUCTION Fluid flow and heat transfer in porous media have received a pronounced attention in the past few decades. The reason behind this attention is the extensive growing applications of this field in engineering such as solar collectors, thermal insulation, grain storage, filtering and draying process, nuclear research (Ingham et al 2002; Nield and Bejan 2006), production of oil and gas from underground reservoir (Clifford et al 2006). The applications extend also to the biological and medical problems (Kulasiri and Verwoerd 2003; Khaled and Vafai 2003). Most of these applications are simulated by a natural convection process in porous enclosures. Different shaped enclosures were investigated to attain a system of high efficient thermal performance. The early reported studies, generally, were focused on the square or rectangular enclosures as in (Groos et al 1986; Manole and Lage 1992; Goyoeau et al 1996; Baytas and Pop 2002). Nevertheless the rectangular geometry is still investigating but in different boundary conditions as in (Barletta and Lazzari 2005) where a square cavity of isothermal walls and an internal concentric circular heating boundary was studied for different cavity inclination angles. The study of (Oztop 2007) concluded that the inclination of a partially cooled rectangular enclosure is the dominated parameter on heat transfer and fluid flow. Voral et al 2008 studied the effect of amplitude of the sinusoidaly varying temperature profile on the bottom wall of rectangular enclosure for different aspect ratios. Another study on square enclosure subjected to localized heating and salting from one side was done by (Zaho 2008) where the double-diffusive convective flow was considered (double-diffusive: diffusion of matter caused by temperature gradients and diffusion of heat caused by concentration gradient). In (Braga and Lemos 2009) a numerical tool was developed to study the natural convection in a composite cavity formed by three distinct regions, clear, porous, and solid region. Al-Qadisiya Journal For Engineering Sciences Vol. 4 No. 3 Year 2011 However, a limited group of non-rectangular porous enclosures was investigated basing on Darcian and non-Darcian assumptions of porous model. The trapezoidal enclosure has been studied by many researchers; in (Marafi and Vafai 2001) a numerical parallel analysis was developed to study the effect of Darcy modified Rayleigh number Ram, Grashof number, and the inclination of the trapezoidal walls. Extensive numerical studies, penalty finite element with bi-quadratic element analysis was performed to study the wall inclination, uniform and non uniform heating of the trapezoidal base and wide range of Rayleigh, Prandtle and Darcy parameters are edited in (Basak et al 2009) and in another work (Basak et al 2009). A right angle trapezoidal enclosure with heated vertical wall and partially cooled inclined wall was studied by (Voral et al 2009). Three aspect ratios and three locations of cooling part with a range of Ram were examined. In addition, the maximum density effect on buoyancy-induced flow and heat transfer was studied in (Voral et al 2010). Right-angle triangular enclosures were studied numerically for various boundary conditions of the walls, various aspect ratios and ranges of Ram (Voral et al 2006). In addition, the effect of adding a thin fin inside the right-angle triangular enclosure was investigated in (Voral et al 2007). The effect of inserting a square body subjected to four different boundary conditions inside the right-angle triangular enclosure was investigated in (Voral et al 2007). Because of its complexity, fewer published works about wavy wall(s) enclosures filled with porous media were found. The studies reported in (Mahmud and Andrew 2004, Adjlout et al 2002, Mahmud et al 2002, Mahmud et al 2003) deal with wavy wall enclosures but limited to clear fluid only. The works of (Kumar and Shalini 2003, Misirlioglu et al 2005, Khanafer et al 2009, Sultana and Hyder 2007) are concerned with wavy wall(s) porous enclosures based on non- Darcian assumption. In all these studies, besides to the classical parameters, Ra, Gr, etc., the effect of phase wavelength and amplitude of the wavy wall(s) were examined. A wavy porous enclosure heated by an internal body has not been seen previously investigated. Therefore, the present study try to contribute in this field of investigation by studying the natural convection in a wavy wall enclosure filled with fluid-saturated porous medium heated by inside circular cylinder basing on Darcian assumption. 2. MATHEMATICAL MODELING 2.1 Physical model Generally, the low velocity flow which is the case in porous media flow obeys Darcy's law (Clifford et al 2006), gU P K (1) Where K is the permeability of the porous media, P is the pressure, U=(U,V) is the velocity vector and g=(0,-g) is the gravitational acceleration. The effects of inertia and boundary permeability are ignored in this study i.e. basing the modeling according to Darcy assumption which is can be safely applied when Darcy number Da < 10-4 (Clifford 2006). The properties of fluid are considered constant everywhere except the density in the buoyancy term of momentum equation which is Boussineq approximation according to the following linear equation of state: oo TT1 (2) Hence the equations governed the problem are: continuity, Darcy-momentum, and energy equations. The Darcy momentum is obtained by differentiating the x-component of equ.(1) with respect to y and that of y- component with respect to x and subtracting the resulting two equations. The following dimensional set of equations is obtained. Al-Qadisiya Journal For Engineering Sciences Vol. 4 No. 3 Year 2011 0 y V x U (3) x TKg x V y U (4) 2 2 2 2 y T x T y T V x T U (5) y U , x V and scaling the different variables as follow: W x X , W y Y , , ch c TT TT where W is the average width of the enclosure shown in Figure 1. Thus, the governing equs.(3-5) can be written in the following non-dimensional form: X Ra YX m2 2 2 2 (6) 2 2 2 2 YXYXXY (7) Where Ram is the Darcy modified Rayleigh number: WTTKg W KWTTg DaRaRa chchm 2 3 . The local Nusselt numbers for both wavy surfaces and the inside circular surface are calculated from: k Lh Nu LLL , k Lh Nu RRR , k Lh Nu iii (8) Where h is the local heat transfer coefficient given by: n T T k h (9) k is the thermal conductivity, LL and LR are the left and right wavy lengths respectively, Li is the n is a vector normal to the wall. Hence, the local Nusselt numbers are as follows: nL Nu , nR Nu , ni Nu (10) The average Nusselt number over the surfaces can be written as L SL L,av dSL 1 Nu L n , R SR R,av dSL 1 Nu R n , i Si i,av dSD 1 Nu n (11) Al-Qadisiya Journal For Engineering Sciences Vol. 4 No. 3 Year 2011 Due to symmetry about X=0, and steady-state assumption, the average Nusselt numbers along both the wavy walls and the inner cylinder are related according to the overall heat balance; DNuLNuLNu i,avRR,avLL,av (12) 2.2 Boundary conditions The enclosure is made of two horizontal flat walls and two vertical wavy walls with inside circular cylinder (Figure 1). The space between the inside cylinder and the enclosure walls is filled with fluid-saturated porous medium. The flat horizontal walls are kept adiabatic while the wavy walls and the inner cylinder are isothermal but kept at different temperatures (the temperature of inner cylinder is higher). The position of the inside cylinder can be moved along a vertical line X=0. The definitions of the boundaries and their conditions are as follow: 0,0) Y i at Y=0 and at Y=A, with X ranges: 2 1 2 1 X for both Y's 0,0)ii on the left wavy wall A Y XAY 2 2 sin1 2 1 ,0 and on the right wall A Y XAY 2 2 sin1 2 1 ,0 0,1)iii on the inside circular cylinder circumference. waviness =a/W, and A is the aspect ratio = H/W =1.5= 4D, W is the average width of the enclosure. 3. NUMERICAL SOLUTION AND VALIDATION 3.1 Software overview The non-dimensional equations system (equs.6-7) is to be solved numerically by mean of the software package FlexPDE Professional Version 5.0.20 3D. The FlexPDE is a scripted finite element model builder and numerical solver written by user (Gunnar 2005). It performs the operations necessary to turn a description of a partial differential equations system into a finite element model, solve the system, and present graphical and tabular output of the results. It has no pre-defined problem domain or equation list, i.e. the domain geometry and the partial differential equations are totally achieved by the user. The FlexPDE is combined of several modules to provide a complete problem solving system, these are: A script editing module with syntax highlighting provides a full text editing facility and a graphical domain preview. A symbolic equation analyzer expands defined parameters and relations, performs spatial differentiation, and symbolically applies integration by parts to reduce second order terms to create symbolic Galerkin equations. It then symbolically differentiates these equations to form the Jacobian coupling matrix. A mesh generation module constructs a triangular or tetrahedral finite element mesh over a two or three-dimensional problem domain. In two dimensions, an arbitrary domain is filled with an unstructured triangular mesh. In three-dimensional problems, an arbitrary two- dimensional domain is extruded into the third dimension and cut by arbitrary dividing surfaces. The resulting three-dimensional figure is filled with an unstructured tetrahedral mesh A Finite Element numerical analysis module selects an appropriate solution scheme for steady-state, time-dependent or eigenvalue problems, with separate procedures for linear and nonlinear systems. The finite element basis may be either quadratic or cubic. An adaptive mesh refinement procedure measures the adequacy of the mesh and refines the mesh wherever the error is large. The system iterates the mesh refinement and solution until a user-defined error tolerance is achieved. Al-Qadisiya Journal For Engineering Sciences Vol. 4 No. 3 Year 2011 A dynamic time step control procedure measures the curvature of the solution in time and adapts the time integration step to maintain accuracy. A graphical output module accepts arbitrary algebraic functions of the solution and plots contour, surface, vector or elevation plots. A data export module can write text reports in many formats, including simple tables, full finite element mesh data, CDF or TecPlot compatible files. An example of customized graphical outputs is shown in Figure 2 and the edited script of the present work is presented in appendix 1. The results are to be represented by stream lines, isotherms and local and average Nusselt numbers. The contours plots are directly captured from the output screen but the data of Nusselt number were exported to another Grapher tool to be able of showing more than one curve in single plot. 3.2 Software validation To check the validity of the present used software, and enhance the reliability of the obtained results, the average Nusselt number along the perimeter of the inner cylinder Nuav,i was evaluated at five grid densities and for three values of Darcy modified Rayleigh number Ram. The grid density in this package is controlled by imposing a relative error together with a mesh refinement times. The relative error was varied from 10-2 to 10-6 and the mesh refinement times were left to be free. The results are presented in Figure 3. It can be seen that at low Ram, the results are slightly affected by the grid size. But at higher Ram (Ram>500), there is a noticeable sensitivity of the results to the grid size especially at Ram=700. However, this figure indicates that when the relative error 10-4 i.e. higher grid density, the results of Nuav,i are approximately fixed for all Ram's. Accordingly, a relative error of 10-4 was chosen in this study as a compromise between the accuracy of the results and the time consumed in each run. Some cases take a run time of about 900 seconds. The girded domain of 10-4 relative error is shown in Figure 4 which indicates that FlexPDE is adaptively refining the mesh wherever it detects that there are strong curvature in the solution domain. 4. RESULTS AND DISCUSSION To enhance the validation of the present numerical results, a comparison with other works was conducted for two different physical domains. These two domains have flat walls ( =0) with (A=1) i.e. square enclosure, and filled with fluid-saturated porous media according to Darcian assumption. The first case is that of (Barletta and Lazzari 2005) where there exist an internal concentric circular cylinder subjected to a uniform heat flux while the flat walls are kept isothermal at =0. The results are presented in Table1 by average Nusselt number on the internal circumference. The second case which was investigated in more than one published works is free of any internal body and adiabatic horizontal walls and isothermal (but with different temperature) vertical walls (Misirlioglu et al 2005). The results are presented in Table2 by average Nusselt number on the isothermal walls. It seen that there is a very good agreement between the present results of Nuav,i and the two different cases. Therefore, the confidence in using the present numerical package is highly enhanced. The effects of Darcy modified Rayleigh number Ram isotherms, and local and average Nusselt numbers. The examined ranges of Ram is from 100 to 1000; : 0, 0.1, 0.25, 0.35; and : 0.45, 0.55, 0.75 (enclosure center), 0.95, and 1.05. The results are figured and gathered in the following two categories. 4.1 Stream and isothermal lines Figure 5 shows the streamlines inside the enclosure of =0.1 and Ram=100 (upper row) and Ram=1000 (lower row) for three different values of (1.05, 0.75 and 0.45). In all these cases, a double symmetric circulation cells with different directions is noticed. The clockwise direction is denoted by negative sign of stream function while the anticlockwise one is denoted by the positive sign. However, the double circulation cells are formed as the hot fluid moves up from the Al-Qadisiya Journal For Engineering Sciences Vol. 4 No. 3 Year 2011 surroundings of inner cylinder towards the adiabatic upper flat wall and turn to the neighboring wavy wall (in both sides) where it mixes with cold fluid stream falling down along the wavy (cold) walls. Thus, the center of these circulation cells is over the inner cylinder center in all these subfigures. The offset between the two centers increases wi resulting in a more strength circulation occupying large space of enclosure as in the case of =0.45 where the fluid circulation covers most of the enclosure space. When Ram is increased, the strength of stream function is increased also where its value at Ram=1000 is about five times that when Ram=100. This is due the domination of the convection over conduction which in turn leads to a good circulation inside the enclosure Figure 6 shows the isotherms for the previous case ( =0.1, Ram=100, 1000, =1.05, 0.75, 0.45). At low Ram(=100), the upper row of Figure 6, a parallel lines are observed around the inner cylinder except at its upper portion where a plume like distribution is observed and this distribution is grow when the inner cylinder is lowered. It is seen that the temperature gradient between the inner cylinder surface and the porous media is minimum at the source regime of the plume like distribution. For higher Ram and especially at Ram=1000, as presented in the lower row of Figure 6, the isothermal lines take the shape of tangent function or what is called climbing sine function (± sin x ± x). The plume like distribution at the upper portion of the inner cylinder is sharper here and its regime is limited. This complex distribution of isotherms is due to the enhanced convection. It is seen that lowering the position of the inner cylinder also enhance the convection heat transfer due to the increased contribution space of the enclosure filled with porous medium. Generally, a thermal boundary layer appears on the upper part of the wavy walls. This thermal boundary layer becomes thinner at high Ram as it appears from Figure 6. The parameters of Figure 7 are same as that of Figures 5 and 6 except =0.35. There is no specified difference between the distribution and the magnitude of stream function except in the case of Ram=100 and =1.05 where the center of the double circulation cells is below that of the inner cylinder. This because the geometry of the enclosure at =0.35, where the flat adiabatic wall length is significantly small and the upper (and lower) part of the wavy cold walls is semi- horizontal which concentrates the faller cold stream far away from the inner cylinder. On the other hand, the isotherms of Figure 8 which are refer to the same previous case of =0.35 could be classified into two patterns. In the first pattern Ram=100 (upper row of Figure 8), the isotherm behavior is similar to that when =0.1 but the temperature is less at the upper part of the enclosure. While the isotherms of Ram=1000 (lower row of Figure 8) behaves exactly as same as that of =0.1. 4.2 Variations of Nusselt number The local Nusselt number is examined along the circumference of the inner cylinder and the right wavy wall. Figure 9 illustrates the distribution of the local Nusselt number along the circumference of the inner cylinder for three Ram (100, 500, 1000) and three values of (1.05, 0.75, 0.45) and for =0.25. For all these six parameter combinations, the distribution of Nusselt number is symmetric (with one peak) about the vertical line X=0 which appears in figures at Si=0.5. The values of Nui increase with increasing Ram due to the added heat transfer which is seen clearly from the stronger circulation at high Ram shown in Figures 5 and 7. The peaks appear at Si=0.5 i.e. at the lowest point of the circumference because of the temperature gradient here is maximum as it previously clarified on the isothermal lines of Figures 6 and 8. It is seen also that the curvature of the two sides of the bell shaped distribution of Nui diminishes with lowering the inner cylinder position. It becomes completely straight line at =0.45, thus resulting in an increase in the overall heat transfer (the area under curves). The virtue of this increase in heat transfer refers to that at cylinder. Finally the peak of Nui distribution of the case =0.45 and Ram=100 becomes more flat along the segment 0.25