IJCPE Vol.11 No. 1 (March 2010) Iraqi Journal of Chemical and Petroleum Engineering Vol.11 No.1 (March 2010) 29-45 ISSN: 1997-4884 Numerical study of the mixed convection flow over a square cylinder Sajida Lafta Ghashim Jassim Mechanical Engineering Department - College of Engineering - University of Baghdad – Iraq Abstract In this work, a numerical study is performed to predict the solution of two – dimensional, steady and laminar mixed convection flow over a square cylinder placed symmetrically in a vertical parallel plate. A finite difference method is employed to solve the governing differential equations, continuity, momentum, and energy equation balances. The solution is obtained for stream function, vorticity and temperature as dependent variables by iterative technique known as successive over relaxation. The flow and temperature patterns are obtained for Reynolds number and Grashof number at (Re= -50,50,100,-100) (positive or negative value refers to aidding or opposing buoyancy , +1 assisting flow, -1 opposing flow) and (102 to 105) , respectively. The results displaced that the recirculation length above the cylinder increases with the increase in Gr number and the average Nu number is the highest at the lower surface of the cylinder, while is the lowest at the top of the cylinder surface. A comparison between the obtained results and the published computational studies has been made and it showed a good agreement. Keywords: CFD, Square cylinder; Mixed convection; Blockage ratio; Nusselt number Introduction The combined forced and free convection heat transfer from cylinders (of circular or square cross-section) has numerous industrial applications such as cooling towers, oil and gas pipelines, tubular and compact heat exchangers, cooling of electronic components, flow dividers in polymer processing applications and so on. The analysis of mixed convection is more complicated than that of the pure forced or free convection alone. Most of the important work on mixed convection past bluff bodies deals with the heat transfer around a heated circular cylinder. An account of many of the experimental results and correlations for forced convection heat transfer from a circular cylinder has been given by Morgan [1]. There are, however, comparatively fewer numerical studies on this topic than for isothermal flow past a circular cylinder. Dennis et al. [2] obtained the numerical solution for the steady laminar forced convection from a circular cylinder. Karniadakis [3] used the spectral element method to obtain forced-convection heat transfer from an isolated circular cylinder in cross flow, for both steady and unsteady periodic flow, for Re≥200. Yang et al. [4] obtained a numerical solution for the transient laminar forced convection from a circular cylinder at Re=100, 200, and 500. Saha [5] studied the effects of the heat transfer and flow analysis of free convective flow past a square cylinder placed symmetrically in a vertical parallel plate channel. The flow is found to be unstable when the Grashof number crosses the critical value of 3.0x104. Farouk and Güceri [6] computed the mixed and free convection from an isothermal circular cylinder in a two- dimensional vertical channel with adiabatic walls in the steady flow regime. Biswas et al. [7] investigated the cross-flow of air over a horizontal square cylinder confined in a channel (isothermal walls) of the blockage ratio of 0.25 for a range of values of the Reynolds and University of Baghdad College of Engineering Iraqi Journal of Chemical and Petroleum Engineering Numerical study of the mixed convection flow over a square cylinder 30 IJCPE Vol.11 No.1 (March 2010) Grashof numbers. They reported the maxima and minima in the velocity and temperature profiles close to the bottom walls of the channel, respectively. They also reported that the periodicity of flow and asymmetry of the wake can occur at lower Reynolds numbers than that in pure forced convection. Turki et al. [8] have numerically investigated the mixed convection from a horizontal square cylinder to air for 120≥Re≥200 for Ri up to 0.1 for a blockage ratio of 0.25. The value of the critical Reynolds number (onset of periodic flow) decreases while the Strouhal number increases with the increasing Richardson number. They also proposed correlations for Nusselt number at different values of the Richardson number (Ri=0, 0.05 and 0.1). Atul Sharma and Eswaran [9] studied the effects of the flow structure and heat transfer characteristics of an isolated square cylinder in cross flow numerically for both steady and unsteady periodic laminar flow in the two-dimensional regime, for Reynolds numbers of 1 to 160 and a Prandtl number of 0.7. The effect of vortex shedding on the isotherm patterns and heat transfer from the cylinder is discussed. Heat transfer correlations between Nusselt number and Reynolds number are presented for uniform heat flux and constant cylinder temperature boundary conditions. The current study as shown in Figure(1) is focused on the effect of aiding and opposing buoyancy over a square cylinder, placed symmetrically in a vertical parallel plates. We focused on the calculation of mean and average Nusselt number along the position on the surfaces of cylinder. Proplem description The problem studied in this paper is a two dimensional laminar buoyancy flow and heat transfer in flow past a square cylinder placed symmetrically in a vertical parallel plate see Figure (1). The considered Grashof number ranged from (10 2 – 10 5 ), the working fluid was air with (Pr=0.71).A square cylinder with side B, the dimensionless length scale, is heated or cooled to a constant temperature Tw. It is exposed to a constant free- stream upward velocity and temperature represented by u∞ and T∞, respectively.Lu the dimensionless distance between the top surface of the cylinder and exit plane, Ld the dimensionless distance between the inlet plane and the bottom surface of the cylinder and L the total dimensionless height of the computational domain. Figure (1) Computational domain for the flow around a square cylinder and the associated parameters. Mathematical model The steady, conservative, dimensionless from of the Navier- Stockes equations in two dimensions for the incompressible laminar flow is given as follows [5]: Continuity: 0      Y V X U (1) X- momentum:                           2 2 2 2 Re 1 Y U X U X P VU Y UU X (2) Y- momentum:      22 2 2 2 ReRe 1 Gr Y V X V Y P VV Y UV X                        (3) Energy:                    2 2 2 2 ReP r 1 YXY V X U  (4) With   u u U ,   u v V , B x X  B y Y  , 2   u p P  ,    TT TT w  Sajida Lafta Ghashim Jassim IJCPE Vol.11 No.1 (March 2010) 31 Vorticity transport and stream function:                       2 2 2 2 2 Re 1 Re )()( YXX Gr V Y U X   (5) The only non- zero component of the vorticity is: Y U X V       (6) From the definition of stream function which verifies the continuity equation, vertical and horizontal components can be written as: X V     (7) Y U     (8) By substituting equation (7) and (8) into equation (6) to obtain the following stream equation:    2 2 2 2 2        YX (9) Boundary conditions Solid surface of the cylinder U=0, V=0, θ=1 Left and right walls (boundary) U=0, V=0, 0   Y  Bottom boundary (inlet) V=0, U=1, θ=0 Top boundary (out let) 0   Y  Numerical Solution The numerical analysis of steady state convection problem may be based on a system of heat convection equation of the form [10]:         2 2 2 2 YX (10) X Gr V YY U XX                              2 ReRe 1 Re 1 (11) 0 RePr 1 RePr 1                           V YY U XX (12) The numerical solution of this system may be obtained by solving its difference system with some available iteration procedure. Consider the second order conservative, monotonic finite difference scheme approximating system of equations (10), (11) and (12). The scheme has been used for steady convective problems involving wide ranges of process parameters and has given good results. It employes the integrointerpolation method. A system of difference equations is obtained by integrating the original system (10), (11) and (12) over a mesh Di ( xi-1/2 ≤ x ≤ xi+1/2 , yi-1/2 ≤ y ≤ yi+1/2 ) shown in figure ( 2 ) [10]. Following is the procedure in [10], the governing finite difference equations for (, , and ) can be written in the standard five point formula form. These finite difference equations which subject to appropriate boundary conditions are solved by an iterative method known as successive substitution. If ( s,s, and s) denote functional values at the end of sth iteration, the value of ( ,, and ) at (s+1)th iteration level are calculated from the following expressions: )()1( 1 1,, 1 ,1,1, 1 ,       s ji s ji s ji s ji s ji s ji dcba A F F      (13)                        1 ,1 1 ,12 1 1,1, 1 ,1,1 1 , 1 , Re 5.0)()1( s ji s ji s ji s ji s ji s ji s ji s ji Gr dcba A F F      (14) )( 4 )1( 1 1, 21 1,1, 1 ,1,1 1 , 1 ,         s ji s ji s ji s ji s ji s ji s ji hdcba F F     (15) Numerical study of the mixed convection flow over a square cylinder 32 IJCPE Vol.11 No.1 (March 2010) Where (s) is the iteration number, and  FF , and F are the relaxation parameters which depend on the mesh size and fluid mechanical parameters. For all cases in present study, an overrelaxation parameters  FF , and F equal to 0.7, 0.5, and 0.7 [10]. A converged solution was defined as one that meets the following criterion for all dependent variables. 4 ),( 11 ),( 10/)(   ji sss jiError  (16) Figure (2) the mesh Di Numerical calculation of the vorticity values on the walls 2 *3 ,1 2 1,, , wijiji wi h        Vorticity equation on vertical wall 2 *3 ,1 2 ,1, 1, jijiji ji h         Vorticity equation on horizontal wall 2 *3 1, 2 1,, 1,       jijiji ji h   Calculation of average Nusselt number The average Nusselt number at the wall is calculated from the temperature gradient at the wall by the following form [10]: k hB Nu  Where: N T kTTh w     )( Then )(      TT N T B Nu w The above equation can be written in dimensionless form as follows: N Nu jL     )( Using the third order polynominal function (θ) , the final form of ( )( jL Nu ) equation is [ 10 ]:  jijijijL h Nu ,2,1,)( 43 2 1    Calculation of Mean Nusselt Number The mean Nusselt number can be calculated from equation below[ 10 ]:  B Nudx B Nu 0 1 The integration can be evaluated by using numerical integration (simpsions rule) as below [10]:           )( 2 )2( 1 )1()1( 24 3 nL j jL j jLL NuNuNuNu h Nu Where: j1= 2,4,6,…….., Y-1 j2= 3,5,7,…….., Y-2 Sajida Lafta Ghashim Jassim IJCPE Vol.11 No.1 (March 2010) 33 Results and Discussion In this work, the results were obtained for heat transfer by laminar mixed convection flow over a square cylinder. The results obtained using the numerical solution, computer program (Fortran 90) and tecplot program. Figures (3-6) show the distribution of the average Nusselt number versus position along the cylinder surface (for different Gr values). Has been shown for the left half of the square cylinder (the other half is symmetric), for various Gr number the figures show that the Nusselt number along the left half of the bottom face (ab) of the cylinder increases .For the left surface (bc), it has a minimum, this is due to turning. In the separation region (cd), there is a local minimum. When the buoyancy force is in the direction of the flow (buoyancy aided), it can be seen that the average Nusselt number increases with the increase of Gr number, an opposite observations are obtained for the case, when the buoyancy force is opposite of the flow direction (buoyancy opposed). Figures (7-10) to (10-13) show the isotherms and streamlines contour for different Gr number values in the vicinity of the cylinder, under the influence of buoyancy. The mechanism of the flow occurs when the fluid near the hot wall is heated causing the density to be decreased and the fluid will be start to move towards the cold wall. It can be seen that the values of isotherms at the cylinder surface increased when Gr number increased. When decreasing (Gr/Re 2 ) < 0 , increases the width of the wake, and the separation point moves from the downstream corners to the upstream end corners, due to the increased negative buoyancy force of the colder air. It also reduces the velocity of the fluid in the wake region, opposite observations are obtained when (Gr/Re 2 ) > 0. So to discern the effect of the recirculation region on the isotherms and heat transfer from the cylinder. The maximum crowding of the isotherms is seen on the bottom face, indicating the highest Nusselt number, as compared to other surfaces of the cylinder, since the thermal boundary layer growth starts from this face. Figure (15) shows the variation of mean Nusselt number versus the Gr number with different values of Re number. This figures show that when the buoyancy force is in the direction of the flow the value of mean Nusselt number increase with increased Re number, but an opposite effect was found for the case when the buoyancy force is opposite of the flow direction (buoyancy opposed). Figure (16) shows average Nusselt number versus position along the cylinder surface at different Gr values (Refernce[5] ) . It has been seen that the Nusselt number is at its highest at the lower surface of the cylinder. On the side walls it is lower because fluid from bottom of the cylinder becomes heated and flows past the cylinder. The Nusselt number is the lowest at the top of the cylinder surface because of the low velocity fluid that forms the vortices (bubbles). The present numerical results are compared with the published results as shown in Figure (17). From the figure, the comparison indicated a good agreement. Conclusions 1. The strength of the re-circulating vortices is increased with the increase of Gr number. 2. The average Nusselt number is the highest at the lower surface of the cylinder, while is the lowest at the top of the cylinder surface because of the low velocity fluid. 3. The maximum crowding of the isotherms is seen on the bottom face, indicating the highest Nusselt number, as compared to other surfaces of the cylinder, since the thermal boundary layer growth starts from this face. Nomenclature Symbol Description Unit B width of the square cylinder p c Specific heat J/Kg.K g Gravitational acceleration m/sec 2 Gr Grashof number ( 2 3 )(      TTBg Gr w ) H Dimensionless width of the computational domain k Thermal conductivity of fluid W/m.K Lu Dimensionless distance between the top surface of the cylinder and exit plane Ld Dimensionless distance between the inlet plane and the bottom surface of the cylinder L Dimensionless height of the computational domain N Cylinder surface normal direction L Nu Average Nessult number Nu Mean Nessult number p Air pressure N/m2 P Dimensionless air pressure (   upP  ) Numerical study of the mixed convection flow over a square cylinder 34 IJCPE Vol.11 No.1 (March 2010) Pr Prandtl number ( kc p /Pr  ) Re Reynolds number (  /Re Bu   ) Ri Rechardson number (Ri= Gr/Re 2 ) T Temperature K  T Free stream temperature K U Dimensionless velocity component in x-direction u∞ Free stream velocity m/s V Dimensionless velocity component in y-direction x Horizontal axis m X Dimensionless horizontal axis )( BxX  y Vertical axis m Y Dimensionless vertical axis )( ByY   Dimensionless temperature (     TT TT w  )   Dynamic viscosity of the fluid Pa.s   Kinematic viscosity of the fluid m 2 /sec   Density of the fluid Kg/m3  Dimensionless stream function  Vorticity References 1. Morgan, V. T., “The Overall Convective Heat Transfer from Smooth Circular Cylinders", Adv. Heat Transfer, Vol. 11, pp. 199–264, (1975). 2. Dennis, S. C. R., Hudson, J. D and Smith, N. " Steady Laminar Forced Convection from a Circular Cylinder at Low Reynolds Numbers", J. Phys. Fluids, Vol. 11, pp. 933–940, (1968). 3. Karniadakis, G. E," Numerical Simulation of Forced Convection Heat Transfer from a Cylinder in Cross.Flow", Int. J. Heat Mass Transfer, Vol. 31, pp. 107–118, (1988). 4. Yang, Y. T., Chen, C. K. and S. R. Wu, "Transient Laminar Forced Convection from a Circular Cylinder Using a Body-Fitted Coordinate System", J. Thermophys., Vol. 6, pp. 184–188, (1992). 5. Saha, A. K. "Unsteady free convection in a vertical channel with a built- in heated square cylinder", J. Numerical Heat Transfer, 38, Part A, P.P 795-818,( 2000). 6. Farouk, B. Güceri, S.I. " Natural and mixed convection heat transfer around a horizontal cylinder within confining walls",J. Num. Heat Transf. Part A, 5,PP. 329–341, (1982) . 7. Biswas, G. , Laschefski, H. , Mitra, N.K. ,and Fiebig, M."Numerical investigation of mixed convection heat transfer in a channel with a built-in square cylinder", J. Num. Heat Transf. Part A, 18 , PP. 173–188 (1990). 8. Turki, S. Abbassi, H. and Nasrallah, S.B., " Two- dimensional laminar fluid flow and heat transfer in a channel with a built-in heated square cylinder", Int. J. Therm. Sci., 42,PP. 1105–1113, (2003). 9. Atul Sharma and V. Eswaran "Heat and fluid flow across a squre cylinder in the two dimensional laminar flow regime".J.Numerical heat transfer , Part A , pp. 247- 269, (2004). 10. Nogotov, E.F. "Applications of Numerical heat transfer", Hemisphere publishing corp. Washington. (1978). Sajida Lafta Ghashim Jassim IJCPE Vol.11 No.1 (March 2010) 35 Fig.(3) Average Nusselt number versus position along the cylinder surface (different Gr values at Re=50) Fig.(5) Average Nusselt number versus position along the cylinder surface (different Gr values at Re=100) Fig.(4) Average Nusselt number versus position along the cylinder surface (different Gr values at Re=-50) Fig.( 6) Average Nusselt number versus position along the cylinder surface (different Gr values at Re=-100) Numerical study of the mixed convection flow over a square cylinder 36 IJCPE Vol.11 No.1 (March 2010) Fig. (7( Isotherm contour maps at Re=50 for different values of Gr number Sajida Lafta Ghashim Jassim IJCPE Vol.11 No.1 (March 2010) 37 Fig.( 8) Isotherm contour maps at Re=-50 for different values of Gr number Numerical study of the mixed convection flow over a square cylinder 38 IJCPE Vol.11 No.1 (March 2010) Fig. (9) Isotherm contour maps at Re=100 for different values of Gr number Sajida Lafta Ghashim Jassim IJCPE Vol.11 No.1 (March 2010) 39 Fig.(10) Isotherm contour maps at Re=-100 for different values of Gr number Numerical study of the mixed convection flow over a square cylinder 40 IJCPE Vol.11 No.1 (March 2010) Fig.) 11) Streamline contour maps at Re=50 for different values of Gr number Sajida Lafta Ghashim Jassim IJCPE Vol.11 No.1 (March 2010) 41 Fig. (12) Streamline contour maps at Re=-50 for different values of Gr number Numerical study of the mixed convection flow over a square cylinder 42 IJCPE Vol.11 No.1 (March 2010) Fig. (13) Streamline contour maps at Re=100 for different values of Gr number Sajida Lafta Ghashim Jassim IJCPE Vol.11 No.1 (March 2010) 43 Fig.(14) Streamline contour maps at Re=-100 for different values of Gr number Numerical study of the mixed convection flow over a square cylinder 44 IJCPE Vol.11 No.1 (March 2010) (a) Re=50 (b) Re=-50 (c) Re=100 (d) Re=-100 Fig. (15)Variation of mean Nu number with Gr number at different values of Re number 0 1 2 3 4 5 6 0 20000 40000 60000 80000 100000 120000 Gr N u 0 1 2 3 4 5 6 7 0 20000 40000 60000 80000 100000 120000 Gr N u 2 2.5 3 3.5 4 4.5 5 0 20000 40000 60000 80000 100000 120000 Gr N u 2 2.5 3 3.5 4 4.5 5 5.5 6 0 20000 40000 60000 80000 100000 120000 Gr N u Sajida Lafta Ghashim Jassim IJCPE Vol.11 No.1 (March 2010) 45 Fig.(16) Average Nusselt number versus position along the cylinder surface (different Gr values ) Refernece [5] Fig.(17) Comparsion of average Nusselt number versus position along the cylinder surface (at Gr=1x10 4 ) with reference [5]. 0 5 10 15 20 25 0 5 10 15 20 25 position along surface cylinder N u x Gr=9x10000 Gr=4x10000 Gr=2.25x10000 Gr=1x10000 Gr=6.4x1000 Gr=1.6x1000 0 2 4 6 8 10 12 14 0 5 10 15 20 position along surface cylinder N u x Present study Refernce [4]