Research Journal of Recent Sciences _________________________________________________ ISSN 2277-2502 Vol. 4(1), 9-16, January (2015) Res.J.Recent Sci. International Science Congress Association 9 Investigation of Sphericity Assumption for Bubbles in the Numerical Simulation of Bubbly flow in a Vertical CylinderSeyed Emad Hosseini¹, Farshad Abbasi² and Ahmad Farmahini Farahani³ Deputy Director of Research, Institute for International Energy Studies (IIES), IRAN Department of Engineering, Islamic Azad University, Marivan Branch, IRAN Director of Research, Institute for International Energy Studies (IIES), TEHRAN, IRAN Available online at: www.isca.in,www.isca.me Received 27th September 2013, revised 13th January 2014, accepted 4th April 2014Abstract To investigate the sensitivity of the numerical results to the sphericity of the bubbles, numerical simulation of bubbly flow caused by the axial gas injection into a vertical cylinder is studied. The mathematical model solves transport equations for the variables of each phase including the interphase mass and momentum exchange. Both spherical and non-spherical bubbles are investigated. Therefore, two methods are used to calculate the drag coefficient. The first is the empirical correlation obtained from the standard drag curve. This method is suitable when the bubbles remain spherical in shape and the surface-tension effects are negligible. Another approach is the "dirty water" model. This model covers the complete range of Reynolds numbers, Weber numbers, and is suitable for the various shapes of bubbles. The predicted results are in good agreement with experimental data available in the literature. It is found that the numerical results computed by the "dirty water" model are closer to those of the experimental data. Keywords: Gas injection; bubbly flow; two-fluid model; non-spherical bubbles; cylinder. Introduction Bubble-driven flows occur widely in industrial applications such as bubble columns, gas absorption, oxidation processes, hydrogen production, chlorine process, electroplating, metal purification, etc. In order to enhance the rate of transport between the phases in many direct contact heat and/or mass exchange processes, a condensable gas is often injected into a liquid through a submerged orifice or nozzle. They are also used commonly in liquid metal containers during steel making and in nuclear cooling equipment. Despite the extensive applications of this basic flow phenomenon, the design of these systems has been mainly based on trial and empirical methods. Hence, a detailed examination on the mass, momentum and energy transfer is essential for improving and optimizing the design of such industrial processes. However, due to the complexity of two-phase flow systems, standard formulation of the governing equations for bubble-driven flow has not yet been established. There is extensive research in the literature on the bubbly flow modeling and two-phase flow simulation. Some single-fluid models were proposed to take into account the existence of bubbles in liquid flows. These models considered the gas-liquid mixture as a single fluid with variable density caused by the distribution of bubbles in liquid. Since the properties of two fluids are combined into those of one mixture, these models cannot be used to predict the flow behavior of each phase. In the present state-of-the-art, the two-fluid model can be considered as the most detailed and accurate macroscopic formulation of the two-phase flow systems. In the two-fluid model, the gas and liquid phases are regarded as two interpenetrating phases, and each phase has its own set of conservation equations of mass, momentum and energy, coupled with some phase interaction terms. Such models receive the most attention in the literature since they require less knowledge about the gas distribution compared with single-fluid models, although some empirical correlations are also needed to estimate the interacting properties between the two phases.There are extensive experimental and computational studies available in the literature investigating bubbly flow hydrodynamics. In the recent study of Sherman et al. a two-fluid model coupled with population balance approach is used to numerical simulation of bubbly flow in a vertical isothermal channel¹. In their study, interfacial momentum transfer embraced various interfacial forces including drag, lift and wall lubrication. The turbulent dispersion force has also been accounted for appropriately. Durst et al. performed experimental studies on bubble-driven laminar flows by investigating liquid circulation and bubble's street with Laser-Doppler system. In their experiments, the flow was established in a glass cylinder of 100 mm i.d. with bubbles generated at a 0.5 mm dia. nozzle located at the centre of the bottom wall of a container. The cylinder was filled with castor oil of which the kinematics viscosity was sufficiently large to ensure liquid flow to be in the laminar regime. Air was passed through a well-controlled pressure regulator and the resulting bubbles left the nozzles at a steady rate of one for every 0.55 s with diameter db of about 5.0-6.5 mm. Two depths of liquid (H=0.098 m and H=0.278 m) were investigated to determine the profiles of liquid and bubble velocities. Durst et al. also conducted numerical simulations on Research Journal of Recent Sciences _____________________________________________________________ ISSN 2277-2502Vol. 4(1), 9-16, January (2015) Res.J.Recent Sci. International Science Congress Association 10 the bubble-driven flows at which both the single-phase model and the two-fluid model were used. In the two-fluid model, momentum equations for the liquid phase were derived from the Navier-Stokes equation, while equations for the gas-phase were simplified by neglecting the viscous forces³. Celik and Wang used a modified drift-flux model to simulate the problem with the same operating parameters as the experiments by Durst et al.It was reported that the liquid-phase circulation pattern is not sensitive to the actual shape of the void fraction profiles. Detailed predictions on liquid phase are given in their studies, while the corresponding ones for gas phase are not complete. Johansen et al studied the fluid dynamics in bubble-stirred ladles by employing a Laser-Doppler system to measure the axial and radial mean velocities of liquid phase. Air was supplied through a porous plug placed in the bottom wall of a cylindrical Perspex-water model of a ladle. However, no measurement on the gas phase was conducted. Johansen and Boysan proposed a mathematical model to describe the fluid flow in the bubble-stirred ladle at which Eulerian and Lagrangian methods were used to analyze the liquid and gas phases, respectively. A limited number of comparisons were conducted between simulations and experimental results of the liquid phase, but no information on the gas phase was reported. Castillejos and Brimacombe studied the plume formation and liquid circulation induced by gas injection at the bottom of cylindrical water containers. They developed an electro-resistivity probe to measure the detailed bubble characteristics such as gas volume fraction, bubble velocity and size. In addition, the axial and radial components of the liquid velocity surrounding the plume were measured by means of a Laser-Doppler system. Jinsong Hua et al. developed a modified two-fluid model and applied it to numerical investigation of bubbly flow in a cylinder with bottom gas injection, assuming rigid-surface spherical bubbles. In addition, they used the same operating parameters as the experiments by Durst et al.. In a range of two-fluid model, Mat et al. and Kemal Aldas used a two-fluid mathematical model for hydrogen production in a forced flow of the electrochemical system and found that the model successfully captures the main characteristics of the electrolysis process9,10. Tomomi Uchiyama simulated the bubbly flow around the cylinder by an incompressible two-fluid model11,12. A finite element method is applied to solve the model in conjunction with the Arbitrary Lagrangian–Eulerian (ALE) scheme considering the moving boundary of Uchiyama13. He found that the increase of the damping ratio in the bubbly mixture is attributable to the phase lag of the drag force acting on the cylinder behind the cylinder displacement. A numerical study is carried out to comprehend growth and movement of spherical bubble during foam processing by Gregory Rosebrock et al The numerical model is based on a level set technique for capturing the phase interface. It was found that the increase of initial bubble radius causes an increase in both the bubble growth rate and movement14. The dynamical processes of a Newtonian spherical drop rising freely through shear-thinning fluids expressed by the generalized Cross-Carreau model were considered experimentally and computationally by Mitsuhiro Ohta et al. The local effects of shear-thinning on the drop motion, which are hard to evaluate from an experimental approach, are clearly revealed by the numerical results15. Chen et al.applied a non-spherical model for bubble formation coupled with phase change at a submerged nozzle in a flowing sub cooled liquid16. The model was attempted to calculate the instantaneous shape of the steam bubble during its formation and to determine the bubble size at detachment, frequency of bubble formation, total steam flow rate as well as steam condensing rate. The major aim of present study is to investigate the sensitivity of the numerical results to the sphericity assumption for the bubbles in bubbly flow regimes. Therefore, we use a two-fluid mathematical model to allow the mass and momentum transfer between two phases in the numerical simulation of bubbly flow in a vertical cylinder. Experimental studies available in the literature have shown that typical sizes of the bubbles formed in the cylinder are in the millimeter regime and their shapes, during formation are generally non-spherical2,16. Typically, a bubble starts off as a hemisphere, and then becomes approximately spherical as it expands, and finally attains the shape of an irregular ellipsoid with a neck. Due to the large equivalent diameters of bubbles (db=6mm), they are of high Weber number, meaning that they are deformable and their dynamics is dependent upon their surface tension. The large bubbles have a high tendency of breaking up due to their low surface tension 17, but we neglect bubble break-up in the present work. Both the spherical and non-spherical bubbles are investigated. Therefore, the two methods are applied to calculate the drag coefficient. The first method, which has been used in many researches, is the empirical correlation of Simonin et al., obtained from the standard drag curve18. This method is suitable for rigid surface bubbles that remain spherical in shape and their surface-tension effects are negligible, such as small bubbles (db1mm) moving in contaminated fluid, or fluid containing surfactants 19,20. An alternative method that assumes a non-rigid surface is the "dirty water" model of Wallis et al. 21. This model covers the complete range of Reynolds numbers, Weber numbers, and the various shapes of bubbles22,23. Finally, the predicted results have been compared with the experimental results of Durst et al. It is found that the numerical results computed by the "dirty water" model are closer to those of the experimental data2.. Methodology Physical problem: The schematic sketch of problem considered is given in Fig. 1. In this figure, r and z represent the radial and axial coordinates, respectively; v and w stand for the velocity components in radial and axial directions, respectively. In such a system, air is injected into a liquid bath (with a depth Research Journal of Recent Sciences _____________________________________________________________ ISSN 2277-2502Vol. 4(1), 9-16, January (2015) Res.J.Recent Sci. International Science Congress Association 11 of H and a radius of R) from an orifice located at the bottom wall of a cylindrical vessel at a steady state flow rate. The cylinder is filled with castor oil of which the kinematic viscosity 10699.0) is sufficiently large to ensure liquid flow to be in the laminar regime. The actual dimensions of the liquid container and gas injection rates are specified according to the available experimental data of Durst et al.. The configuration parameters for the simulation are listed in able1. Figure-1 Schematic diagram of the problemTable -1 Physical configuration for the laminar bubble-driven flow Radius of cylinder (R) 0.1 m Bubble diameter (db) 6*10 - 3 m Liquid depth (H) 0.098 m Gas injection rate (Q) 2.6*10 - 7 m3/s Density of liquid ( r ) 960.3 kg/m3 Gas density ( r ) 1.29 kg/m3 Viscosity of liquid ( m ) 0.6712 Ns/m2 Viscosity of gas ) 1.64*10 - 5 Ns/m2 Governing equations: For numerical analysis, the following assumptions were utilized: i. The mass transfer between the liquid and gas phases exist. ii. The liquid and gas phases are treated as incompressible fluid. iii. The coalescence and break-up of the bubbles are not considered. iv. No turbulence model is employed. v. The bubbles are assumed to be in the same size. In order to represent the flow behavior in the system, a two-phase mixture of the liquid and gas was considered. The phases were assumed to share some spaces in proportion to their existence probabilities in such a way that their sum of the volume fractions would reach to a unity in the flow field. This assumption can be expressed mathematically as: (1) ,1 Where are the volume fractions of the liquid and the gas, respectively. The zone-averaged quantities are obtained through the solution of separate transport equations for each phase. Within this framework and above assumptions, the governing equations for the two-phase flow system can be expressed in cylindrical coordinates and conservative form as follows: Mass conservation: (2) int Where subscript, i represents the phases and takes the value of l and g refer to the liquid and the gas phases, respectively, in this and subsequent formulations. The term on the right hand of the above equation represents mass diffusion between the two phases at the liquid–gas interface which is calculated as (3) int Whererepresents the molecular diffusivity of phase i and may be expressed as 9,10,11: (4) eff m = Where eff is the effective viscosity of phase i and expressed as 9,10,11: (5) eff m m = Where m is the thermodynamic viscosity of phase i. r-momentum: (6) 2 ()()() ()()() ()(),vvvwiiiiiiiiiitrz i pFvvijiieffrrrvvviiiiieffeff ii zzrr r rararaaamamam¶¶¶ ++= ¶¶¶¶¶ -+-++ ¶¶¶¶¶+- ¶¶¶ z-momentum: (7) ()()() ()()() ()(),iiiiiiiiiiiiirjiiiieffieffieffwvwwtrzpFwwg zrr ww zzrz rararaaraamamam¶¶¶++=¶¶¶¶¶ -+-+++ ¶¶¶¶¶ ¶¶¶ Research Journal of Recent Sciences _____________________________________________________________ ISSN 2277-2502Vol. 4(1), 9-16, January (2015) Res.J.Recent Sci. International Science Congress Association 12 is the buoyancy force where g being the gravity vector. in both momentum equation is interface friction term and represents momentum exchange between the phases per unit volume. The phase interaction forces may include drag forces, lift forces and virtual mass forces. Only the effect of drag forces is considered in this study. Interface friction term can be expressed as: (8) Re, Where is the drag coefficient. A large number of formulations exist for estimating the drag coefficient. In this study both spherical and non-spherical bubbles are investigated. Therefore, the two methods are used to calculate the drag coefficient. The first of these methods, which has been used in many researches, is the empirical correlation of Simonin et al. 18, obtained from the standard drag curve.Spherical bubbles assumption: (9) Re425001(42.0Re15.01(Re0.2416.1687.0 The method is suitable for rigid surface bubbles that remain spherical in shape and their surface-tension effects are negligible, such as small bubbles (db1mm) moving in contaminated fluid, or fluid containing surfactants19,20An alternative method, that assumes a non-rigid surface and allows for the various shapes of bubbles, uses a range of equations dependent on the Reynolds number, Re, and the Weber number, We. These non-dimensional numbers are calculated from the equations:Non-spherical bubbles assumption: (10) ReWe Where is the slip velocity vector between the two phases, g is the bubble surface tension and being the equivalent spherical diameter of a non-spherical bubble defined as: (11) ,)6(31bbVdp= Where is the volume of a bubble. The equations associated with the model are: (12) 100Re,8,3/81.2065Re,8,3/1.2065Re,8Re/3.6100ReRe68.2049.0ReRe,166.26.2385.0643.0WeWeWeWeWeWe The above equation is the "dirty water" model of Wallis et al.. The drag coefficient for this method is calculated using21(13) 1(7.1 here is the volume fraction of bubble swarm. The multiplier of in the above equation, termed the swarming coefficient, accounts for the reduction in drag due to bubble-bubble interactions. Under the assumption that the bubbles retain their spherical shape at all times then, eq. 9 may be used to calculate the drag coefficient. But, for the case of non-spherical bubbles the dirty water model (eq. 12) is used to calculate the drag coefficient. Boundary conditions: The conditions for the dependent variables need to be specified along the boundaries of the computational domain. Due to the constraints of axis-symmetry, along the centerline of vessel (r=0), the following conditions are applied:(14) ,0,0 (15) ,0 (16) ,0 For laminar bubble-driven flows, we can use the no-slip condition for the boundary conditions along the solid walls. In addition, the gradients of volume fractions of the gas and the liquid along the walls vanish because no mass flux can pass through the walls. These conditions are given by: (17) ,0,0 (18) ,0 Where new points to the direction normal to the solid wall. For simplicity, we neglect the formation of waves on the liquid surface by assuming that the liquid surface is flat. In addition, Research Journal of Recent Sciences _____________________________________________________________ ISSN 2277-2502Vol. 4(1), 9-16, January (2015) Res.J.Recent Sci. International Science Congress Association 13 no liquid flow is allowed to pass through the top boundary, while the gas leaves the surface at a rate given by the natural boundary condition:(19) ,0;0 For the laminar flow, the numerical results of Durst et al. shows that the case with zero liquid velocity both along and normal to the surface [as given in equation (19)] agrees better with the experimental results. They attributed this fact considering that bubbles could cause an air-liquid interface to act as a partially rigid surface. The boundary conditions on the gas injection orifice and the gas injection rate are given by the available experimental data of Durst et al.. Numerical method: The governing equations for the two-fluid model are solved using the appropriate algebraic equations derived by integrating over the meshed control volume. For this purpose, the computational domain is overlaid with an orthogonal uniform grid with the nodal points of velocity variables for each phase being staggered with respect to those of pressure and void fraction, located at the regular grid nodes. The governing equations are integrated over a control volume by which flow domain is divided and the following algebraic equation is obtained:(20) In which a’s represent convection-diffusion coefficient and subscripts W, E, S, N represent west, east, south and north of node P and also b represents the volumetric fluid interaction coefficient. In addition, the power-law scheme is employed to approximate the convection-diffusion terms. The whole set of equations is solved by incorporating the SIMPLE algorithm. Because of the nonlinearity of the governing equations and the coupling of variables, iterative numerical procedures are conducted until the convergence is reached. Furthermore, grid sensitivity is examined to ensure sufficient spatial resolution. The results reported are based on the mesh of 30×15 uniform grids. Results and Discussion Predicted results include the void fraction and velocity fields, computed for both the spherical and non-spherical bubbles, which are compared with experimental data measured by Durst et al.The simulated axial gas velocity along the centerline of vessel (normalized by the maximum bubble axial velocitymax), which is compared with the experimental data of Durst et al. , is shown in Fig. 2. Near the gas entrance, the bubbles are driven by the buoyancy force and move upwards. After a certain distance from the leading point, bubbles reach a terminal velocity and move up at this speed until they reach the vicinity of the free surface where the bubble velocities are reduced due to the increase of drag force caused by increasing relative motion between bubble and liquid phases. It should be noted that the liquid flow is a result of friction between two phases. Therefore, when the liquid approaches the free surface, the axial velocity decreases while the radial velocity increases. As a result, the liquid shifts in the radial direction and may cause the liquid circulation in the container. It is seen that the gas velocity computed by the non-spherical assumption is closer to that of the experimental data. Figure-2 Gas phase axial velocity along the centerline of vessel. 0.20.40.60.81.200.20.40.60.81z/H/Wg,max Experimental datas of Durst et al [2] Results for non-spherical bubbles Results for spherical bubbles Research Journal of Recent Sciences _____________________________________________________________ ISSN 2277-2502Vol. 4(1), 9-16, January (2015) Res.J.Recent Sci. International Science Congress Association 14 Figure-3 Predicted vertical component of liquid velocity at the vertical section of z/H=0.36 Figure-3 shows the comparison between the simulated profiles of axial liquid velocity at the vertical section of) z/H=0.36) and the experimental data of Durst et al. . Here, the maximum bubble velocitymax, is used to normalize the liquid velocity. It is seen that the liquid velocity computed by the non-spherical assumption is closer to that of the experimental data's. The figure shows that the maximum liquid velocity is located on the centerline of the plume where the maximum gas volume fraction is formed. The significant decrease of axial liquid velocity along the lateral direction may be attributed to the decrease of bubble void fraction due to the weakness of bubble diffusion24,25The predicted gas void fraction profiles at two vertical sections in the cylinder, compared with the experimental data of Durst et al. are presented in figures. 4, 5. It is seen that the maximum gas volume fraction is located in the centre of the cylindrical container. Most of the bubbles are located in the central region of the container, and bubble volume fraction near the gas injection orifice is higher than that far away from the orifice. The significant decrease of gas void fraction along the lateral direction of cylinder may be attributed to the decrease of bubble diffusion. It is seen that the gas void fraction computed by the non-spherical assumption is closer to that of the experimental data. It should be noted that the gas volume fractions for 2.0 R are two orders of magnitude lower than the central regime and hence their value cannot be observed on a linear scale shown in these figures. Figure-4 Predicted gas void fraction distribution at the vertical section of z/H=0.06. -0.10.10.20.30.400.20.40.60.81r/R/Wg,max Experimental datas of Durst et al [2] Results for non-spherical bubbles Results for spherical bubbles 0.050.10.150.20.2500.050.10.150.2r/Raa Experimental datas of Durst et al [2] Results for non-spherical bubbles Results for spherical bubbles Research Journal of Recent Sciences _____________________________________________________________ ISSN 2277-2502Vol. 4(1), 9-16, January (2015) Res.J.Recent Sci. International Science Congress Association 15 Figure-5 Predicted gas void fraction distribution at the vertical section of z/H=0.36 Conclusion Numerical simulation of bubbly flow caused by the axial gas injection into a vertical cylinder is studied in this paper. It is found that the radial components of bubble velocities are much smaller than the corresponding axial ones, and the bubbles move up almost vertically. The maximum liquid velocity is located on the centerline of the plume where the maximum gas volume fraction is formed. The maximum gas void fraction is located in the centerline of the container. Most of the bubbles are located in the central region of the container, and gas volume fraction near the gas injection orifice is higher than that far away from the orifice. The major aim of this study was to investigate the sensitivity of the numerical results to the sphericity assumption for the bubbles. It is found that the numerical results computed by the "dirty water" model are closer to those of the experimental data. Therefore, taking into account the fact that the gas-liquid; inter-phase surface is mobile, the "dirty-water" drag model can be more suitable than the standard drag curve originally obtained for fixed inter-phase surfaces of liquid-solid or gas-solid systems. References 1.Sherman C.P. Cheung, Yeoh G.H. and Tu J.Y., On the modeling of population balance in isothermal vertical bubbly flows-Average bubble number density approach, International journal of Chemical Engineering and Processing, 46, 742–756, (2007) 2.Durst F., Taylor A.M.K.P. and Whitelaw J.H., Experimental and numerical investigation of bubble-driven laminar flow in an axisymmetric vessel, Journal of Multiphase Flow,10, 557-569 (1984) 3.Durst F., Schonung B., Selanger K. and Winter M., Bubble-driven liquid flows, International Journal of Fluid Mechanics,170, 53-82 (1986)4.Celik I. and Wang Y.Z., Numerical simulation of circulation in gas-liquid column reactor : Isothermal, bubbly, laminar flow, International Journal of Multiphase Flow,20, 1053-1070 (1994)5.Johansen S.T., Robertson D.G.C., Woje K. and Engh T.A., Fluid dynamics in bubble stirred ladles : Part I. Experiments, Metallurgical Transactions,19B, 745-754 (1988) 6.Johansen S.T. and Boysan F., Fluid dynamics in bubble stirred ladles : Part II, Mathematical modeling, Metallurgical Transactions,19B, 755-764 (1988) 7.Castillejos A.H. and Brimacombe J.K., Measurements of physical characteristics of bubbles in gas-liquid plumes : Part II., Local properties of turbulent air-water plumes in vertically injected jets, Metallurgical Transactions,18B,659-671 (1987) 8.Jinsong Hua and Chi-Hwa Wang, Numerical Simulation of bubble-driven liquid flows, International Journal of Chemical Engineering Scince,55, 4159-4173 (2000)9.Mahmut D. Mat, Kemal Aldas and Olusegun J. Ilegbusi, A two-phase flow model for hydrogen evolution in a electrochemical cell, International Journal of Hydrogen energy,29, 1015-1023 (2004)10.Mahmut D. Mat and Kemal Aldas, Application of a two phase flow model for natural convection in an electrochemical cell, International Journal of Hydrogen energy,30, 411-420 (2005)11.Kemal Aldas, Application of a two-phase flow model for hydrogen evolution in an electrochemical cell, Applied Mathematics and Computation,154, 507-519 (2004)12.Tomomi Uchiyama, Numerical prediction of added mass 0.050.10.1500.050.10.150.2r/Raa Experimental datas of Durst et al [2] Results for non-spherical bubbles Results for spherical bubbles Research Journal of Recent Sciences _____________________________________________________________ ISSN 2277-2502Vol. 4(1), 9-16, January (2015) Res.J.Recent Sci. International Science Congress Association 16 and damping for a cylinder oscillating in confined incompressible gas–liquid, two-phase mixture, Nuclear Engineering and Design,22, 68–78 (2003) 13.Uchiyama T., ALE finite element method for gas–liquid; two-phase flow including moving boundary based on an incompressible two-fluid model, Nuclear Engineering and Design,205, 69–82 (2001) 14.Gregory Rosebrock, Ahmed Elgafy, Thomas Beechem and Khalid Lafdi, Study of the growth and motion of graphitic foam bubbles, International Journal of Carbon,43, 3075–3087 (2005)15.Mitsuhiro Ohta, Eiji Iwasaki, Eiji Obata and Yutaka Yoshida, A numerical study of the motion of a spherical drop rising in shear-thinning fluid systems, International Journal of Non-Newtonian Fluid Mechanics,116, 95-111 (2003)16.Chen W.B., Reginald B.H. Tan, A model for steam bubble formation at a submerged nozzle in flowing subcooled water, International journal of heat and fluid flow,22, 552-560 (2001) 17.Thomas Nierhaus, David Vanden Abeele and Herman Deconinck, Direct numerical simulation of bubbly flow in the turbulent boundary layer of a horizontal parallel plate electrochemical reactor, International Journal of Heat and Fluid Flow,28, 1–10 (2007) 18.Study of Ionospheric Perturbations during Strong Seismic Activity by Correlation Analysis Method Gwal A.K., Jain Kumar Santosh, Panda Gopal and Gujar Y.S., Res. J. Recent Sci.,1(1), 2-9 (2012) 19.Giusti A., Lucci F. and Soldati A., Influence of the lift force in direct numerical simulation of upward/downward turbulent channel flow laden with surfactant contaminated micro bubbles, International journal of Chemical Engineering Science,60, 6176–6187 (2005)20.Empirical Correlation of Various Inclusions on the Effect of Primary and Secondary Parameters for Estimation of (ETC) of Two Phase Materials Senthil Kumar A.P., Karthikeyan P., Prabhu Raja V., Ramu M., Somasundharam S. and Vasudevan V., Res. J. RecentSci., 1(1), 22-32 (2012) 21.Wallis G.B., The terminal speed of single drop or bubbles in an infinite medium, International Journal of Multiphase Flow,1, 491, (1974)22.An Overview of Green Supply Chain Management in India Nimawat Dheeraj1 and Namdev Vishal, Res.J.Recent Sci.,1(6), 77-82(2012) 23.Effective Factors on Determination of Audit Fees in Iran Khani Hamid and Yazdani Qanbar Ali, Res.J.Recent Sci.,1(6), 38-44(2012)24.Estimation of Global Solar Radiation at Onitsha with Regression Analysis and Artificial Neural NetworkModels Agbo G.A., Ibeh G.F. and Ekpe J.E., Res.J.Recent Sci.,1(6), 27-31(2012)25.An Analysis of Malaysian Renewable energy target using Simulation Modelling Approach Akhwanzada, Salman Ahmad, Mat Tahar and Razman bin, Res. J. Recent Sci.,3(1),38-44(2014) Nomenclature drag coefficient molecular diffusivity coefficient, m 2 /s D equivalent spherical diameter, m b volumetric inter-fluid friction, kg/m 2 s 2 Reynolds number Re Weber number We r-velocity vector, m/s v z-velocity vector, m/s w buoyancy forces, kg/m 2 s 2 gravity vector, m/s 2 g volume, m 3 V liquid depth, m H radius of cylinder, m R Greek letters a r m g Subscript g l nb eff P void fraction density, kg/m 3 viscosity, N.s/m 2 bubble surface tension, N/m bubble gas phase liquid phase neighbor effective unknown node