Maejo Int. J. Sci. Technol. 2010, 4(01), 43-52 Maejo International Journal of Science and Technology ISSN 1905-7873 Available online at www.mijst.mju.ac.th Communication Modelling of natural convection along isothermal plates and Yin channels using diffusion velocity method AR Ademola A. Dare* and Moses O. Petinrin R Department of Mechanical Engineering,University of Ibadan, Nigeria IB * Corresponding author e-mail: ademola.dare@mail.ui.edu.ng L Received: 19 January 2009/ Accepted: 5 February 2010 / PublisAhed:N 17 February 2010 Abstract: The diffusion velocity method (DVM), aD version of vortex element method (VEM), was used to model the steady staBte, Alaminar natural convection flows along isothermal vertical plates and in isothermal vertical channels. For each case, numerical models were developed using DVM from Ithe vorticity transport equation and the energy equation. This study shows that the difFfusion velocity method is a viable numerical tool at modelling not only fluid flow pro bOlems but also the heat transfer problems. Y Keywords: natural convTection, isothermal plates, isothermal channels, diffusion velocity method I RS Introduction E NManIy V engineering problems are represented by non-linear, partial differential equations. In order to simplify these problems which may be difficult to solve analytically, different numerical meUthods are used to analyse the problems. Some of these numerical methods are: the finite difference method, the finite element method and the Monte Carlo method. Another numerical method that has been used successfully to study and analyse very complex, unsteady fluid flows and thermal engineering problems is the vortex element method (VEM). This method is simple to implement as it requires simple mathematics compared to other numerical methods that may require many mathematical operations such as variational calculus and matrix inversions [1]. An engineering problem involves flows of gases or liquids over solid bodies. Examples include: air flows over cars and aeroplanes; wind blowing over bridges and buildings; and sea waves 44 Maejo Int. J. Sci. Technol. 2010, 4(01), 43-52 slashing against the supporting columns of an off-shore oil rig. Often these flows do not follow the contour of the solid surface completely, but may be formed separate from it, e.g. a wake behind a ship. Such separated flows are difficult to handle by conventional numerical schemes. The main reason for basing the numerical method on the vorticity is that typically only a small portion of the flow contains vorticity. This can lead to significant savings in storage and computational effort [2]. A number of numerical schemes model diffusion by vortex methods without using a mesh, which is based on the Lagrangian approach. However, the VEM is grid-free in implementation [2-3]. Recent applications of the vortex methods based on the Biot-Savart law have been extenYded to numerical prediction of unsteady and complex characteristics of various flows related tRo difficult engineering problems concerning flow-induced vibration, off-design operation of fluAid machinery, automobile aerodynamics, biological fluid dynamics, etc. Kamemoto [4] described Rhow VEM, with its simple algorithms based on physics of flow, has been used to find the vIirBtual operation of fluid machinery (pumps and water turbines), and the calculated flows around bluff bodies, an oscillating airfoil and a swimming fish. L Cheng et al. [5] developed a hybrid vortex method to simulate two-dimensional viscous incompressible flows over a bluff body. It was based on a combinatNion of the diffusion-vortex method and the vortex-in-cell method whereby the flow field is divided Ainto two regions. In the region near the body surface the diffusion-vortex method is used to solvDe the Navier-Stokes equations while the vortex-in-cell method is used in the exterior domaBin. TAhey compared the results obtained with those from the finite difference method and other vor teIx methods and experiments. Which showed that the method is well adapted to calculate two-dimensional external flows at high Reynolds number. Ghoniem and Oppenheim [6] applFied random-walk vortex method to an assortment of problems of diffusion of momentum a nOd energy in one dimension as well as heat conduction in two dimensions in order to assess its vaYlidity and accuracy. The numerical solutions obtained were found to be in good agreement with theT exact solution except for a statistical error introduced by using a finite number of elements. They claIimed further that the error could be reduced by increasing the number of elements or by using enseSmble averaging over a number of solutions. The concept oRf the diffusion velocity method (DVM), a version of the VEM, was extensively discussed by Ogami [7]. The velocity is defined in order that the vorticity is conserved in the transfer of diffusion pVrocEess as it is so in the convection process. Unlike the other vortex element methods, this technique hIandles vorticity equation in a deterministic manner by calculating the diffusion velocity to accountN for diffusion in the flow [8]. Ogami [8] used the DVM to simulate the diffusion of vorticity andU temperature from one-dimensional vorticity transport equation and compared the results with the analytical solutions. The method was successfully used to treat the diffusion equation (Re = 0), the boundary layer and two-dimensional flows around a circular cylinder (Re = 0.1 ~ 107), aerofoil, the Burger equation, and the equations of incompressible fluid. This study employs the diffusion velocity technique to model the natural convection along isothermal plates and in isothermal channels. The channels, consisting of two parallel vertical plates, are asymmetrically heated at uniform wall temperature. The Nusselt numbers are obtained with the velocity and the temperature distributions. 45 Maejo Int. J. Sci. Technol. 2010, 4(01), 43-52 T he Governing Equations The continuity equation and the Navier-Stokes equations for Newtonian, two-dimensional, incompressible flow are presented as follows: u v   0 (1) x y u u u v u 1 P 2 f ( u  2u      x  2  2 ) (2) t x y  x x y Y v u v 2 v v 1 P f ( v  2v      y  2  2 ) (3) Rt x y  y x y A Also, the energy equation, assuming no viscous dissipation or thermal generation, is T u T v T ( 2T  2T R     2  2 ) I B (4) t x y x y Variables u and v are the velocity components in x and y directioLns respectively, fx and fy are the components of the gravitational acceleration in the x and y directio ns respectively, T is the fluid temperature, P is the fluid pressure, t is the time, α is the fluid therNmal diffusivity,  is the kinematic viscosity, and ρ is the fluid density [9-10]. A Following the Lagrangian scheme, the alternative Dexpression of the governing equations of viscous and incompressible flow gives the vorticity tranAsport equation as w w w  2w  2w T  u  v  ( 2  2 )  g IB (5) t x y x y y where vorticity, w v u  F (6) x y Considering the diffusive term oOnly in the vorticity equation, then w 2 ( w  2w  2  2 ) T Y (7) t x y The solution to equation (7) wIas given by Batchelor [11] as r 2 w(r, t)  Re( S4t) (8) 4t where r  (x2 Ey2 ) and Γ is the vortex strength or circulation [12] V NumerNical IFormulations UThe initial velocity,  , is induced by a buoyancy effect created by the temperature difference between the plate and the fluid in contact with its surface, and is given as  g T  ts (9) y where s is the elemental length, t is the time step, β is the volumetric thermal expansion coefficient, and g is acceleration due to gravity. Velocity and temperature vortices with strengths Γ and Θ respectively, which are created on one elemental surface of m elements into which a plate is divided are given by Γq = γΔs (10) 46 Maejo Int. J. Sci. Technol. 2010, 4(01), 43-52 q  TuTqt T  t (11) y y where m L ;   t ; q = 1, 2, 3,...m; and L is the plate length. s The elemental length has a great influence on the strength of the vortices. The vortices are initially distributed and separated at distance Δs before diffusing. The vorticity of each velocity vortex and the temperature of each temperature particle are respectively given by   r 2 Y 1 r 2 2 w(r, t) i e( 4t ) e ( 4t)      (12) 4t   R   (r21 r22 AT (r,t)  i 4t ) ( 4t ) e  e  (1R3) 4t   For flow on a flat plate, each velocity vortex or temperature partiIcBle has a corresponding negative image. The distances between a vortex and the inducing vortex and its corresponding image can be deduced respectively by L r1 (x j , y j )  (x 2 2 j  xi )  ( y j  yi ) and r2 (x j , y j )  (x 2 j Nxi )  (y 2j  yi ) Here i = j = 1, 2, 3,…N, and N is the number of vortices on theA plate surface and in space. However, i is not equal to j most of the times. Figure 1 shows an induciDng vortex with its image and how r1 and r2 are determined from a vortex of interest. A The divergences of equations (12) and I(1B3) are respectively .w and .T , which can be written in the forms: .w w w    F (14) xj y j and T T O .T   Y (15) xj y j The diffusion velocity Tof each velocity vortex is induced by all the velocity vortices and their corresponding images. ASlsoI, the diffusion velocity of temperature particle is induced by all the temperature particlesR and their corresponding images. As discussed by Ogami [8], the diffusion velocities of eachE velocity vortex and temperature particle are now respectively given as Nu wj V[ .w] (16) i1 w NI Nu [ Tj   .T ] (17) U i1 TThe diffusion distances of each velocity vortex and temperature particle can be expressed respectively as d vj  uwj t (18) and dTj  uTj t (19) These distances are added to the original positions to move the vortex and particle to other positions. For convection, each vortex is repositioned to a new location using the induced velocities (u, v) by neighbouring velocity vortices. 47 Maejo Int. J. Sci. Technol. 2010, 4(01), 43-52 AR Y R IB N L A Figure 1. Interaction of vortices onD a vertical flat plate The slip flow condition is met when the veBlociAty on each elemental surface is approximately zero and the numerically calculated temperature isI approximately equal to the initial temperature of the plate. To check for the slip flow condition, thFe v elocity on the elemental surface can be determined as N   v  i i  e    O (20) i1 2r1 2r2  The corresponding temperature onY the surface is determined as N r 21 r 22 T   i e( 4ITt) e( 4t ) e      TB (21) i1 4t   The practical substituRtionS for the given boundary condition is yTB  TwerEfc( j ) (22) 2 t where Tw isI tVhe wall temperature. At each time step, the strength of a vortex is increased by adding  g TNts to equation (9) and Te is substituted for T in equation (11) until the slip flow condition is Uymet. When a vortex is very close to the wall or an inducing vortex, i.e. when r s (r can be r1 or  r  r 2), then i 2 replaces i where it is appropriate. The minimal distance, r0, is equal to s . 2r0 2r   r  Likewise for temperature particles, i replaces i . 2r 20 2r 48 Maejo Int. J. Sci. Technol. 2010, 4(01), 43-52 The velocity and temperature distributions at specific locations in the hydrodynamic and thermal boundary layers are determined when the slip flow condition is meant for the velocities and temperatures on the plate surfaces. The two components of the velocity distribution are obtained as N   y  y   y  y  v (r)   i  j i   imnx    j i  (23) i1 2r1  r1  2r2  r2  N  i  x j  xi  i  x j  x vmny (r)      i     (24) i1 2r2  r2  2r1  r1  Y The temperature at any point can be obtained from R N      T i imn (r)     TB (25) A i1 2r1 2r2  R where m and n are the positions of the velocity or temperature along aInBd normal to the plate respectively. L Natural Convection in Channels N The type of channels described here are two parallel pAlates placed vertically. The number of images of a vortex or a particle is infinite, but for easy computation the number of the images is reduced to eight; therefore, we have one positive vorteAx wDith four positive images and four negative images [1]. As discussed by Petinrin [13], the vortic ityI Bof each velocity vortex and the temperature of each temperature particle can then be respectively given as w(r, t)   2 2 2 2 2 (r1 ) (r2 )F(r3 ) r4 r5   i e 4t  e O4t  e 4t  e ( 4t) e( 4t ) 4t  2 2 2 (r6 ) Y(r7 ) (r8 ) (r 29  e 4t  e 4t  e 4t e 4t )T   (26)   2 2 2 2 2T (r, t) i Se( rI1 r2 r3 r4 r5  4t ) e( 4t ) e( 4t ) ( ) (   e 4t  e 4t) 4t R 2 2(r6 ) (r7 ) (r 28 ) (r 29 E  e 4t  e 4t  e 4t  e 4t)   (27) AlsoI, V  the distances of each vortex from the inducing vortex and its images can be determined as r1 to r9,N where h is the distance of the plates apart, i.e. 2 2U r1 (x j , y j )  (x j  xi )  ( y j  yi ) , r2 (x j , y j )  (x j  x 2 2i )  (y j  yi ) , r3 (x j , y j )  (x j  xi ) 2  (y j  ( yi  2h)) 2 , r 24 (x j , y j )  (x j  xi )  ( y j  ( yi  4h)) 2 , r5 (x j , y j )  (x j  xi ) 2  ( y j  yi  2h) 2 , r6 (x j , y 2 2 j )  (x j  xi )  ( y j  yi  4h) , r7 (x j , y 2 j )  (x j  xi )  (2h  y 2 j  yi ) , 49 Maejo Int. J. Sci. Technol. 2010, 4(01), 43-52 r8 (x j , y j )  (x 2 j  xi )  (4h  y 2 j  yi ) , r9 (x j , y j )  (x j  xi ) 2  (2h  y  y )2j i (28) Results and Discussion The results of the numerical analysis, which was solved with Visual Basic programming language, are automatically displayed on Microsoft Excel Workbook. The velocity and tempeYrature distributions are displayed on the Sheet 1 and Sheet 2 of the workbook respectively. The input parameters used to simulate natural convection on a single plate lying veRrtically are listed in Table 1. The fluid (air) thermophysical properties of the fliud are taken at fiRlm tAemperature. Table 1. The input parameters for natural convection over a single plaIte B Length of the plates 0.5 m L Fluid (air) temperature 100C Plate wall temperature 300C N Coefficient of thermal expansion 0.00341 A(K-1) Kinematic viscosity of air at 200C 0.00D00157 m2/s Thermal diffusivity of air at 200C 0A.000022 m2/s The logarithmic plot of Nusselt number aIgBainst Rayleigh number is presented in Figure 2 by varying the plate wall temperature, Tw, fromF 40 0C to 1200C while keeping the fluid (air) temperature constant at 100C. The fluid properties are taken at film temperature. It can be deduced that the Nusselt number increases with the Rayleigh nu mOber for a fixed plate length as in Table 2. The slope of the plot is 3.47 while the intercept is -28.22Y on the log scale and 6.08E-29 on the normal scale. Therefore, the relatIioTnship between Nusselt number and Rayleigh number is Nu = (6.08E-29)Ra3.47. Comparing this relationship with Nu = 0.59Ra0.25 as reported by Incropera and Dewitt [10], there is mSuch deviation in the two correlations which may be due to convergence difficulties at theE platRe surface. The correlation coefficient of the plot is 0.9065. V32.5 UN I 2 1.5 1 0.5 0 8.5 8.55 8.6 8.65 8.7 8.75 8.8 8.85 8.9 Log(Ra) Figure 2. Logarithmic plot of Nusselt number (Nu) against Rayleigh number (Ra) for a vertical plate Log(Nu) 50 Maejo Int. J. Sci. Technol. 2010, 4(01), 43-52 Table 2. Distribution of Rayleigh number with Nusselt number for natural convection over a single plate Tw (0C) Ra Nu Log(Ra) Log(Nu) 40 336123348 27.0125 8.5265 1.4316 60 479235019 84.0570 8.6805 1.9246 80 594247932 144.6929 8.7740 2.1604 100 655562931 297.4379 8.8166 2.4734 120 706467195 362.2841 8.8491 2.5590 Y R The Nusselt number increases with the wall temperature as presented in Figure 3A. However, the relationship between them is not linear owing to the Nusselt number being also depRendent on some of the fluid properties, e.g. the thermal conductivity, which keeps changing as the wall temperature changes. IB 140 N L 120 100 80 AD A 60 40 B 20 I 0 F 0 50 100 O150 200 250 300 350 400Nusselt Number, Nu Figure 3. EYffect of wall temperature on Nusselt number T The input parameStersI for the symmetrically heated, isothermal plates listed in Table 3 were used to simulate natuRral convection in channels. The channels are two parallel plates placed vertically. The graph of logarithm of Nusselt number against logarithm of Rayleigh number in Figure 4 is obtained by varyEing the magnitude of the space, S, between the channel plates from S = 0.1m to 0.3m at a step ofI 0V.05m. Input parameters are fed from Table 3. The slope of the plot is 0.51 while the intercepNt is -0.75 on the log scale and 0.18 on the normal scale. UTherefore, the relationship between Nusselt number and Rayleigh number is Nu = 0.18Ra 0.51 (Figure 4). The Figure also shows a deviation in the Nusselt-Rayleigh relationship from other workers. Ofi and Hetherington [14] with finite element method obtained a relationship Nu = 0.699Ra0.26; Bodoia and Osterle[15] obtained Nu = 0.56Ra0.25; and Ogundare[1] obtained Nu = 0.43Ra0.21. Temperature (degree Celcius) 51 Maejo Int. J. Sci. Technol. 2010, 4(01), 43-52 Table 3. Input parameters for natural convection in vertical channel Length of plates 0.5 m Gap between plates 0.1m Fluid (air) temperature 100C Wall temperature of first plate 600C Wall temperature of second plate 600C Coefficient of thermal expansion 0.00325 Y Kinematic viscosity of air at 350C 0.0000171 m2/s R Thermal diffusivity of air at 350C 0.0000241 m2/s RA B 3.5 LI 3 N 2.5 DA Model2 A Ofi and Heterington1.5 B Bodoia & OsterleI Ogundare1 0.5 OF 0 5 5.5 6 6.5 7 7.5 8Y Log Ra IT Figure 4. LogarithmR ploSt of Nusselt number (Nu) against Rayleigh number (Ra) for vertical channel ConclusionIs V E NModelling of natural convection in isothermal vertical plates and channels has been sucUcessfully carried out with the diffusion velocity method, a version of the vortex element method. However, a large deviation recorded for correlation of Nusselt number and Rayleigh number for both the plate and the channel with existing correlations may stem from convergence difficulties encountered at the plate surface. From the results obtained, it is established that as the wall temperature increases while keeping the mainstream fluid temperature constant, the thermal boundary layer thickness increases. The study has also established that the diffusion velocity method is a viable numerical tool capable of modelling fluid and heat transfer problems. Log Nu 52 Maejo Int. J. Sci. Technol. 2010, 4(01), 43-52 References 1. M. A. Ogundare, “Vortex element modelling of fluid and heat flow in vertical channels and ducts”, PhD. Thesis, 2006, University of Ibadan, Nigeria. 2. S. A. Subramanian, “New mesh-free vortex method”, 1996, Retrieved 20th Feb, 2008 from http://www.eng.fsu.edu/~dommelen/papers/subram/style_a/node78.html 3. L. Liu, J. Feng, J. Fan and K. Cen, “Recent development of vortex method in incompressible viscous bluff body flows”, J. Zhejiang Univ. SCI, 2005, 6A, 283-288. Y 4. K. Kamemoto, “On contribution of advanced vortex element methods toward virtualR reality of unsteady vortical flows in the new generation of CFD”, J. Braz. Soc. Mech. ScRi. &A Eng, 2004, 26, 368-379. 5. M. Cheng, Y. T. Chew and S. C. Luo, “A hybrid vortex method for flows ovBer a bluff body”, Int. J. Numer. Meth. Fluids, 1997, 24, 253-274. I 6. A. F. Ghoniem and A. K. Oppenhim, “Random element metho d Lfor numerical modeling of diffusional processes”, Proceedings of 8th International ConfNerence on Numerical Methods in Fluid Dynamics, 1982, Aachen, Germany, pp. 224-232. 7. Y. Ogami, “Simulation of heat-fluid motion by the vortex mAethod”, JSME Int. J. Ser. B, 2001, 44, 513-519. D 8. Y. Ogami, “Simulation of heat-vortex interacBtion Aby the diffusion velocity method”, ESAIM Proceedings of 3rd International Worksho p Ion Vortex Flows and Related Numerical Methods, 1999, Toulouse, France, Vol. 7, pp 314-324. 9. F. Kreith and M. S. Bohn, “Principles oFf Heat Transfer”, 5th Edn., West Publishing Co., St. Paul, 1993. 10. F. P. Incropera and D. P. DewYitt, “ OFundamentals of Heat and Mass Tranfer”, 5th Edn., John Wiley & Sons, New York, 2005. 11. G. K. Batchelor, “An InItrToduction to Fluid Dynamics”, Cambridge University Press, Cambridge, 1967. 12. R. I. Lewis, “VortexS Element Methods for Fluid Dynamic Analysis of Engineering Systems”, University PressR, Cambridge, 1991. 13. M. O. PVetinrEin, “Diffusion velocity modelling of forced and natural convection in isothermal plates aInd channels”, MSc. Thesis, 2008, University of Ibadan, Nigeria. 14. O.N Ofi and A. J. Hetherington, “Application of the finite element method to natural convection Uheat transfer from the open vertical channel”, Int. J. Heat Mass Transf., 1977, 20, 1195-1204. 15. J. R. Bodoia and J. F. Osterle, “The development of free convection between heated vertical plates”, J. Heat Transf., 1962, 84, 40-44. © 2010 by Maejo University, San Sai, Chiang Mai, 50290 Thailand. Reproduction is permitted for noncommercial purposes.