Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 51 - NUMERICAL SIMULATION OF NATURAL CONVECTION IN A LAMINAR TWO – DIMENSIONAL FLOW THROUGH A VERTICAL RECTANGULAR DUCT. Abstract This research uesd to compute the flow field variables of natural convection air flow based on finite difference computational fluid dynamic methods. The problem considered deals with a two- dimensional internal, laminer , isothermal flow over a vertical rectangular duct. In this work, governing equations were solved using finite different forward explicit technique. The flow characterstics are evaluated for Prandtl number at 0.7 and dimensionless inlet velocity equals to 0.005. The results explain that the heat transfer rate is a strong function of flow velocity also the hydrodynamic properties such as pressure, temperature and velocity are increased when the flow moves upward the duct.The results showed a good agreement with other published results. Key word: Natural convection , CFD , Vertical duct , Buoyancy-induced flow. المحاكاة العددیة للحمل الحر لجریان طباقي ثنائي البعد خالل نفق عمودي مستطیل المقطع :الخالصة التي تم دراستها تتعلق المسألةفي هذا البحث تم حساب متغيرات جريان الهواء بالحمل الحر اعتمادا على طرق الفروق المحددة في هذا البحث المعادالت الحاكمة تـم . رارة فوق مجرى عمودي مستطيل الشكلبجريان ثنائي البعد داخلي طباقي متساوي الح وسـرعة 0.7 خواص الجريان تم حسابها لرقم براندل يـساوي (Space Marching).حلها بواسطة تقنية المتغير المكاني دالة لـسرعة الجريـان وأيـضا النتائج أوضحت بأن معدل انتقال الحرارة هو . 0.005الجريان الالبعدية عند المدخل تساوي النتـائج . كالضغط ودرجة الحرارة والسرعة تزداد عندما الجريان يتحرك باتجـاه أعلـى المجـرى الخواص الهيدروديناميكية .أوضحت تطابق جيد مع النتائج األخرى المنشورة Dr. Ahmed Kadhim Hussein College Of Engineering Babylon University Mechanical Engineering Dep. Dr. Hayder Shakir Abdulla College Of Engineering Al-Qadisiya University Mechanical Engineering Dep. Dr.AbbasA.S. AL-Jeebori College Of Engineering Al-Qadisiya University Mechanical Engineering Dep. عباس علیوي الجبوري. د القادسیة جامعة –كلیة الھندسة قسم الھندسة المیكانیكیة حیدر شاكر عبد اهللا. د القادسیة جامعة –دسة كلیة الھن قسم الھندسة المیكانیكیة أحمد كاظم حسین. د جامعة بابل–كلیة الھندسة قسم الھندسة المیكانیكیة Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 52 - List of Symbols: Symbol Description Dimension cP Specific heat at constant pressure. Kj /kg.K Gr Grashof number. Grw Grashof number based on duct width. g Acceleration due to gravity. m/s2 e Specific internal energy per unit mass. J/kg H Dimensionless duct height. h Duct height. m k Thermal conductivity . W/m.K P Dimensionless pressure. p pressure. N/m2 Pr Prandtl number. Q Dimensionless heat transfer rate. q Heat transfer rate per unit area. W/ m2 T Temperature. C U Dimensionless velocity component in z-direction. u Velocity component in z-direction. m/s Uav Dimensionless average velocity at any section of the duct. uav Average velocity at any section of the duct. m/s V Dimensionless velocity component in y-direction v Velocity component in y-direction. m/s w Duct width m y Coordinate in horizontal direction. m z Coordinate in vertical direction. m Greek Symbols β Coefficient of thermal expansion. K-1 θ Dimensionless temperature. µ Dynamic viscosity. N.sec/ m2 υ Kinematic viscosity. m2/s ρ Density. Kg/m3 ZY ∆∆ , Spatial steps in computational domain. Subscript i, j Node symbols indicates position in z and y directions respectively. ∞ Conditions at free stream. w Conditions at surface. av average. Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 53 - Introduction A natural convection flow field is a flow driven by the presence of a temperature gradient.As a result of the temperature difference, the density field is not uniform also.Buoyancy forces will induce a flow current due to gravitational field and the variation in the density field. In general , a natural convection heat transfer is usually much smaller compared to a forced convection heat transfer .It is therefore important only when there is no external flow exists . For example in fluid flow through ducts , solar energy collectors, enclosures and spheres Martynenko and Khramtsov (2005).In the present study ,natural convection due to the laminer flow through a rectangular vertical duct was investigated . Such a case occurs in a number of situations involving the cooling of electrical and electronic equipments, also in the flow through a certain types of fin arrangement.During the past, the experimental and analytical methods were used to simulate the free convection over a limited number of shapes.Ortega and Moffat (1986) studied natural convection from an array of cubical elements in a channel flow.They concluded that the superposition technique could be used to predict temperatures in a free convection channel flow.Pu etal. (1999) presented an experimental results of mixed convection heat transfer in a vertical packed channel with a symmetrical heating of opposite wall.A correlation equation for the Nusselt number in terms of Peclet number and Rayleigh number was obtained from the experimental data.Cadafalch etal. (2002) used a finite volume numerical computation in order to obtain a correlation for the heat transfer in a large air channels considering radiative heat transfer between the plates with different inclination angles.Guimaraes and Menon (2004) studied free and forced convection in an inclined rectangular channel. The system of governing equation was solved using the finite element method . They observed that the inclination angle has a stronger effect on the flow and heat transfer for low Reynolds number.Bain and Armfield (2004) applied a collocation scheme to vertical natural convection problem.Marginal stability curves and critical heat fluxes are obtained for Prandtl numbers from 0 to 1000. Comunelo and Guths (2005) examined the heat transfer coefficient of an isothermal vertical plate.They used novel technology to measure heat flux.The results expected that an increasing of heat transfer coefficient was very useful in heat exchange devices.Abid etal.(2006) studied the free and forced convection in a horizontal vertical duct of transversal aspect ratio close to 2.0 uniformally heated from below and thermally insulated else where.They noticed that the heat flux supplied to the wall induces a secondary flow which manifests itself through natural convection rolls.In this work, a numerical simulation gives the results with a short-time and an accurate computation and the computer program may be changed easily to deal with any other complex shapes such as enclosures, sphere and circular ducts.In the numerical simulation, the governing differential equations are overcame by replacing it with differences, calculated from a finite number of values associated with the computational nodes, which are distributed on a suitable grid over the solution domain. In the present study, the explicit space marching finite difference method was adopted to predict the free convection characteristics of two-dimensional internal laminer isothermal flow over a vertical rectangular duct to compute the primitive variables such as the velocity, pressure and temperature at each grid.In the next section, the mathematical analysis will be given in detail and the style, which used to solve the governing equations was described. Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 54 - Mathematical Analysis The studied case is a natural convection flow of air which passes through a vertical rectangular duct with height ( h ) and width ( w ) and the problem considered is described schematically as shown in Figures ( 1 and 2 ). The inlet and exit sections of the duct are open to ambient. The flow creats since the walls of the duct are at different temperature from that of the fluid which surrounds the duct , so it is described by continuity,momentum and energy equations. These equations are written in a dimensionless form by dividing all dependent and independent variables by suitable constant terms. The solution is obtained using explicit space marching finite difference technique and the following assumptions are considered:- 1. The velocity at inlet section of the duct is considered uniform and equals to the flow average velocity uav. 2. The flow enters the duct through a smooth inlet so the viscous forces are neglected. 3. The flow is considered steady , two-dimensional , moving upward through the duct and symmetrical about the duct center line ( solution from y = 0 to y = w / 2 ). 4. No heat transfer occurs to the entering fluid and the fluid is at atmospheric temperature and the duct walls are heated to the same uniform temperature also all the viscous dissipation and radiative transport are neglected. 5. The duct width ( w ) is considered very small compared with its height ( h ). 6.The flow is considered laminer, since the duct is short and the driving temperature is not too high. 7. Density variations are related only to the buoyancy terms of the momentum equations ( Boussinesq approximation) and all the thermophysical properties are assumed to be constant. 8.Air flow is entering the vertical duct from below ,such that only the buoyancy force is aiding the flow. The governing equations for a two-dimensional, steady, laminer, internal flow expressed in conservation form are:- Continuity Equation: 0 yz = ∂ ∂ + ∂ ∂ vu (1) The Conservation of Momentum Equation in Z-Direction : g y u dz dpu v u u ρµρρ − ∂ ∂ +−= ∂ ∂ + ∂ ∂ )( yz 2 2 (2) Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 55 - The conservation of energy equation is : ))(( 2 2 y T c k y T v z T u p ∂ ∂ = ∂ ∂ + ∂ ∂ ρ (3) The average heat transfer equation rate Oosthuizen and Naylor ( 1999) is dyTTuc z q w pa v )( 1 2/ 0 ∞−= ∫ ρ (4) The boundary conditions used in the present study can be arranged as follows:- 1.Since the flow is symmetrical about the center line,so at y = 0.0 0.0 y = ∂ ∂u and v = 0.0 , 0.0 y = ∂ ∂T (5) 2.Near the wall of the duct ( w / 2) at y = w / 2 then u = 0.0 v = 0.0 and T= T w (6) 3.At inlet plane of the duct ( z=0.0) at z = 0.0 then u = u av ∞= TT and 0.2 2 a vuPP ρ −=− ∞ (7) 4.At exit plane of the duct ( z = h ) at z = h then ∞= PP (8) It is suitable to put these equations in a dimensionless form before applying a numerical scheme to it. The governing equations of natural convection and the boundary conditions can be non- dimensionalized by dividing all dependent and independent variables by suitable constant quantities. The following dimensionless variables are used in the present study:- ∞ ∞ − − = TT TT w θ , Y = y / w and )(* )( 4 2 h w G h wTTg G wr w r = − = ∞ υ β ( 9) ))(( rhG wuw U υ = , V = υ wv * and 2222 4*)( υρ rGh wpp P ∞ − = (10) Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 56 - so, the governing equations (Eqs.(1-4)) and the boundary conditions (Eqs.(5-8)) becomes in a dimensionless form: 0 YZ = ∂ ∂ + ∂ ∂ VU (11) θ+ ∂ ∂ +−= ∂ ∂ + ∂ ∂ )( YZ 2 2 Y U dZ dPU V U U (12) )( 1 YZ 2 2 YP VU r ∂ ∂ += ∂ ∂ + ∂ ∂ θθθ (13) ∫= 2/1 0 1 dYU Z Q a v θ (14) at Y = 0.0 0.0 Y = ∂ ∂U and V = 0.0 , 0.0 Y = ∂ ∂θ (15) at Y = 1 / 2 then U = 0.0 V = 0.0 and 1=θ (16) at Z = 0.0 then U = U av 0.0=θ and 2 2 a vU P −= 0.0=θ (17) at z = H then P=0.0 (18) The dimensionless average velocity (a vU) and the dimensionless height ( H ) of the duct are given by: - )(*)( r a v a v hG wwu U υ = and (19) rhG h H = (20) Numerical Scheme: Explicit space-marching finite difference method is used to solve the set of governing equations. The numerical scheme begins starting from known conditions on the inlet enterance, and marches forward in the z-direction from grid line to another and so on. This method is very effective finite difference technique for laminer and turbulent flow, specially for steady flow. By using this technique a computer code is developed to predict the velocities, pressure and temperature .The continuity , momentum and energy becomes in a finite difference form as follows: Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 57 - 1. Continuity Equation ][ 2 1,11,,1,1,, −−−−− −+− ∆ ∆ −=− JIJIJIJIjiji UUUUZ Y VV ( 21 ) 2. Momentum Equation ji ji ji jijijijiji ji ji ji i ji i JIJI U Z U Z Y UUU Z Y UU U V U P U P UU ,1 ,1 ,1 2 ,11,11,11,11,1 ,1 ,1 ,1 1 ,1 ,1, * 2 *) 2 (* − − − −−−+−−−+− − − − − − − ∆ + ∆ ∆ −+ +∆ ∆ − −+−= θ (22 ) 3. Energy Equation )] 2 (*) Pr 2 [( 1,11,1,12 ,11,11,1 ,1 ,1, Y V YU Z jiji ji jijiji ji jiji ∆ − − ∆ −+∆ += −−+−− −−−+− − − θθθθθ θθ ( 23 ) 4. Average Heat Transfer rate It is computed from the following balancing :- Total heat transfer from intet = Rate of enthalpy crosses duct section – Rate of enthalpy enter the duct, numerically it becomes:- YUUU U Z Q niniiiii ii i av ∆++++= −− *]......0.2 [ 1 1,1,3,3,2,2, 1,1, θθθ θ ( 24 ) Since the fluid velocities associated with free convection computation are typically less than 1.0 m/sec Gubaidullin ( 2002).The calculations are performed at dimensionless inlet velocity = 0.005 and Prandtl number = 0.7.The use of an iterative solution method needs the definition of a convergence and stopping criteria to terminate the iteration process.To make the present solution is stable and to obtain a reasonable accuracy ,so that the solution is not diverage the following relation is used for stability Oosthuizen and Naylor ( 1999): Z∆ < SC 2Y∆ Ui-1,n-1 ( 25 ) where SC is a constant with the range from ( 0 ) to ( 0.5 ) Results and Discussion: Figure (3) shows a dimensionless temperature variations for a two-dimensional laminer flow in a vertical rectangular duct for Prandtl number at 0.7.This figure indicates that the temperature distribution occurs at the region between the line of symmetry and the vertical duct surface. Also, the temperature increases gradually as the flow moves upward through the duct , since buoyancy force continuously increasers through the duct. The higher dimensionless Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 58 - temperature can be noticed at the surface of the duct, fitting to the boundary condition at this region.This result is in good agreement with Sohail (2006) . Figures (4) and ( 5 ) shows a dimensionless velocities variation in longitudinal and lateral directions respectively over a two-dimensional vertical rectangular duct for Prandtl number at 0.7..The figure refers that the velocity values equals to the average velocity at the inlet section of the duct.By comparsion between the two figures, it shows that in natural convection in a vertical rectangular duct , the lateral velocity component is very much smaller than the longitudinal velocity component. Figure (6) shows the dimensionless pressure variation over a two-dimensional laminer flow in a vertical rectangular duct for Prandtl number at 0.7. This figure shows that near the inlet of the duct, the viscous forces effect are high and the pressure drops. Further up the vertical duct the effects of the density changes due to the temperature changes becomes dominant and the dimensionless pressure values increases. The figure explains that the pressure on the inlet section of the duct is below ambient. This result is in good agreement with McBain (1999). Figure (7) shows the dimensionless heat flux variation over a two-dimensional laminer flow in a vertical rectangular duct for Prandtl number at 0.7. This figure shows that the heat flux increases with increasing the height of the duct . This clear increase in the heat flux is due to the increase in the velocity also the other reason of this increasing is due to buoyancy forces which causes an increase in the mass flow rate close to the wall and accelerate the fluid flow so as a result causes an increasing in the amount of the heat flux. Conclusions: The following conclusions can be drawn from the results of the present work :- 1. for capturing the flow field parameters ; a more mesh points are required near the surface of the vertical rectangular duct. 2. Because the flow leaving the duct is parallel to the walls of the duct, the pressure on the exit plane is uniform and equals to the ambient pressure. 3. since the mesh generation has been separated from explicit solver any mesh type can be used. 4. the rise in the value of temperature and the pressure far away from the inlet section of the duct is due to the increasing in the effect of buoyancy forces. 5. the space-marching explicit solution which is used to deal with a two-dimensional vertical rectangular duct proves that it can be used to predict the natural convection flow field parameters accurately and in a short time. 6. from the results obtained, the hydrodynamic properties such as pressure, temperature and velocity are increased when the flow moves upward the duct due to the increasing in the buoyancy force effect. 7. The heat transfer rate is a strong function of velocity. The higher the velocity, the higher the heat transfer rate. 8. No buoyancy force effect occurs in the outer domain of the duct, since the effects of viscosity are negligible. Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 59 - References: Martynenko, O. and Khramtsov,P. “Free convectiove Heat Transfer”, Springer-Verlag, Publiciation, Germany , 2005. Ortega, A. and Moffat ,R.," Experiments on Buoyancy-Induced Convection Heat Transfer From an Array of Cubical Elements on a Vertical Channel Wall " ,Rept. No. HMT-38, 1986,Department of Mechanical Engineering , Stanford University, Stanford,CA. Pu, L. ,Cheng,P and Zhao,T."An Experimental Study of Mixed Convection Heat Transfer in a Vertical Packed Channels ", AIAA Journal of Thermophysics and Heat Transfer , Vo.13 , No.4,1999, pp:: 517-521. Cadafalch,J.,Oliva,G.,Graaf,G.and Albets,X." Natural Convection in a Large Channel with Asymmetric Radiative Coupled Isothermal Plates ", Journal of Heat Transfer , ,2002, pp:: 111-130. Guimaraes, P. and Menon,G."Combined Free and Forced Convection in an Inclined Channel with Discrete Heat Sources ", Mecanica Computacional , Vol.XXIII,2004, pp:: 1-18. Mc Bain, G. and Armfield,S."Linear Stability of Natural Convection on an Evenly Heated Vertical Wall ", 15 th Australasian Fluid Mechanics Conference, University of Sydeny , 2004, Australia. Comunelo, R. and Guths,S. " Natural Convection at Isothermal Vertical Plate: Neighbourhood influence ", 18 th International Congress of Mechanical Engineering, Ouro Preto , MG, 2005. Abid, C. ,Medale,M.,Cerisier,P.and Papini,F. " Mixed Convection in a horizontal Rectangular Duct Heated from Below ", International Journal of Low Carbon Technologies, Vo.1 , No.3,2006, pp:: 236-244. Oosthuizen, P. and Naylor,D. “Introduction to convectiove Heat Transfer Analysis”, McGraw-Hill Company, U.S.A. , 1999. Gubaidullin, A. “Natural Convection Heat Transfer in Two Fluid Stratified Pools with Internal Heat Source”, PhD thesis, Royal Institute of Technology ,Sweden,2002. Sohail, A. "Natural Convection Flow in Parallel Plate Vertical Channels ", MSc. Thesis, Department of Mechanical Engineering , King Fahd University, Dhahran , Saudi Arabia, 2006. McBain , G.D. ," Fully Developed Laminer Buoyant Flow in Vertical Cavities and Ducts of Bounded Section ", Journal of Fluid Mechanics , Vo.401 ,1999, pp:: 365-377. Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 60 - Dimensionless Duct Height (H) D im en si on le ss D uc t W id th ( W ) Figure ( 2 ) Vertical rectangular duct which is studied in the present work Gravity h L eft W all R ight W all Z y w Duct inlet Duct exit Figure(1) Mesh Generation Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 61 - 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless Temperature theta Y-AXIS Figure ( 3 ) Dimensionless temperature variation over a two -dimensional vertical rectangular duct for Prandtl number = 0.7 and dimensionless inlet velocity = 0.005 . 0 0.0005 0.001 0.0015 0.002 0.0025 0.003 0.0035 0.004 0.0045 0.005 0.0055 0.006 0.0065 U U Y-AXIS Z-A XIS Figure ( 4 ) Dimensionless velocity variation over a two -dimensional vertical rectangular duct for Prandtl number = 0.7 and dimensionless inlet velocity = 0.005 . Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 62 - -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 V V Y-AXIS Z-AXIS Figure ( 5 ) Dimensionless velocity variation over a two -dimensional vertical rectangular duct for Prandtl number = 0.7 and dimensionless inlet velocity = 0.005. -1.2E-005 -1.1E-005 -1E-005 -9E-006 -8E-006 -7E-006 -6E-006 -5E-006 -4E-006 -3E-006 -2E-006 -1E-006 PP Y-AXIS Z-AXIS Figure ( 6 ) Dimensionless pressure variation over a two -dimensional vertical rectangular duct for Prandtl number = 0.7 and dimensionless inlet velocity = 0.005 . Al-Qadisiya Journal For Engineering Sciences Vol. 3 No. 1 Year 2010 - 63 - 0 5 E-0 05 0 .0 001 0 .0 001 5 0 .0 002 0 .0 002 5 0 .0 003 0 .0 003 5 0 .0 004 0 .0 004 5 0 .0 005 0 .0 005 5 0 .0 006 0 .0 006 5 0 .0 007 0 .0 007 5 0 .0 008 0 .0 008 5 0 .0 009 0 .0 009 5 Y-AXIS Z-AXIS Q Q Figure ( 7 ) Dimensionless heat flux variation over a two -dimensional vertical rectangular duct for Prandtl number = 0.7 and dimensionless inlet velocity = 0.005 .