UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 THERMOFLUID ANALYSIS OF TWO PHASE SLUG FLOW IN MICROCHANNELS WITH APPLICATION S TO LEAKAGE MODELING IN ROTARY VANE TURBOMACHINERY By AYYOUB MEHDIZADEH MOMEN A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY O F FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2010 PAGE 2 2 2010 A yyoub M ehdizadeh M omen PAGE 3 3 To my loving m other and f ather family and friends, without whose enduring patience and s upport this work would not have been possible PAGE 4 4 ACKNOWLEDGMENTS I would like to thank all those who have helped and supported me along the way. Special thanks go to my mentors (Dr. S. A. Sherif and Dr. William E. Lear). Their relentless time and guidance allowed me to scale to greater heights. I would like to thank my committee members for their guidance and input: Dr. James Klausner, Dr. Skip Ingley and Dr. Arnoldo Valle Levinson. I have been fortunate to have a diverse, knowledgeable committee that had answers to my questions, solutions to my problems, and great advice and encouragement. I would like to thank all the people working in Dr. Sherifs and Dr. Lears Labs, especially M rs. Fotouh Al Raqom and Dr. Ahmad Mahmoud, that have helped me throughout my years in graduate school, as well as make my long hours in the lab enjoyable. I would also like to thank all the staff that has assisted me greatly at UF (Department of Mechanical and Aerospace Engineering). I would like to thank my friends who have made my years in Gainesville more enjoyable: Miss. Romina Mozafarian, Mr. Babak Mahmoudi, Mr. Ashkan Behnam, Dr. Masi Rajabi, Miss. Moojan Daneshmand, Mr. Daniall Sabri, Miss. Sahar Mirshamshiri, Mr. Richard Parker, Dr. Farzad Fani, Miss. Sydni Credle, Mr. Cheng Chan Kuo Miss. Cinthia Perez, to name a few. Finally, I would like to thank the U niversity of Florida and the C ollege of Engineering for providing such a rich and competitive graduate program in which I am v ery lucky to have been a part. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ...................................................................................................... 4 LIST OF TABLES ................................................................................................................ 8 LIST OF FIGU RES .............................................................................................................. 9 NOMENCLATURE ............................................................................................................ 12 ABSTRACT ........................................................................................................................ 17 CHAPTER 1 GENERAL BACKGROUND AND OBJECT IVES ...................................................... 19 Two Phase Flows in Microchannels .......................................................................... 19 Objective ..................................................................................................................... 20 2 INTRODUCTION A ND LITERATURE REVIEW ........................................................ 22 3 SIMPLIFIED MODEL OF TWO PHASE SLUG FLOW IN MICROCHANNELS ....... 30 Outlook ........................................................................................................................ 30 Numerical Simulation .................................................................................................. 31 Single Phase Simulations .................................................................................... 31 Analytical Model of Gas Liquid Interface............................................................. 35 Introducing Moody Friction Factor for General TwoPhase Slug Flow .............. 37 First Attempt for Two Phase Simulations ............................................................ 40 Slugs Formation, Pressure Fluctuations and Main Frequencies ....................... 46 Conclusions .......................................................................................................... 49 4 CFD MODELING OF TWO PHASE SLU G FLOW IN MICROCHANNELS .............. 51 Outlook ........................................................................................................................ 51 Numerical Simulation .................................................................................................. 51 Problem For mulation and Solution Procedure ........................................................... 52 Governing Equations ........................................................................................... 52 Differencing Schemes .......................................................................................... 56 Interface Tracking ................................................................................................ 57 Model Geometries and Studied Test Cases ....................................................... 59 Results and Discussion .............................................................................................. 60 Bubble Shape ....................................................................................................... 60 Velocity and Pressure Field ................................................................................. 64 Volume Fraction ................................................................................................... 65 Gas Slug Velocity and Length ............................................................................. 66 PAGE 6 6 Liquid Film Thickness .......................................................................................... 67 Wall Shear Stress and Pressure Drop ................................................................ 68 Conclusions .......................................................................................................... 71 5 A COMBINED ANALYTICAL NUMERICAL MODEL FOR A TWO PHASE FLOW THROUGH A SUDDEN AREA CHANGE IN MICROCHANNELS ................ 72 Outlook ........................................................................................................................ 72 Analysis ....................................................................................................................... 72 Simplified Model of Two Phase Slug Flow Pressure Drop in Microchannel ...... 72 Flow Area Expansion Analysis, Single Phase .................................................... 75 Flow Area Contraction Analysis, Single Phase ................................................... 78 Two Phase Pressure Change across Area Change (Conventional Models) ..... 79 New Analytical Model for Void Fraction in Contraction ....................................... 83 New Pressure Drop Model for Slug Flow through a Contraction in Microchannels ................................................................................................... 92 Conclusions ........................................................................................................ 102 6 NUMERICAL SIMULATION OF THERMOFLUID CHARACTERISTICS OF TWO PHASE SLUG FLOW IN MICROCHANNELS ............................................... 103 Outlook ...................................................................................................................... 103 Introduction ............................................................................................................... 104 Problem Formulation and Solution Procedure ......................................................... 105 Governing Equations: ........................................................................................ 105 Test Cases Examined ........................................................................................ 109 Model Verification and Validation ...................................................................... 110 Results and Discussion ...................................................................................... 112 Conclusions ........................................................................................................ 123 7 APPLICATIONS OF TWO PHASE SLUG FLOW IN MICROCHANNELS TO ESTIMATING INTERNAL LEAKAGE IN ROTARY VANE TURBOMACHINERY .. 125 Outlook ...................................................................................................................... 125 General Background ................................................................................................. 126 Refrigeration Cycle............................................................................................. 126 Thermodynamic Overview on Refrig eration Cycle ........................................... 126 Potential of Increasing Coefficient of Performance .......................................... 128 Microscale Leakage Paths in Rotary Vane Expander ...................................... 130 Analytical Model ........................................................................................................ 131 Geometrical Analysis ......................................................................................... 131 Vane Friction ...................................................................................................... 134 Torque and Power .............................................................................................. 137 Ideal Volumetric Flow Rate ................................................................................ 140 Leakage .............................................................................................................. 141 Efficiencies ......................................................................................................... 147 PAGE 7 7 Two Phase Internal Leakage Estimation .......................................................... 150 Conclusion .......................................................................................................... 151 8 RECOMMENDATIONS AND FINAL CONCLUSIONS ............................................ 153 LIST OF REFERENCES ................................................................................................. 156 BIOGRAPHICAL SKETCH .............................................................................................. 162 PAGE 8 8 LIST OF TABLES Table page 3 1 Coefficients in equation (35) ................................................................................. 40 4 1 Test cases examined ............................................................................................. 60 4 2 Comparison of velocity of gas slugs for different test cases with the available empirical correlation and previous numerical simulations .................................... 67 4 3 Comparison of liquid film thicknesses with those calculated by available empirical correlations and Gupta et al .s numerical simulation ............................ 68 4 4 Comparison of pressure drop with available empirical correlations and previous numerical simulations .............................................................................. 70 6 1 Property values employed ................................................................................... 109 6 2 Test cases examined ........................................................................................... 110 7 1 Leakage flow rate from the clearance between the rotor stator cylinders for different pressures ................................................................................................ 144 7 2 Leakage flow rate (mean) from the clearance between the side of the vanes and end plates ...................................................................................................... 151 PAGE 9 9 LIST OF FIGURES Figure page 1 1 Schematic of different possible twophase flow regimes in a hor izontal pipe ...... 19 2 1 Gas liquid multiphase flow pattern in a microchannel, Serizawa et al. (2002) on 20m ID pipe (top). Triplett et al. (1999) on 1.097 mm ID pipe (bottom) ........ 23 3 1 Experimental image of a typical slug flow in a microchannel (Serizawa et al. 2002) ....................................................................................................................... 31 3 2 Axial velocity contours and streamlines for a 2mm pipe flow for two liquid slug aspect ratios: Top is 1 mm liquid slug length and bottom is 5mm liquid slug length .............................................................................................................. 33 3 3 Non dimensional wall shear stresses over the liquid slug lengths and the averaged she ar stresses for a 2 mm diameter pipe, Re=2000 (slug moves to the left) .................................................................................................................... 34 3 4 Calculated air water interface profile, Z(r) for a 2 mm diameter microchannel, and ........................................................................ 37 3 5 Axial velocity contours and streamlines (top). Nondimensional wall shear stress over slug length and the average shear stress for a 2mm di ameter pipe (bottom) .......................................................................................................... 38 3 6 Averaged non dimensional shear stress for different slug lengths and Reynolds numbers ................................................................................................. 39 3 7 Axisymmetric computational domain (top half is used in numerical simulations, while the lower half is a mirror image of the top half) ....................... 42 3 8 Predicted slug flow regime in a microchannel for a water air system .................. 47 3 9 (left) Time sequence of slug formation close to the T junction (superficial velocities. (right) Pressure fluctuations in the mixing region ................................. 49 4 1 Repea ted). Gas liquid multiphase flow pattern in a microchannel, Serizawa et al. (2002) on 20m ID pipe (top); Triplett et al. (1999) on 1.097 mm ID pipe (bottom) .................................................................................................................. 53 4 2 Schematic showing the test ca ses and boundary conditions (not to scale) ......... 5 9 4 3 Volume fraction contours for the course meshes. ................................................. 61 PAGE 10 10 4 4 Comparison of volume fraction contours with the experimental images of Triplett et al. (1999) on a 0.5 mm pipe with refined mesh i n the vicinity of the interface .................................................................................................................. 63 4 5 Simulation results for Run 12 (dynamic refined mes h at the interface employing 15,639 elements) [top]; Simulation results of Gupta et al. (2009) (employing 448,000 elements) [bottom] ................................................................ 65 4 6 Velocity vectors and pressure contours in a frame of re ference moving with the gas slug (Run 7) ............................................................................................... 65 4 7 Pressure distribution on wall and axis of the microchannel for Run 12 (t=0.07887 s). ......................................................................................................... 69 5 1 Left: Experimental image, Serizawa et al. (2002). Right: Schematic picture of two phase slug flow in a straight microchannel .................................................... 73 5 2 Schematic of single phase flow through an expansion ......................................... 75 5 3 Schematic of single phase flow in Contraction ...................................................... 79 5 4 Schematic of twophase slug flow in a microchannel facing a contraction .......... 84 5 5 Summary of algorithm used for modeling the average void fraction in the vicinity of the flow area contraction ........................................................................ 89 5 6 Gas slug pressure, liquid center of mass velocity and liquid front velocity .......... 90 5 7 Two phase flow regime map for a 1.1 mm diameter tube, Triplett et al. (1999) Bold big circles are Abdelall et al. s (2005) data w hich were used for the analytical modeling ................................................................................................. 95 5 8 Relationship between homogeneous void fraction and void fraction for different models in the vicinity of the flow area contraction .................................. 96 5 9 Comparison of six different void fraction models a) New analytical model (Equations (534) to (5 45)) b) Armand correlation c) Kawahara et al.s correlation d) Homogeneous model e) Pure curve fit of experimental data f) Zivis correlation. .................................................................................................... 98 5 10 Nondimentional pressure drop versus twophase Reynolds number for the experimental data of Abdelall et al. (2005) .......................................................... 100 6 1 Schematic of thermal boundary layer development for plug flow (left) and Poiseuille flow (right) ............................................................................................ 104 6 2 Comparison of simulated wall temperature with experimental data ................... 111 6 3 Temperature rise along the microchannel (Test Case 1) ................................... 113 PAGE 11 11 6 4 Temperature contours, velocity vectors and dimensionle ss wall heat flux for a train of a liquid and gas slug at the stationary frame of reference ..................... 114 6 5 Comparison of the averaged and local Nusselt numbers of a slug flow in a microchannel ........................................................................................................ 117 6 6 Wall temperature variation of the solid heater for several times, left. Convective heat transfer coefficient in different locations, right ......................... 121 6 7 Dimensionless temperature in the solid wall for various dimensionless times .. 123 7 1 A schematic of a typical refrigeration system ...................................................... 126 7 2 Temperature vs. Entropy diagram for a typical refrigeration cycle ..................... 127 7 3 Various components of a rotary vane expander machine .................................. 129 7 4 Schematic of twophase flow passing through a microscale leakage path at the tip of a vane .................................................................................................... 130 7 5 The tested rotary vane air motor (a and b) and the schematic of the device .... 131 7 6 The calculated variation of cross sectional area of a chamber during one revolution ( m2) ...................................................................................................... 134 7 7 Schematic of forces acting on a v ane (not to scale) ........................................... 135 7 8 Variation of normal forces acting on a vane as a function of (deg) ................ 136 7 9 Single van e torque during a revolution ................................................................ 138 7 10 Measured and the calculated torque and power for 2.8 bar inlet pressure ........ 139 7 11 Measured and calculated (ideal) volumetric flow rate at the outlet at Pi=2.8 bar ......................................................................................................................... 142 7 12 Schematic of leakage paths ................................................................................. 143 7 13 The Measured and the modeled volumetric flow rate and leakage (Pi=2.8 bar) ........................................................................................................................ 145 7 14 The calculated leakage ratio for different inlet pressures ................................... 146 7 15 The c alculated efficiencies at Pi=2.8 bar ............................................................. 149 PAGE 12 12 NOMENCLATURE A pipe cross sectional area m2 A a rea of throat m2 2 1, C C e mpirical constants Ca Capillary number Co Courant number CpG gas specific heat J/Kg.K CpL liquid specific heat J/Kg.K D microchannel diameter, m d eccentricity m E energy term J f Moody friction factor SFF Continuum surface force vector, N/m3 ff wall shear stress Pa G total mass flux kg/m2.s g gravitational acceleration, m/s2 Gr Graetz number h thickness of vane housing, m GJ g as superficial velocity, m/s LJ l iquid superficial velocity, m/s k s urface curvature PAGE 13 13 Kc contraction loss coefficient Kd momentum correction factor Ke expansion loss coefficient L length of vane, m l length of the leakage path, m LR leakage ratio gm g as slug mass kg lm l iquid slug mass kg tm. t otal mass flow rate kg/s n number of vanes n s urface normal unit vector p static pressure P a Pe Peclet number Pi inlet pressure P a PO ambient pressure, P a Pr Prandtl number ps perimeter of microchannel channel m2 Qi ideal volumetric flow rate (no leakage) m3/s QL leakage of rotor stator c learance gap, m3/s Ql leakage through the blades side, m3/s ave wq a verage heat flux W/m2 x wq l ocal heat flux W/m2 PAGE 14 14 R radius of microchannel, m R gas constant j/K.mol ReD Reynolds number Ri stator radius m RO rotor radius m S slip ratio t time, s T temperature, K Tm mean temperature, K T0 inlet temperature, K Tv torque N.m t t ime step, s LSU l iquid slug velocity m/s GSU g as slug velocity m/s gu g as slug velocity m/s l u l iquid slug velocity m/s fu f ront slug velocity m/s u(x,y) velocity of flow in the leakage path, m/s v v elocity vector m/s W width of the chamber m x vapor mass quality PAGE 15 15 Xc liquid slug center of mass void fraction L l iquid phase void fraction G gas volume fraction homogeneous g as void fraction l iquid film thickness, m r o ughness m c oefficient of surface tension, N/m a a rea ratio L d ensity of liquid phase, kg/m3 G d ensity of gas phase, kg/m3 m ixture density, kg/m3 m ixture viscosity, kg/m.s l l iquid phase viscosity kg/m.s g g as phase viscosity kg/m.s c ontact angle cir c irculation time scale r stator curve radius to origin, m a ngular velocity of rotor rad/s m m echanical efficiency PAGE 16 16 s i sentropic efficiency o verall (to tal) Efficiency PAGE 17 17 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy THERMOFLUID ANALYSIS OF TWO PHASE SLUG FLOW IN MICROCHANNELS WITH APPLICATIONS TO LEAKAGE MODELING IN ROTARY VANE TURBOMACHINERY By A yyoub M ehdizadeh M omen May 2010 Chair: S. A. Sherif Co chair: W illiam E Lear Major: Mechanical Engineering Understanding the characteristics of twophase slug flow inside microchannels is a challenging problem that has many applications in chemical mechanical and biological processes. Even though many experimental studies focused on two phase flow s in microchannels, considerable discrepancies among the d ata still exist, mainly due to the difficulties in setting up the experiment s and/or in conducting the measurements. Similarly, numerical modeling of these flows has failed to accurately capture the flow physics especially those pertaining to the characte ristics of the slug flow regime. The presence of a thin liquid film on the order of 10 m around the bubble may be a contributing factor to the above difficult ies This study describes slug flow hydrodynamics and heat transfer in microchannels and tries to overcome the aforementioned modeling difficulties. T he effects of the surface tension on hydrodynamics and dimensionless wall shear stress are investigated and a general Moody friction factor for a typical twophase slug flow in microchannels is define d. PAGE 18 18 Formation of the slugs for different superficial velocities, c apillary numbers, and gas volume fractions are investigated with the help of Volumeof Fluid (VOF) method Innovative analytical models f o r the calculat ion of two phase slug flow pressure dr op in microchannels encountering a sudden contraction are introduced. I n contrast to the existing empirical correlations which were tuned for a specific geometry, the new analytical models are capable of taking geometrical parameters as well as flow condit ions into account A fundamental study of heat transfer characteristics of twophase slug flow in microchannels showed that the convective heat transfer coefficient in the film region is significant. A strong coupling between the conductive heat transfer i n the solid wall and the convective heat transfer in the flow field is observed and characterized. Local and time averaged Nusselt numbers for slug flows are reported for the first time. Significant improvements in the heat transfer coefficient could be ac hieved by short slugs where the Nu number was found to be 610% higher than in single phase flows. Understanding the physics of twophase flows in micro scale leakage paths in any devices is one of the applications of this study. This study demonstrates th at large internal leakage that inherently exists in rotary vane expanders undermines its functionality. The developed models in this study employed to estimate the internal leakages in a rotary vane expander. PAGE 19 19 CHAPTER 1 GENERAL BACKGROUND AND OBJECTIVES Two Phase Flows in Microchannels Pipes with hydraulic diameters smaller than the Laplace constant, ) (G Lg g is gravitational acceleration and L and G are the liquid and gas densities, respectively) are typically classified as microchannels. The above threshold for an air water system in the room temperature is about 2.5 mm. Figure 11. Schematic of different possible twophase flow regimes in a horizontal pipe Figure 11 shows different schemes of two phase flow regimes in a horizontal pipe. Different regimes occur based on the flow conditions and fluids properties In PAGE 20 20 typical large pipes the characteristics of any regimes are governed by contributions of g ravitational, surface tension and frictional forces However, such characteristics in microchannels are mainly governed by the surface tension force with the gravitat ional acceleration playing only a minor role on the flow field. Slug flow is a typical twophase regime that occupies a large area on the flow regime map in microchannels. Two phase flow in microchannels has wide applications in many technologies. Examples include electronic micro chips, capillary micro reactors biology microelectromechanical systems (MEMS) and leakage modelings Objective The main objective of this study is to investigate the characteristics of two phase slug flows in microchannels. Sev eral features of this study help for a better understand ing of the flow physics which was not already discovered by other researchers. This fundamental study has wide applications in many technologies. In order to achieve to the above objective e xtensive studies carried out to understand twophase flow physics in microchannels. A detailed list of the tasks done in this dissertation is as follows: Develop a simplified model of two phase slug flow in microchannels. Investigate the effect of the air water interface on the hydrodynamics of the slugs by solving Young Laplace equation. Introduce dimensionless form of wall shear stress and defin e a new general Moody friction factor for a typical twophase slug flow. Find the main frequencies of pressure fluctuation generated by slug breakup at the axisymmetric T junction using the Fast Fourier Transform. Utilize CFD simulation using Volumeof Fluid method to captur e slug flow physics in different conditions PAGE 21 21 Investigate the flow pattern in the thin liquid film (in the order of 10 m ) around the bubble Introduce the strategy of dynamic mesh generation in order to have a successful simulation by which the liquid film could be captured utilizing the minimum number of meshes Develop innovative analytical models to calculate the two phase slug flow pressure drop in microchannels encountering a sudden contraction. Investigate heat transfer characteristics of twophase slug flow in microchannels. Introduce an existence of a strong coupling between the conductive heat transfer in the solid wall and the convective heat transfer in the flow field. Introduce local and timeaveraged Nusselt numbers for slug flows in the microchannels. The study of the slug flow in microchannel helps to investigate internal leakage of two phase flows in micro scale leakage paths. The internal leakage and the characteristics of a rotary vane expander are introduced in Chapter 7 as an application of this study. This task is done by combining experimental, analytical and numerical modeling. A detailed list of the tasks done in this part of study is as follows: Develop thermodynamic and fluid mechanics model of a rotary vane expander. Conduct a series of experiments on the rotary vane expander in order to understand the working principles and operating constrains of the device. Apply the knowledge obtained from the work done on two phase flow in microchannel to model the internal leakage in the rotary vane expander. PAGE 22 22 CHAPTER 2 INTRODUCTION AND LITERATURE REVIEW Slug flow can be found in many applications involving smallscale twophase fluid flows. It generally occurs when two streams form alternating pockets of fluid (called slugs) in a channel under the action of surface tension forces. These slugs may merge or break up because of surface tension, pressure fluctuations, and wall adhesion as they move through the channel. The study of flow in microchannels has become of greater interest in recent years mainly due to their presence in a broad array of applications such as microelectromechanical systems (MEMS), electronics cooling, chemical process engineering, medical and genetic engineering, and bioengineering. Pipes with hydraulic diameters smaller than the Laplace constant, ) (G Lg liquid and gas densities, respectively) are typically classified as microchannels. In macrochannels, Taylor instabilities control many critical interface processes, while in microchannels hydrodynamic interface proce sses are not governed by the above instabilities (Triplett et al 1999). The characteristics of two phase flows in microchannels are mainly governed by the surface tension force with the gravitational acceleration playing only a minor role on the flow fiel d (Galbiati and Andreini, 1994; Triplett et al 1999; Brauner, 1990; Fukano and Kariyasaki, 1993). The absence of the stratified twophase flow regime in microchannels was observed in many previous experiments (Galbiati and Andreini, 1994 and Triplett et al 1999). Slug flow (sometimes called Taylor bubble flow, bubble train flow, intermediate flow or plug flow) is a typical PAGE 23 23 two phase regime that occupies a large area on the flow regime map in microchannels [Triplett et al 1999 and Serizawa et al 2002 ( s ee Figure 2 1 )]. Figure 2 1. Gas liquid multiphase flow pattern in a microchannel, Serizawa et al (2002) on 20m ID pipe (top). Triplett et al (1999) on 1.097 mm ID pipe (bottom) The flow circulation in the liquid slug significantly enhances radial heat transfer. Also, the presence of a gas slug between two consecutive liquid slugs reduces axial liquid mixing (Angeli and Gavriilidis, 2008). Two phase hydrodynamic characteristics in microchannels have been found to be different than those in large channels. Recent experimental studies by Ghiaasiaan and Abdel Khalik (2000), Abdelall et al (2005), Toufik et al (2008), Kawahara et al (2002), PAGE 24 24 Serizawa et al (2002) and many others tried to formulate and monitor flow in microchannels. However, there are still discrepancies in the data, mainly due to the difficulties in experimental setup and measurement. On the other hand, there are very few analytical or numerical studies in twophase flows in microchannels mainly due to the lack of th e detailed experimental data or robust physical based models for simulation. He and Kasagi (2008) simulated a single bubble in a micro tube. Fukagata et al. (2007) numerically simulated two phase flow in a micro tube and found that the gas liquid slip rati o is approximately 1.2 and this is in accordance with the Armand correlation that is valid for twophase flows in micro sized channels. De Schepper et al. (2008), utilizing CFD simulation, investigated the performance of the existing numerical tools and approaches for modeling of twophase flows. Their qualitative comparisons between computed contours and the experimental photos showed that simulation could capture two phase flow regimes except for slug flow. They tried all different available two phase mod els to overcome this problem; however, their simulation failed in capturing slug flow regime both qualitatively and quantitatively. They presumed that this fact is due to the presence of a small region of slug flow in the Baker chart and because of this fa ct, unpredictable transition between regions could affect their simulations. Understanding slug hydrodynamics and its internal structure has significant engineering value. Even though slug flow appears to be a relatively simple physical phenomenon, there are still considerable unresolved questions about the physics and boundary conditions. These include the dynamics of a moving contact line among three materials (two fluids in contact with a solid boundary, Darhuber and Troian 2005). Also, PAGE 25 25 there is uncertai nty as to the applicability of the no slip boundary condition at the micro scale (Tretheway and Meinhart 2002; Granick et al. 2003). The complexity of the flow physics increases as the geometry becomes more complicated. For example, many unsettled problems are associated with the slug flow pressure drop in a microchannel with an abrupt area change (Abdelall et al. 2005; Mehdizadeh et al. 2009a) To date there are no slug flow studies that address all of the a forementioned problems. Even if the above difficulties were overcome, analytical solutions for slug flow with many simplifying assumptions have only recently become available. For example, Bla ck et al. (2007) analytically solved a fourthorder ordinary differential equation of the stream function within the slug (moving wall axisymmetric cavity flow) at low Reynolds numbers. Mehdizadeh et al. (2009b) analytically showed that the slug length in a microchannel is an important parameter that affects the pressure drop in the sudden contraction. In Chapter 3 the conservation equations within a single rectangular liquid slug (as a simplified geometry of slug) are solved numerically. Then, the surface tension effects on the air water interface are analytically investigated employing the second order Young Laplace differential equation. The calculated interface shape is then used to define a more realistic geometry and a set of boundary conditions for numerical simulation purposes. Wall shear stress results are then compared to those in steady laminar duct flows by introducing a non dimensional form of the wall shear stress. Employing many numerical simulations, a correlation that predicts the Moody fric tion factor is defined. This correlation can in principle be used to calculate the slug flow PAGE 26 26 pressure drop in microchannels under different conditions. Moreover, the formation of slugs through the channel and their merging and breaking up are numerically predicted using the two phase volumeof fluid (VOF) method. Finally, the frequencies of the pressure fluctuations (generally difficult to measure experimentally in microchannels) in the mixing region are reported. Using Fast Fourier Transform (FFT) techniqu es, the main frequencies of the pressure fluctuations are extracted in order to correlate them to the slug lengths. Recently, discrepancies in connection with the presence of a liquid film surrounding the gas bubble in slug flows have been observed in vari ous CFD simulations. While some researchers admitted that the thin liquid film was not captured because of the need for a finer mesh (Qian and Lawal, 2006; Yu et al., 2007), others reported that the film simply does not exist. The latter group then had to study the effect of the threephase contact angle on the characteristics of the slug flow (He et al., 2007; Kumar et al. 2007). Gupta et al. (2009) reviewed the presence of wall dry out during experimental studies and compared the observed regime with tho se described in numerical simulations. They showed that a finer mesh is needed as nonphysical results could easily be produced if care was not exercised when modeling slug flows in microchannels. They compared five mesh resolutions and concluded that if m esh size were to be approximately of the liquid film thickness, capturing the liquid film can in fact be achieved. However, achieving this goal is extremely computationally costly (450,000 elements for their case). In Chapter 4, the evolution of slugs t hrough a microchannel and their merging and breaking up are numerically predicted using the twophase Volume of Fluid (VOF) PAGE 27 27 method. Employing a dynamic mesh adaption method with interface tracking one could significantly reduce the number of elements requi red for capturing the liquid film. Being able to reduce computational time (in compared to Gupta et al. (2009) simulation) allowed investigating the effects of different Capillary numbers and superficial velocities on the slug characteristics. This was ach ieved employing 12 test cases, with the results being in a good agreement with available experimental data. Two phase flow in microchannels with abrupt area change is among the least studied aspects of this type of flow. The work reported in Cha pter 5 provides a new analytical void fraction model in the vicinity of the area transition ( venacontracta location) for twophase slug flow and to apply this model to estimate the pressure drop due to the abrupt area change for twophase slug flow in mic rochannels. Compared to previous models that provide a static value for the void fraction without geometrical or flow constraints, the new model provides an analytical expression for the void fraction that considers both dynamic and static variables. Also, this model introduces a geometrical parameter that has historically been ignored in previous experimental work, which is the average liquid slug length and/or the gas slug length. The ratio of these two lengths can be calculated by knowing the void fracti on, but the actual length of each slug can be varying depending on the flow condition. Since the slug flow through the contraction behaves like pulsing flow, the actual length of each slug (in addition to the gas liquid slugs length) plays an important role on the calculated pressure drop. Moreover, it is expected that the actual length of the slugs controls whether or not the boundary layer is located in the developing region in each slug. PAGE 28 28 Heat transfer study of slug flow in microchannel is another aspec t of this flow which was rarely studied before. Walsh et al (2009) seem to be the only group of researchers who had performed experimental studies on heat transfer of slug flows in microchannels. In Chapter 6 we attempted to simulate the aforementioned ex periment employing the Volumeof Fluid (VOF) method. In Chapter 6 Heat transfer characteristics of the slug flow will be compared to that in the conventional Plug and Poiseuille flows where a constant heat flux boundary condition is applied. Local and tim e averaged Nusselt numbers for slug flows are reported for the first time. Significant improvements in the heat transfer coefficient could be achieved by short slugs where the Nu number was found to be 610% higher than in single phase flows. The study repo rted in Chapter 6 reveals new findings related to slug flow heat transfer in microchannels with constant wall heat flux. Understanding the physics of twophase flows in micro scale leakage paths in any devices is one of the applications of this study. The study reported in Chapter 7 demonstrates that large internal leakage that inherently exists in rotary vane expanders undermines its functionality. We have used models developed in the previous chapters to model the internal leakages in a rotary vane expan der. There are many studies on rotary vane turbomachinaries. Over the past two centuries, a large number of engine and compressor concepts have been offered. Rotary vane motors or compressors were offered in recent years due to their small size, low cost, mechanical simplicity and reliability, Peterson and McGahan (1972). Rotary vane expanders have been recommended as alternatives to turbines in low power (1050 horsepower) Rankine cycle applications, Robertson and Wolgemuth (1978). However, experimental re sults PAGE 29 29 have shown that losses due to leakage, friction and heat transfer must be considered to realistically judge the suitability of the vane expander to produce work, Eckard (1975), Robertson and Wolgemuth (1975). Even though rotary vane expanders are not yet promising alternatives for conventional power generation devices, rotary vane compressors are commonly used for automobile air conditioning systems, Li (2003). In Chapter 7 the advantages and the barriers to utilizing a rotary vane expander as a throttle valve replacement in vapor compression refrigeration systems is investigated. We have developed an analytical model to investigate frictional losses and performances in the device. Finally employing the knowledge we have gained from previous chap ters we will investigate internal leakage in a rotary vane turbomachine. PAGE 30 30 CHAPTER 3 SIMPLIFIED MODEL OF TWO PHASE SLUG FLOW IN M ICROCHANNELS Outlook In this chapter the Navier stokes equations for a single liquid slug have been solved in order to predict the circulation patterns within the slug. Surface tension effects on the air water interface have been investigated by solving the Young Laplace equation. The calculated interface shape has been utilized to define the liquid slug geometry at the fr ont and tail interfaces of the slug. Then the effects of the surface tension on the hydrodynamics of the twophase slug flow have been compared to those where no surface tension forces exist. The importance of the complex flow field features in the vicinit y of the two interfaces has been investigated by defining a nondimensional form of the wall shear stress. The latter quantity has been formulated based on nondimensional parameters in order to define a general Moody friction factor for typical two phase slug flow in microchannels. Moreover, the hydrodynamics of slug flow formation has been examined using computational fluid dynamics (CFD). The volumeof fluid (VOF) method has been applied to monitor the growth of the instability at the air water interface The lengths of the slugs have been correlated to the pressure fluctuations in the mixing region of the air and water streams at an axisymmetric T junction. The main frequencies of the pressure fluctuations have been investigated using the Fast Fourier Tr ansform (FFT) method. Key words : Microchannels, Surface Tension, CFD, Slug Flows, Volume of Fluids (VOF) PAGE 31 31 Numerical Simulation Single Phase Simulations In this section we will calculate the flow field within a single liquid slug. The simplest geometr y for one liquid slug in a microchannel may be assumed as a rectangle. Figure 3 1 shows an experimental image of a typical slug flow in a microchannel. Figure 3 1. Experimental image of a typical slug flow in a microchannel (Serizawa et al 2002) If the gravitational acceleration is neglected, an axisymmetric boundary condition can be applied to the centerline of the microchannel. Slugs are stationary in the frame of reference that moves with the slug velocity, us, while the pipe (microchannel) wall is moving with the same velocity in the opposite direction. The interfaces at the front and tail of the liquid slug can be modeled as a zerostress boundary condition if the air viscosity is neglected compar ed to that of the liquid. For simplicity purposes (in this section) it is assumed that both interfaces are planar and are perpendicular to the pipe wall. The length of the slugs, Ls, can be significantly varied depending on the flow conditions. Figure 3 2 shows the calculated flow field for two different slug aspect ratios, Ls/D, in a twomillimeter diameter pipe This figure shows that the effect of the interfaces on the wall shear stress is dominant for the short slugs while it becomes less important PAGE 32 32 as s lugs lengthen. The center of flow circulation within the liquid slug moves closer to the tail of the slug as it increases in length. Figure 3 3 shows a nondimensional form of the wall shear stress over the length of the liquid slug employing the fully de veloped laminar shear stress, l for a straight pipe. The nondimensionalized wall shear stress can thus be expressed as 28 Rem w l wu where is the liquid density and um is the mean liquid slug velocity. The streamwise distance x has been non dimensionalized by the slug length in Figure 3 3 Figure 33 shows that the wall shear stress is much greater than the laminar fully developed shear st ress in the short slugs while for the long slugs the nondimensional wall shear stress value approaches unity. The large average shear stress in the short slugs is due to the local boundary layer growth near the corner of the tail and front interfaces. The shear stress for all cases approaches infinity at the leading and trailing edges of the slugs. It is interesting to compare the entrance length of the liquid slug to the general entrance length of a laminar pipe flow. For the latter, the hydrodynamic e ntry length may be obtained from an expression of the form (xe/D) =0.05Re (Langha r and Lafayette 1942). However, numerical simulations predict (xe/D) =0.0005Re for slug flows, which is one hundred times smaller than the value, obtained from the expression of Langha r and Lafayette (1942). The primary difference is due to the condition of the air water interfaces, which force the flow to change direction (typically by 90 degrees). The change in direction indicates the presence of a large radial velocity in th e vicinity of the front and tail of the slug. That radial velocity provides fast growth or decay for the boundary layer. PAGE 33 33 Figure 3 2. Axial velocity contours and streamlines for a 2mm pipe flow for two liquid slug aspect ratios: Top is 1 mm liquid s lug length and bottom is 5mm liquid slug length For intermediate and long slugs, Figure 33 shows that the majority of the liquid slug length would experience a fully developed condition in which the nondimensional shear stress is close to unity. A minimu m wall shear stress exists in all cases shown in Figure 3, due to the flow circulation within the slug. The location of these minima can be correlated to the flow field circulation within the slug (see Figure 3 2). As the center of PAGE 34 34 the circulation gets clo ser to the tail of the slug, the resulting minimum value of the wall shear stress becomes smaller and closer to the front of the slug. The minimum indicates that the axial velocity gradient is smaller than that in general laminar Poiseuille flow. Figur e 3 3. Nondimensional wall shear stresses over the liquid slug lengths and the averaged shear stresses for a 2mm diameter pipe, Re=2000 (slug moves to the left) The values of the minimum shear stresses in Fig ure 3 3 are more interesting when compared to those obtained by the analytical model of Black et al (2007). The stresses PAGE 35 35 calculated by their model have the same or larger values than those obtained by the general laminar Poiseuille flow model. Since Black et al (2007) assumed that Re<1 in order to neglect the inertial terms in the Navier Stokes equations, it could be concluded that the minima shown in Figure 3 3 are due to the inertial terms in the Navier Stokes equations. The normalized averaged wall shear stress can be defined as 1 0 ) / ( *) / ( ) / (sl x s l w avel x d ( 3 1) The above parameter, ave represents the average nondimensional value of the wall shear stress over the slug length. Figure 33 shows a strong dependency of the above parameter on the slug length. Analytical Model of Gas Liquid Interface The surface tension force creates curvature at the air water interfaces. Therefore, the rectangular domain chosen in the above analysis is merely a simplifying assumption. For a more realistic cal culation, the interface shape needs to be taken into account. For particular conditions, the interface curvature can be evaluated from the experimental images. However, since experimental images are not available for all different microchannel diameters and slug velocities, using an analytical model in order to define the interface shape is necessary. The Young Laplace equation is a nonlinear partial differential equation that describes the capillary pressure difference sustained across the interface betwee n two static fluids due to the surface tension force. It relates the pressure difference to the shape of the surface and it is fundamentally important in the study of capillary surfaces. If the effect of the gravitational acceleration is neglected, PAGE 36 36 the int erface would be axisymmetric. The YoungLaplace equation for an axisymmetric interface, Z(r), can be expressed as ) ) 1 ( ) 1 ( (2 / 1 2 2 / 3 2Z r Z Z Z p ( 3 2) radius and the pressure difference, can be obtained either from CFD simulations or from analytical models. Equation ( 3 2) is a secondorder ODE with the following boundary conditions: 0 0 r for dr dZ ( 3 3) a r for dr dZ ) 2 / tan( (3 4) The finite difference method, applying NewtonRaphson iteration, was utilized to solve Equation ( 3 defined from the experimental data. The calculated interface shape has been used to define a more realistic geometry of the liquid slug for the CFD simulation. Figure 3 4 shows a sample interface shape for an air water system while the contact angle is The flow field does not noticeably change if the more realistic interfa ce shape is included in the CFD simulation. However, including the interface shape affects the wall shear stress. PAGE 37 37 Figure 3 4. Calculated air water interface profile, Z(r) for a 2 mm diameter microchannel, and Introducing Moody Friction Factor for General Two Phase Slug Flow Figure 3 5, in which the interface curvature is included in the simulation, shows the nondimensional wall shear stress for the same slug lengths that are reported in Figure 3 3 Comparisons between Figures 3 3 and 3 5 reveal that the predicted wall shear stresses for the slugs accounting for the interface curvature are always smaller than those if the interface were assumed as a straight surface. A large difference in the predicted wall shear stress can be observed for the short slugs since the interface shape has a local effect on the flow field (47% for the slug length of Ls/D=0.5 while for the long slugs Ls/D>7.5 the differences are less than 8%). The main difference between Figures 3 3 and 3 5 is due to the consi derable difference in behavior of the wall shear stress at the corner of the slug. PAGE 38 38 Figure 3 5. Axial velocity contours and streamlines (top). Nondimensional wall shear stress over slug length and the average shear stress for a 2mm diameter pipe (bottom), Re=2000 (Interface shape is included in the CFD simulation) The infinity value of the wall shear stress at the corner in the case of the curved interface drops faster than that when the interface is assumed as a straight surface. The averaged non dimensional wall shear stress, ave [ Eq uation (3 1 )] for different nondimensional slug lengths as a function of the Reynolds numbers is shown in Figure PAGE 39 39 3 6. This figure shows that generally the average non dimensional wall shear stress increases as the Reynolds number increases. The flow is assumed laminar in all of the studied cases in Figure 3 6 since the Reynolds numbers are smaller than the critical Reynolds number for a typical pipe flow ( Recritical tical Reynolds number is not quite well defined in slug flows for different slug lengths (Mehdizadeh et al 2009). Figure 3 6. Averaged nondimensional shear stress for different slug lengths and Reynolds numbers According to Figure 36, the dependency of the average wall shear stress, ave on the Reynolds number can be expressed as c aveb a Re* for Re < 2100 ( 3 5) PAGE 40 40 where the coefficients a, b and c are defined in Table 31. The complexity o f the coefficients is due to the complexity of the flow field close to the front and tail interfaces. The average error in Equation ( 3 5) due to interpolation is less than 7%. Figure 36 or alternatively Eq. ( 3 5) can be used to calculate the pressure drop of slug flows in microchannels for different slug conditions. Therefore the conventional laminar Moody (or Darcy) friction factor, f, for general twophase slug flow can be introduced as follows: ) 1 )( Re Re ( 64) 1 ( cb a f (3 6) where the (1multiplier indicates the portion of pipe length in contact with the liquid slugs. While numerical simulation coupled with interpolation produce the above form of the Moody friction factor, experimental studies are needed in order to modify or verify the above coefficients. Table 3 1. Coefficients in e quation ( 3 5) a b c 00678 0 08 1) ( 04 1 ) ( 971 0 D L D Ls s 0668 0 0653 0) ( ) ( 365 0 D L Exp D Ls s 11 0 34 3) ( 189 0 ) ( 0344 0 D L Exp D Ls s First Attempt for Two Phase Simulations The focus of this section is on the hydrodynamics of slug flow formation. Al though, there are many numerical modeling studies on different two phase flow regimes, slug flows are one of the least studied regimes primarily due to modeling complexities. De Schepper et al. (2008) compared the performance of the existing numerical tools and approaches for modeling of twophase flows. Their qualitative comparisons with the PAGE 41 41 experimental images showed that, except for the slug flow regime, the numerical simulation could more or less capture two phase flow regimes. They examined the effect of several tuning parameters in their simulations to overcome the failure of capturing the slug flow regime, but they had no success. They presumed that the failure was caused by the relatively small region of the slug flow regime in the Baker chart (Baker 1954). The small region may have caused an unpredictable transition among regimes which could have affected the simulations. The geometry and the mesh used for the numerical simulations are shown in Figure 3 7. This geometry is the same as that of the test rig employed by Abdelall et al. (2005) who studied the pressure drop caused by abrupt flow area changes in a microchannel. Mehdi zadeh et al. (2009) analytically examined the same problem with the same geometry and proposed a new approach for computing the pressure drop and gas volume fraction in the vicinity of an abrupt contraction. Figure 3 7 shows air and water being introduce d into an axisymmetric T junction with specific constant mass flow rates. The effect of gravitational acceleration is neglected, thus the flow can be assumed axisymmetric. This fact helps convert a 3D geometry into a 2D axisymmetric geometry. Kashid et al (2007) successfully simulated liquid liquid slug flows in a microchannel. They reported that the interfa ce may undergo severe deformations and break ups in flows with large viscosity and density differences (e.g. air water flow). This makes simulations more complicated The Eulerian approach, which formulates interface tracking by two methods: front trackin g and volume tracking, has been used. The volume of fluid (VOF) method was developed by Hirt and Nichols (1981) In this method, two or more phases are not PAGE 42 42 allowed to interpenetrate each other. For each additional phase the volume fraction of the phase can be added in the computational cell. The summation of the volume fractions of all phases in each control volume is set to unit y. Figure 3 7. Axisymmetric computational domain (top half is used in numerical simulations, while the lower half is a mirror image of the top half) The VOF method is probably one of the most convenient methods relative to other two phase flow modeling m ethods. The method has reasonable accuracy, is less computationally intensive (compared to the Eulerian model), is relatively simple, and can solve highly complex free surface flows. In this chapter the liquid phase is assumed PAGE 43 43 as an incompressible fluid wh ile the gas phase is assumed to behave like an ideal gas. The VOF method is built to solve a single momentum equation throughout the domain, while the resulting velocity field is shared among the phases. The momentum equation, shown below, is dependent on the volume fractions of all phases through the properties and SF TF g v v p v v v t )] ( .[ ) .( ) ( (3 7) w here SFF and in the above equation can be calculated from the following eq uations: g g l l (3 8) and g g l l g g g l l l (3 9) and where l and g are the volume fractions of the liquid and gas phases, respectively. The term SFF in Eq uation ( 3 7) represents surface tension and wall adhesion effects on the flow filed. The surface tension is due to the strong intermolecular attractive forces among molecules in a fluid. These forces contribute to minimizing the surface area of slug. The wall adhesion mentioned earlier in the chapter is due to the stronger attractive forces among liquid molecules relative to their attraction to the wall. This causes the fluid to create some contact angle with the wall. The continuum surface for ce (CSF) model proposed by Brackbill et al (1992) is used to PAGE 44 44 model surface tension. The surface curvature is computed from the local gradients in the surface normal at the interface. Thus, the source term, SFF in the momentum equation can be specified as follows (see Brackbill et al 1992): ] ) ( 5 0 [ g l g g l l SFn k F (3 10) where n is the surface normal vector and k is the surface curvature ln (3 11) )] ( ) [( 1 ) ( n n n n n n k (3 12) A wall adhesion angle in conjunctio n with the surface tension model is also included in this simulation. The model is again taken from the work of Brackbill et al (1992) The contact angle that the fluid makes with the wall is used to adjust the surface normal in cells near the wall. This so called dynamic boundary condition results in the adjustment of the curvature of the surface near the wall (the same approach used in generating Figure 3 4). If is the contact angle at the wall, then the surface normal at the live cell next to the wall is sin cos w wt n n (3 13) w here wn and wt are the unit vectors normal and tangential to the wall, respectively. The tracking of the interfaces between the phases is achieved through solving the continuity equation for the volume fraction of the second phase as follows (see Brackbill et al 1992): PAGE 45 45 )] ( ) .( ) ( )[ 1 (lg .m m v tgl l l l l l l (3 14) The term glm. is the mass transfer rate from the gas phase to liquid phase and lg .m is the mass transfer rate fr om liquid phase to the gas phase. T he mass transfer between gas and liquid phases has been neglected. The volume fraction equation [ Eq uation ( 3 14) ] will not be solved for the primary phase; the primary phase volume fraction will be computed based on the f ollowing constraint: 1 g l (3 15) The volume fraction equation is solved through explicit time discretization. The geometric reconstruction scheme (see Youngs 1982) can be used to represent the interface between fluids using a piecewiselinear approach. This approach assumes that the interface between two fluids has a linear slope within each cell, and uses this linear shape for calculation of the advection of fluid through the cell faces (Youngs 1982). The above equat ions were solved using the finite volume method which is employed in the commercial fluid flow solver FLUENT 6.3. The geometry of the two phase simulations is shown in Figure 3 7. A structured rectangular grid was used for the simulation. The grids were re fined in the vicinity of the axisymmetric T junction in order to capture the hydrodynamics of the slugs generation and break up. The number of elements was selected as 28,057 which is 10 times larger than what is needed for a single phase simulation. In t he next chapters we will show that not only the number of meshes is important but also the aspect ratios of elements significantly affect on the PAGE 46 46 obtained results. For all discretizations, the thirdorder QUICK scheme has been used. A pressure boundary condition was used at the outlet of the channel with a value equal to the atmospheric pressure. The time steps were fixed at 2106 seconds. For this time step the maximum Courant number at the domain was found to be less than 0.42. In order to simulate liquid liquid two phase slug flow, Kashid et al (2007) used time steps that are fifteen times larger than the on es employed in this investigation. Since interface deformations and fluctuations for air water two phase flows are much faster than those in liquidliquid twophase slug flow, very small time steps are necessary. Simulations were performed on a 2.61 GHz mi croprocessor. Unsteady simulation took almost 8 weeks for each test case (for 20 seconds of flow time). Slug s Formation, Pressure Fluctuations and Main Frequencies Figure 38 shows a predicted slug flow regime in a microchannel in which many liquid slug lengths can be observed. This figure shows that instability in the air water interface (vortex shedding) produces short and long waves close to the axisymmetric T junction. Small waves grow shortly downstream of the T junction and the accumulated water i n the waves creates liquid slugs. A closer look at the instabilities in the vicinity of the junction reveals that when air and water streams cross paths, several short and long vortices form at the air liquid interface. Kuru et al. (1995) used linear stability theory in order to study transition from the separated flow regime to the slug flow regime. Comparison of the results obtained by Kuru et al. (1995) with experimental images reported in their work indicates that the approximate models do not correctly predict the point of neutral stability and that models differ substantially from the predictions of the full Navier Sto kes equations. They have concluded that the non linear terms that they have neglected play a significant role in instability behavior when the waves are PAGE 47 47 growing, thus a linear instability theory is inherently incapable of predicting transition from the sep arated flow regime to the slug flow regime. Figure 3 8. Predicted slug flow regime in a microchannel for a water air system In a more sophisticated study of the slug formation phenomenon, the pressure fluctuations and the generated wave lengths (slug len gths) have been monitored by solving the full Navier Stokes equations. Figure 3 9 (left) shows the time sequences of slug formation. This figure shows one period in which one slug forms close to air water mixing region. In Fig ure 3 9 (left) it can be seen that air penetrates into the next water pocket. As air penetrates, a small wave forms at the horizontal air water interface. The formed wave grows and creates a neck between two consecutive air slugs. The air neck becomes progressively thinner due to the growth of the water wave at the air water interface. When the height of the wave reaches the radius of the channel, water molecules in both sides pull each other and break the air neck and produce an air slug. The wave length of the generated wave is the s ame as the length of the formed air slug. PAGE 48 48 The generated air and water slugs change their shape and attach to the channel wall shortly downstream of the T junction due to the wall adhesion mentioned earlier (not shown here). Figure 3 9 (right) shows pressure frequencies in the T junction in both the time and frequency domains. Since it is difficult to find the main frequencies by merely looking at the pressure fluctuations in the time domain, the Fast Fourier Transform (FFT) method has been applied to the pressure signals. The pressure frequencies in the frequency domain have been calculated by the FFT algorithm using MATLAB. According to Fig ure 3 9 (right) the most amplified spectrum can be split into two regions. The highfrequency region, f >500 Hz, corr esponds to numerical solution noises that have no physical significance. The low frequency region, 0< f <500 Hz, on the other hand, corresponds to fluctuations at the interfaces and slug formations, both of which produce pressure waves that propagate into th e channel. In Figure 3 9 (right) the most amplified frequencies are found in the following frequencies: 0.8, 1.5, 4, 8.4, 12 and 20 Hz. These frequencies correspond to the slug formation phenomenon. Each of the above frequencies is associated with a speci al length of the liquid slug. Specifically, lower frequencies correspond to longer slugs, while higher frequencies correspond to shorter slugs. PAGE 49 49 Figure 3 9. Sequence of slug formation. (left) Time sequence of slug formation close to the T junction (su perficial velocities (right) Pressure fluctuations in the mixing region ( right top); Pressure fluctuations amplitude spectrum of pressure ( right bottom); Number of data points 65,536. (Superficial velocities, Jw=18.75 mm/s and Jg=7.48 mm/s) Conclusions In this chapter, the Navier Stokes equations within a single rectangular liquid slug have been solved. Surface tension effects on the air water interface were analytically investigated employing a secondorder Young Laplace differential equation. The calculated interface shape was then used to define a more realistic geometry and a set Neck PAGE 50 50 of boundary conditions for numerical simulation purposes. A large difference of the predicted average wall shear stress relative to laminar wall shear stress was observed for t he short slugs since the interface shape had a local effect on the flow field (47% for the slug length of Ls/D =0.5 while for the long slugs Ls/D >7.5 the differences were less than 8%). A correlation capable of predicting the Moody friction factor has been introduced. Also, formation of slugs through the microchannel along with their merger and break up has been numerically investigated using the two phase volume of fluid (VOF) method. Results demonstrated that the shape of the interface through which two ph ases are mixing has a significant effect on the geometrical properties of the slug. Finally, the frequencies of the pressure fluctuations (generally difficult to measure experimentally in microchannels) in the mixing region were reported. Using Fast Fourier Transform (FFT), the main frequencies of the pressure fluctuations were extracted in order to correlate them to the slug lengths. PAGE 51 51 CHAPTER 4 CFD MODELING OF TWO PHASE SLUG FLOW IN MICROCHANNELS Outlook Despite of the fact that numerical simulation of t wo phase flows in microchannels has been attempted by many investigators, most efforts seem to have failed in correctly capturing the flow physics, especially the slug flow regime characteristics. The presence of a thin liquid film in the order of 10 m ar ound the bubble (sometimes called gas pocket or gas slug) may be a contributing factor to the above difficulty. Typically, liquid films have a significant effect on the flow field. Thus, there is a strong motivation to employ numerical simulation methods i n order to avoid some of the experimental difficulties. In this chapter, the characteristics of two phase slug flow in microchannels are calculated with the help of the Volumeof Fluid (VOF) method. Formation of the slugs for different superficial velociti es, Capillary numbers, and gas volume fractions are investigated. The minimum mesh resolution required to capture the liquid film surrounding the gas bubble is reported employing a dynamic mesh adaption methodology with interface tracking. Results are show n to be in good agreement with experimental data and empirical correlations. Key words : Microchannels, Capillary number, CFD, Slug Flows, Volume of Fluids (VOF) Method Numerical Simulation Background : Slug flows are one of the least numerically studied reg imes in twophase flows in large part due to modeling complexities. De Schepper et al. (2008) compared the performance of the existing numer ical tools and approaches for modeling of two phase flows. Their qualitative comparisons of numerical results with experimental PAGE 52 52 flow visualization data showed that, except for the slug flow regime, numerical simulations could more or less capture twophase flow regime physics. They examined the effect of several tuning parameters in their simulations to overcome the failure of capturing the slug flow regime, but they had no success. They presumed that the failure was caused by the relatively small region of the slug flow regime in the Baker chart (Baker 1954). The small region may have caused an unpredictable transition among regimes which could have affected the simulations. Kashid et al (2007) calculated the mass transfer between two consecutive slugs in the liquidliquid slug flow microreactor utilizing the CFD simulation. They have simplified their problem to two rectangular slugs by means of neglecting the surface tension effect on the shape of the slugs. Their mass transfer prediction showed a fair agr eement with experimental data at high slug velocities but poor results at low velocities. The multiphase flow pattern in microchannels is extensively studied for different two phase fluids, Capillary numbers and pipe diameters. Gas liquid twophase flow p atterns in microchannels have been studied in detail by Serizawa et al. (2002) and Triplett et al. (1999). Figure 4 1 shows a typical gas liquid two phase flow pattern in microchannels. The absence of the stratified flow regime in Figure 4 1 should be noted. Problem Formulation and Solution Procedure Governing Equations The Volumeof Fluid (VOF) method is employed here to simulate gas liquid twophase slug flows in microchannels. The VOF method was developed by Hirt and Nichols (1981) In this method, two or more phases are not allowed to interpenetrate each other. PAGE 53 53 Figure 4 1 (Repeated). Gas liquid multiphase flow pattern in a microchannel, Serizawa et al. (2002) on 20m ID pipe (top); Triplett et al. (1999) on 1.097 mm ID pipe (bottom) The summation of the volume fractions of all phases in each control volume is set to unity. The VOF method is probably one of the m ost convenient methods relative to other twophase flow modeling methods. The method has reasonable accuracy, is less computationally intensive (compared to the Eulerian model), is relatively simple to use, and can solve highly complex free surface flows. B oth the liquid and gas phase are assumed as incompressible. The VOF method is structured to solve a single momentum equation throughout the domain, while the resulting velocity field is shared among the phases. The momentum equation, shown below, is dependent on the volume fractions of all phases through the properties and PAGE 54 54 SF TF g v v p v v v t )] ( .[ ) .( ) ( ( 4 1) w here SFF and in the above equation can be calculated from the following equations: G G L L ( 4 2) and G G L L G G G L L L ( 4 3) w here L and G are the volume fractions of the liquid and gas phases, respectively. The term SFFin Eq. ( 4 1) represents surface tension and wall adhesion effects on the flow field. The surface tension is due to the strong intermolecular attractive forces among molecules in a fluid. These forc es contribute to minimizing the surface area of the slug. The wall adhesion is due to the stronger attractive forces among liquid molecules relative to their attraction to the wall. This causes the fluid to create some contact angle with the wall. In this chapter the Continuum Surface Force (CSF) model proposed by Brackbill et al. (1992) is used to model surface tension. The surface curvature is computed from the local gradients in the surface normal at the interface. Thus, the source term, SFF, in the momentum equation can be specified as follows (see Brackbill et al. 1992): ] ) ( 5 0[ G L G G L L SFn F ( 4 4) w here n is the surface normal vector and is the surface curvature PAGE 55 55 Ln ( 4 5) )] ( ) [( 1 ) ( n n n n n n (4 6) A wall adhesion angle in conjunction with the surface tension model is also included in this simulation. The model is again taken from the work of Brackbill et al. (1992) The contact angle that the fluid makes with the wall boundary is used to adjust the surface normal in cells near the wall. This so called dynamic bounda ry condition results in the adjustment of the curvature of the surface near the wall. If is the contact angle at the wall, then the surface normal at the live cell next to the wall is sin cos w wt n n (4 7) w here wn and wt are the unit vectors normal and tangential to the wall, respectively. Some researchers believe that the wall adhesion properties, such as wall contact angle, play an important role in slug flow simulat ion. However, this may not be entirely true (Gupta et al. 2009). The role of the wall contact angle is only significant when both of the phases are in contact with the solid wall or the liquid film around the gas bubble is so thin that the van der Waal for ces can act across the film. This may be the case for T or Y type junctions where both phases remain in contact with the wall but once a slug flow is formed with a liquid film at the wall, wall adhesion plays no role in determining the bubble shape. Moreov er, the modeling of wall adhesion requires inclusion of the contact angle hysteresis effect to obtain meaningful results (Gupta et al. 2009). The tracking of the interfaces between the phases is achieved through solving the continuity equation for the volu me fraction of the second phase as follows (see Brackbill et al. 1992): PAGE 56 56 ) ( ) .( ) (. LG GL L L L L Lm m v t ( 4 8) The term GLm. is the mass transfer from the gas phase to the liquid phase and LGm. is th e mass transfer from the liquid phase to the gas phase. T he mass transfer between gas and liquid phases has been neglected. The volume fraction equation, Eq. ( 4 8), will not be solved for the primary phase; the primary phase volume fraction will be computed based on the following constraint: 1 G L ( 4 9) Differencing Schemes The above equations were solved using the finite volume method which is employed in the commercial fluid flow solver FLUENT 6.3 (Fluent 6.3.26 Users' Guide, 2006). For momentum discretization, the third order QUICK scheme has been used. The Pressure Implicit with Splitting of Operators (PISO) scheme is used for pressure velocity coupling (Issa 1986). For VOF multi phase model, upwind schemes are generally unsuitable for interface tracking because of their numerical diffusive nature. Also the central differencing schemes, while generally able to maintain the sharpness of the interface, are boundless and usually give unphysical results (Darvish and Moukalled 2006; Muzaferija et al. 1998). In order to avoid these problems, a modified version of the High Resolution Interface Capturing (HRIC) scheme was used (Muzaferija et al. 1998). The modified HRIC scheme consists of a nonlinear mix of upwind and downwind differencing. The modified HRIC scheme provides improved accuracy for VOF calculations when compared to QUICK and secondorder schemes and is less PAGE 57 57 computationally costly than the GeoReconstruct scheme (Youngs 1982, F luent 6.3.26 Users' Guide, 2006). A first order, noniterative fractional step scheme is used for the time marching of the momentum and continuity equations. The time step in Fluent is controlled by a specified maximum value for the Courant number. The Co urant number ( Co ) is defined as V x t Co / (4 10) w here x is the grid size and V is the fluid velocity. A maximum Courant number of 0.25 was set in the present simulation and a variable time step based on the set Courant number was used. Since the elements are very fine close to the interface, the time steps are typically in the order of 107 s. However, at the instance of the slugs break up, the time steps decrease to the order of 108 s because of fast deformations of the interface. Interface Tracking The VOF method allows tracking of immiscible interfaces with surface tension effects. Kashid et al. (2007) successfully simulated liquid liquid slug flows in a mic rochannel. They reported that the interface may undergo severe deformations and break ups in flows with large viscosity and density differences (e.g. air water flows). Severe interface deformations make simulations more complicated since all the high frequ ency fluctuations at the interface have to be correctly tracked. Therefore, an ultra fine mesh and ultra small time steps are necessary for proper gas liquid two phase flow simulations (order of 107 s). PAGE 58 58 In this investigation, structured squared shapes of grids were used for all simulations. It was observed that grids whose aspect ratios are different from unity produce nonphysical results as can be seen in Figure 4 3a. This fact has also been reported by Gupta et al (2009). A uniform base grid size of D/ 30 (53 and 16 m in different studied test cases) was used in all simulations. Even though this grid size is very fine for common single phase studies, it is not enough for capturing the severe deformations of interfaces and the breaking up of slugs common ly found in these flows. A good overview of the effect of the grid size on slug flow simulations was reported by Gupta et al (2009). They tried different grids having element size ranging from D/20 to D/100. However, physical results were not obtained unt il they refined the mesh size to D/200 in the vicinity of the wall. Due to the ultra fine grid, their simulations were extremely costly. In order to reduce the computational time, in this study reported in this chapter, a dynamic mesh adaption method with interface tracking has been employed. It is to be noted that the gradient of the volume fraction at the interface is in the order of 0.4 to 0.6. This indicates that the thickness of the interface in the computational domain is about 3 (or 1/0.4) elements w ide. Since mathematically the thickness of interface is zero, 3 elements wide should be significantly smaller than the other length scales in the computational domain. In the present work, the base mesh (D/30) was refined close to the interface using coars ening and refining thresholds for the gradient of the volume fraction. The base mesh is refined 8 folds when the gradient of the volume fraction in the vicinity of interface exceeds 0.1. The mesh regenerates every 5 time steps in order to refine elements i n the vicinity of the new location of the interface. Cells with volume fraction gradients below the 0.15 threshold are marked for coarsening. Using this PAGE 59 59 method we have repeated Gupta et al .s (2009) test case and were able to achieve better results employi ng 96% less mesh than their simulations. More on this topic will be presented later in the chapter. Model Geometries and Studied Test Cases Schematic geometries of the test cases investigated are shown in Figure 4 2. Two inlet configurations have been exam ined. In the first, the water stream enters the channel through an axisymmetric T junction, while in the second, two parallel streams enter the channel with an inlet void fraction equal to the homogenous void fraction (see Figure 4 2 bottom). Figure 4 2 Schematic showing the test cases and boundary conditions (not to scale) Figure 42 shows air and water being introduced into the microchannels with specific constant mass flow rates. Water and air densities are 997 and 1.18 kg/m3 while viscosities are 8 .89x104 and 1.83x105 kg/(m.s), respectively. The effect of gravitational acceleration is neglected, thus the flow can be assumed axisymmetric. This fact helps convert a 3D geometry into a 2D axisymmetric geometry. We have tried to cover a wide range of flow parameters and Capillary numbers. PAGE 60 60 Details of the test cases are reported in Table 4 1. Runs 1 to 8 pertain to the geometry shown in the top half of Figure 4 2, while Runs 9 to 12 pertain to that in the lower half. Run 12 is identical to the simulated test case of Gupta et al (2009). The results will be compared to the experimental data and Gupta et al .s (2009) results. Table 4 1. Test c ases e xamined Results and Discussion Bubble Shape For testing the mesh sensitivit y, five different mesh resolutions were tried. The first was a rectangular mesh with an aspect ratio of 5 and a smaller side equal to D/30. Results of this mesh were very poor as can be seen in Figure 43a. Gupta et al. (2009) reported similarly poor resul ts for large aspect ratios. In the second employed mesh, the base mesh was a structured square grid with a size of D/30. In the third mesh, elements close to the interface were refined 4 folds (22). Predicted results of mesh numbers 2 and 3 were not physi cal since they could not capture the thin liquid film around the bubble at all. Figure 43 displays sample results for mesh numbers 2 and 3, which show a dry wall when gas slug passes. Figure 43b is a typical flow configuration commonly Run Number Pipe diameter, D (mm) .) / ( s g mL .) / ( s g mG JL(m/s) JG(m/s) Capillary number ) (G L lJ J Ca Homogenou s Gas Void Surface Tension 1 1.60 1.152 0.00534 0.57 2.242 0.0385 0.795 0.079 2 1.60 1.152 0.00641 0.57 2.694 0.0150 0.824 0.217 3 1.60 1.152 0.00108 0.57 0.455 0.0149 0.442 0.068 4 1.60 4.995 0.01046 2.49 4.393 0.0150 0 .638 0.458 5 1.60 20.079 0.00701 10.01 2.942 0.0150 0.227 0.86 6 1.60 6.293 0.00534 3.13 2.242 0.000074 0.416 73.0 7 1.60 0.427 0.00036 0.21 0.154 0.00503 0.419 0.073 8 1.60 1.218 0.00118 0.60 0.498 0.0151 0.450 0.073 9 1.10 0.201 0.00017 0.21 0.154 0 .0050 0.419 0.073 10 1.10 0.576 0.00056 0.60 0.498 0.0152 0.450 0.073 11 1.00 0.199 0.00022 0.255 0.245 0.0069 0.49 0.073 12 0.50 0.049 0.00005 0.255 0.245 0.0069 0.49 0.072 PAGE 61 61 used for simplifie d models. Mehdizadeh et al. (2009b) and Kashid et al (2007) solved the Navier Stokes (N S) equations for the axisymmetric rectangular liquid slug for a frame of reference moving with the liquid slug. Results of both investigations are similar to those obt ained from this investigation as displayed in Figure 4 3b. The mesh in the vicinity of the interface is refined 8 (i.e. 23) and 16 (i.e. 24) folds for mesh numbers 4 and 5, respectively, as defined earlier. Both of these meshes could precisely capture t he thin liquid film around the bubble. Figure 4 4 shows a comparison between the predicted results and those obtained via experimental flow visualization as reported by Triplett et al. (1999). (a) (b) Figure 4 3 Volume fraction contours for the cour se meshes. Mesh No. 1, mesh size D/30D/6 (a) and mesh No. 2, mesh size D/30D/30 (b). Figure 4 4a shows that the shape of the slugs is very similar to the experimental image of Triplett et al (1999). It is seen that the calculated thickness of the inter face is very thin since the elements in the vicinity of the gas liquid interface is refined to 4.1 m PAGE 62 62 (D/270). This refined grid can capture the liquid film around the bubble. In Run 9 (Figure 4 4a) since the air and water superficial velocities are low, t he interface is pretty smooth with very small fluctuations at the nose. Also, the tail of the gas slug has a smooth convex shape. The small fluctuations are produced by the propagation of pressure waves created by the slugs breakup at the inlet. Thus, far ther away from the inlet, those waves cause the gas liquid interface to vibrate slightly. For higher superficial velocities, the shape of the gas slug becomes bullet like. Also, generation of some waves at the interface can be noticed in Figure 4 4b. The bullet like shape of the gas slug indicates a higher pressure drop in the channel when compared to the previous case (Figure 4 4a). The high pressure in the back of the gas slug pushes the interface in the streamwise direction, thus causing the back side o f the gas slug to become flat. For higher superficial velocities, the pressure drop is even higher and pushes the back of the gas slug harder. This causes the interface to have a concave shape. After reaching a threshold for superficial velocities (or pres sure drop), the concave shape stretches into the gas slug and helps create annular or slug annular regimes. In Figure 4 4b, it can be seen that the numerical simulation is able to successfully capture interface instabilities. For higher superficial velocit ies, the generated waves at the interface grow both temporally and spatially. This may lead to a two phase flow regime transition. It should be noted that the experimental images of Figure 4 4 have been taken vertically. Thus, the thin liquid film is sligh tly thicker at the bottom of the gas slug and thinner at the top. However, gravity has been neglected in the numerical simulations. Thus, the liquid film is perfectly axisymmetric in the predicted results. If the experimental image were taken PAGE 63 63 from the top of the test section, it would have been better for comparison since gravity has the least effect from the top view. Figure 4 5 shows a comparison between the present simulations (Run 12) and those of Gupta et al. (2009), which is an identical test case. (a) (b) Figure 4 4. Comparison of volume fraction contours with the experimental images of Triplett et al. (1999) on a 0.5 mm pipe with refined mesh in the vicinity of the int erface. Run 9 (Figure 4 4a) and Run 10 (Figure 4 4b). In the present study approximately 16,000 elements were used, while Gupta et al. (2009) employed 448,000 elements. It can be seen that both simulations could capture the thin liquid film. The aforementioned fact is significant since the twophase slug flow PAGE 64 64 simulation is inherently very costly. The presence of the liquid film (11m thickness) and the mesh adaption close to the interface can be seen in Figure 4 5. Velocity a nd Pressure Field Figure 46 shows static pressure contours and velocity vectors in a frame of reference moving with the gas slug. In the figure, one full gas slug is shown upstream of a portion of the next gas slug. The interfaces of both gas slugs are shown with bold solid lines. Ma ny flow circulation patterns can be seen inside of gas slugs, while one circulation pattern is shown inside the liquid slug. As can be seen, high velocities exist at several interfacial locations due to fast interfacial vibrations which are generated by pr essure waves created by the slugs breaking up in the upstream section of the c hannel. The flow field in the liquid slug is different than that predicted in the simplified models of Mehdizadeh et al. (2009b) and Kashid et al. (2007). However, the number of eddies and the location of their centroids inside the liquid slug are similar to those predicted by the simplified models. The velocity inside of the liquid film is found to be almost linear, while nonlinearities can be observed in the vicinity of the mini mum clearance between the gas slug and the wall. The pressure inside gas slugs is found to be approximately uniform. Figure 4 6 shows that the pressure of gas slugs downstream is smaller than that upstream, mainly due to pressure drop. Also, it can be seen that the pressure inside a gas slug, is higher than the average pressure in a liquid slug along with that in the liquid film. This is due to the curvature of the interface and surface tension effects. PAGE 65 65 Figure 4 5. Simulation results for Run 12 (dynam ic refined mesh at the interface employing 15,639 elements ) [top]; Simulation results of Gupta et al. (2009) (employing 448,000 elements) [bottom] Figure 4 6. Velocity vectors and pressure contours in a frame of reference moving with the gas slug (Run 7) Volume Fraction The gas volume fraction is defined as x x S xL G L G ) 1 ( (4 11) PAGE 66 66 where S is the slip ratio between phases, x is the gas mass quality and L and G are liquid and gas phase densities, respectively. The homogenous void fraction, is defined if the slip ratio in the above equation becomes unity ( ) (L G GJ J J ). The Armand correlation (Armand and Treschev, 1946) is widely used to predict the gas volume fraction in microchannels as follows: G 833 0 (4 12) The void fraction is calculated by dividing the volume of the gas by that of the mixture wherever the gas slug forms. Employing integration over the length of two successive established slugs the calculated void fraction for Run 12 is found to be 0.4265, while the Armand correlation gives a value of 0.4081 (see Figure 4 5). Thus, the calculated value differs by 4% from that predicted by the Armand correlation. Gupta et al. (2009) reported a 10% error for the same test case. Gas Slug Velocity and Length Even though many local small fluctuations of the velocity and pressure were reported in Figure 4 6, the time averaged interface axial velocity (gas bubble velocity) is constant. The velocity of the gas slug UG based on the empirical correlation of Eq. ( 4 12), can be calculated as G GS G G GU Gx U (4 13) where G is the mass flux. A comparison between the calculated velocity of the gas slug and that obtained fr om the above correlation is shown in Table 4 2. PAGE 67 67 As can be seen, the average error in the reported test cases is less than 5%. This indicates that the macro properties of slug flow physics could be captured by numerical simulations. Table 4 2. Comparison of v elocity of g as s lugs for d ifferent t est c ases with the a vailable e mpirical c orrelation and p revious n umerical s imulations Description Experimental correlation (Eq. ( 4 13)) Gupta et al. s (2009) simulation Present work Run 7 0.440 N/A 0.423 Run 9 0.44 1 N/A 0.427 Run 12 0.6 0.55 0.552 The length of the slugs significantly depends on the inlet configuration. For Run 12, the length of the gas slug is 1.7D and the liquid slug length is 1.12D, defined at the centerline of the microchannel. Liquid Film Thickness In the Introduction, it was explained that Gupta et al. (2009) compared five different mesh resolutions and concluded that if the size of the elements were to be approximately of the liquid film thickness, the liquid film would have been succ essfully captured. However, achieving the above goal was extremely costly (450,000 elements for Gupta et al .s case). Using the dynamic mesh adaption method with interface tracking, the liquid film around the bubble has been successfully captured (see Figures 4 4, 45 and 46) with approximately 16,000 elements employed. This indicates that the number of elements can be reduced by 96% and the computational time can correspondingly be reduced by approximately 90% (including the time required for dynamic mesh adaption). Moreover, using a smaller number of elements in areas that do not need mesh refinement can potentially reduce numerical diffusion and errors. PAGE 68 68 Table 4 3 compares the calculated liquid film thicknesses with those reported by Gupta et al. (2009) and by several empirical correlations available in the literature. Table 4 3. Comparison of l iquid f ilm t hicknesses with those calculated by available empirical correlations and G upta et al .s numerical simulation Description (m) 2 15 0 Ca R (Fairbrother and Stubbs, 1935; Taylor, 1960) 10 3 23 1 Ca R (Bretherton, 1961) 11 R ) 34 1 ( 5 2 1 34 13 / 2 3 / 2Ca Ca (Aussillous and Quere, 2000) 10 Gupta et al. (2009) with 5 m near wall fine grid 15 Gupta et al. (2009) with 2.5 m grid 12 Gupta et al. (2009) with 1.67 m grid 12 Present work for Run no. 12 (mesh no. 4) 11 Present work for Run no. 12 (mesh no. 5) 11 As can be seen, Table 4 3 and Figure 4 5 prove that microscale properties of slug flows can also be captured accurately by the present work. Wall Shear Stress and Pressure Drop There exist many discrepancies in the calculated pressure drop of twophase flows in microchannels, mainly due to difficulties in measurem ents (Triplett et al. 1999, Abdelall et al. 2005, Kreutzer et al. 2005, Black et al. 2007, Hayashi et al. 2007, Mehdizadeh et al. 2009b, and Gupta et al. 2009). Figure 4 7 shows the pressure distribution on the axis and the wall of the microchannel for Run 12. As can be seen, the pressure inside of the gas slug is higher than that in the surrounding liquid. Pressure distribution is roughly uniform inside of any gas slug because momentum of gas is negligible in compared to that in the liquid slug. On the other hand, pressure varies significantly due to large momentum of the liquid inside of liquid slugs. Minimum values PAGE 69 69 of time averaged wall pressure are observed at the tail of liquid slugs. The overall trend of Figure 4 7 shows that the pressure at the wall o f microchannel is smaller than that on the axis. Figure 4 7. Pressure distribution on wall and axis of the microchannel for Run 12 (t=0.07887 s). A similar pressure distribution has been reported by Kreutzer et al (2005) and Gupta et al. (2009). Pre ssure distribution on the wall is complicated however some overall trends can be explained here. Since bubbles travel faster than the liquid slugs, they create small vortexes after themselves that are corresponding to the pressure drop in the liquid slug. As can be seen in Figure 4 7, between two dashed lines, pressure in the liquid film has many oscillations due to instabilities on the gas liquid interface. The pressure drop in the middle of the film region is not noticeable however close to the end PAGE 70 70 of the film the film becomes thinner and the curvature of interface increases and makes a concave. The sharp concave at the end of liquid film creates a small vortex and a big pressure drop; this fact can be seen in the zoom region of Figure 4 7 and also can be seen more clearly in Figure 4 6. A linear interpolation on the pressure distribution is calculated and the gradient of the pressure is found to be 95.9 kPa/m for run 12. Table 4 4 compares the predicted pressure drop to those calculated by several sugges ted correlations. Table 4 4. Comparison of pressure drop with available empirical correlations and previous numerical simulations Different correlations Calculated pressure drop (Pa/m) Description ] ) Re ( 1 [ Re 16333 0 1Ca L D fs 1=0.07 184,000 fo 1 =0.17 Kreutzer et al. (2005) ) 1 )( Re Re ( 64) 1 ( cb a f a, b and c =f(Ls/D) 132,067 Mehdizadeh et al. (2009b) ) ) 1 ( ) 1 ( ( 32 ) (2 2 L G L GD G x P 29,620 Mehdizadeh et al. (2009a) Numerical simulation, Gupta et al. (2009) 70,700 Gupta et al. (2009) Present work 95,960 Present simulation The predicted pressure drop is 1.5% bigger than that obtained from the experimental correlation suggested by Kreutzer et al. (2005). Also, the pressure drop calculated from the suggested correlation of Mehdizadeh et al. (2009b), w hich was obtained from the simplified geometry of the liquid slug, is 37% bigger than that predicted by the present study. Even though this work accurately predicted the pressure PAGE 71 71 drop, a detail explanation of pressure distribution along the channel needs m ore studies. The physics of the instability in the liquid film region and corresponding effects on the flow field will be study in the next steps. Conclusions In this chapter, slug flow in microchannels of different sizes and flow conditions has been simulated using the volume of fluid (VOF) method. The elements in the vicinity of the interface were refined using dynamic mesh adaption while the interface between phases was tracked. Different mesh resolutions were tried to find the minimum number of element s required to capture the liquid film around the gas slug. The former was successfully captured employing 96% less mesh compared to previous numerical investigations. This may prove to be significant since the VOF method is inherently computationally costl y and small time steps (in the order of 107) are necessary to capture high frequency interface deformations and fluctuations. The liquid film was not successfully captured in many previous numerical studies, mainly due to employing course meshes and/or to minimize computational time. Guidelines to set the refinement thresholds for dynamic mesh adaption were developed. The negative effects of using a large elemental aspect ratio have been described. The predicted shapes of gas slugs were found to be in good agreement with the experimental flow visualization data. Also, values of the volume fraction, slug velocities, and liquid film thickness were compared to those obtained from available empirical correlations along with data obtained from a previous numeric al simulation study and were found to be in good agreement. This study can be extended to model heat and mass transfer of slug flows in microchannels for different chemical and mechanical applications. PAGE 72 72 CHAPTER 5 A COMBINED ANALYTICAL NUMERICAL MODEL FOR A TWO PHASE FLOW THROUGH A SUDDEN ARE A CHANGE IN MICROCHANNELS Outlook In this chapter two new analytical models have been developed to calculate twophase slug flow pressure drop in microchannels through a sudden contraction. Even though many studies have been reported on twophase flow in microchannels, considerable discrepancies still exist, mainly due to the difficulties in experimental setup and measurements. Numerical simulations were performed to support the new analytical models and to expl ore in more detail the physics of the flow in microchannels with a sudden contraction. Both analytical and numerical results were compared to the available experimental data and other empirical correlations. Results show that models, which were developed b ased on the slug and semi slug assumptions, agree well with experiments in microchannels. Moreover, in contrast to the previous empirical correlations which were tuned for a specific geometry, the new analytical models are capable of taking geometrical par ameters as well as flow conditions into account. A nalysis Simplified Model of Two Phase Slug Flow Pressure Drop in Microchannel Damianides and Westwater (1988), Fukano and Kariyasaki (1993), Triplett et al. (1999), Zhao and Bi (2001) and Kawahara et al. (2 002) developed overall twophase flow regime maps taking into account a wide range of parameters in microchannels. Kawahara et al. (2002) did not observe any bubbly or churn flow patterns in their developed map for a 100m microchannel. On the other hand, for a wide range of gas and liquid superficial velocities, the slug ring, ringslug, multiple and semi annular flows were observed. PAGE 73 73 Figure 5 1. Left: Experimental image, Serizawa et al. (2002). Right: Schematic picture of two phase slug flow in a straight microchannel For slug flow it can be assumed that pressure drop is due to the additive pressure drop of separate liquid and gas phases. Here we assume that a gas slug and a liquid slug follow each other, that there is no mass transfer between the phases, that the gas phase is incompressible in the channel (This assumption would not be valid for the abrupt area change region since the pressure changes noticeably in a short distance in that region), and that these slugs have the same veloci ty (see Figure 5 1). This is completely different than homogenous flow in which both phases are so mixed to each other to behave as new homogenous flow in which it is normally assumed that both phases have the same velocity. In other words, it is assumed t hat the slip ratio is unity for homogenous flow while the slip ratio definition for slug flow in microchannels when there is periodic dry wall is not clear. The mass flow rate of gas and liquid in a microchannel can be written as G G GAJ m. (5 1) L LAJ mL. (5 2) where JG and JL are the superficial velocities of the gas and liquid phases, respectively. A is the cross sectional area of the circular microchannel and G and L are the densities PAGE 74 74 of the gas and liquid phases, respectively. Therefore the total mass flow rate can be expressed as ) (. L L G G tJ J A m (5 3) Si nce the mass flux across the gasliquid interface is zero in the interfac e frame of reference, the velocity of both phases can be assumed to be the same, uG=uL=u ; therefore the velocity of the slugs in the stationary frame of reference can be written as L G L G L L G GG J J u ) 1 ( ) 1 ( ( 5 4) in which is the gas phase void fraction that is equal to the homogenous gas void G is the total mass flux (kg/s.m2) of the two phase flow. In the slug flow, for each phase, we have assumed that the flow is fully de veloped except close to the interfaces. This assumption may be reasonable only for relatively long slugs. For each phase, assuming a parabolic laminar, fully developed velocity profile, the pressure gradient would have the same form as for laminar pipe flo w, with a time weighting correction. ) 1 ( 32 32 ) (2 2 2 D u D u x PL L G G (5 5) This means that at a certain location of the channel, in portion of the time in which only the gas phase exists, the pressure gradient is due to the gas flow friction and in (1 portion of the time, in which the liquid phase contacts the walls, the pressure gradient is due to liquid friction. Using Equation ( 5 4), assuming uLand uG would be the same, the above equation for slug flow can be rewritten as PAGE 75 75 ) ) 1 ( ) 1 ( ( 32 ) (2 2 L G L GD G x P (5 6) Kawahara et al. (2002) compared six different relations for twophase homogenous viscosity models to find the pr essure gradient. Five models had more than 100% error; however, the homogenous model of Dukler et al. (1964) for viscosity that is the same as that presented here (Eq. 5 6), ( G+(1 L), had the best agreement with the experimental data (within 20% ). However, in the procedure of finding Equation ( 5 6) we did not use any homogenous flow assumption. The fully developed and parallel flow assumptions for each phase is not vali d near the gas liquid interface and because of this fact the experimental pressure drop data show the factor of 30.08 instead of 32.0 in Equation ( 5 6), which is 6% lower than that of the conventional correlation (Kawahara et al. 2002). Flow Area Expansion Analysis, Single Phase Figure 5 2 shows a schematic of the flow through a typical flow area expansion. The target control volume stretches from the area change position to a downstream location at which it is assumed that the flow reaches a fully developed state. Figure 5 2 Schematic of single phase flow through an expansion Applying the one dimensional conservation of energy equation for the control volume shown in Figure 5 2 leads to the following: PAGE 76 76 2 ) 2 1 ( ) 2 1 (2 1 2 3 3 3 3 2 1 1 1 1u k u p u pave e (5 7) So } 2 /{ )} 2 1 ( ) 2 1 {(2 1 2 3 3 3 3 2 1 1 1 1u u p u p kave e (5 8) where ke is the expansion loss coefficient and ave is the average density in the expansion area. 1 and 3 are kinetic energy correction factors, which can be defined as 3 3 aveAu dA u (5 9) where uave is the average velocity. When u=uave that is uniform flow, = 1 0 For laminar flow though a round pipe = 2 0 and for turbulent flow = 1 05 Now, applying the one dimensional conservation of momentum equation 3 1 3 1 p pF F M M (5 10) 3 3 3 1 2 3 3 3 3 2 1 1 1 1A P A P u A k u A kd d (5 11) where 1 and 3 are momentum correction factors, which are defined as 2 2 ave Au dA u d k (5 12) When u=uave that is uniform flow = 1 0 For laminar flow though a round pipe = 1 33 and for turbulent flow ( 1/7 law profile), = 1 02 By dividing both sides of Equation ( 5 11) by A3, ) 1 1 3 3 2 ( 2 1 2 3 2 3 3 2 1 1 1 3 1 d k d k u u d k u d k P P (5 13) where the area ratio is defined as PAGE 77 77 3 1A Aa (5 14) T herefore Equation (5 8) can be written as ) 21 /( ) ( 22 1 2 3 3 1 1 1 1 2 3 3 2 2 1u k k u kave a d a d a e (5 15) For an incompressible flow this equation leads to 1 2 1 ) 3 3 2 ( 2 d k d k a e k (5 16) and if a flat velocity profile is assumed (which is common practice for the definition of loss coefficient), this equation leads to the BordaCarnot relation 2 ) 1 ( a e k (5 17) In this study, to have consistency with the collected data of Toufik (2008), the same simple relation was used to compare our modeling results to the experimental data. Therefore, both modeling and experimental loss coefficients can be found by combining Equations ( 5 8) and ( 5 17) as follows: 1 2 2 1 2 1 a u P e k (5 18) where is either the modeling or experimental pressure dif ference across the expansion area. In this relation several simplifying assumptions were applied. All correction factors of momentum and kinetic energy were assumed to be unity. Moreover, this relation inherently assumes that the flow is incompressible, which is valid for the liquid phase. However, for the gas phase this assumption may not be valid under PAGE 78 78 some conditions. When frictional loss is included, as it must be for a very long and narrow pipe, the incompressible flow analysis previously considered ap plies until the pressure drop does not exceed 10% of the initial pressure (ASHRAE, 2001). Since compressibility makes the analysis very complicated, Toufik (2008) assumed that the gas phase behaves as incompressible and Kawahara et al. (2002) used the average gas density between inlet and the outlet conditions of the pipe to calculate the loss coefficient. However, in the next section we will show that the incompressible assumption is not valid for two phase slug flows through a sudden area change, since t he density ratio of phases is on the order of O(1000), providing a sharper gradient of pressure along the channel. Flow Area Contraction Analysis, Single Phase It is common practice that the converging section of the flow (until the vena contracta) in which deceleration takes place from the vena contracta to the fully developed flow region can be modeled as flow through a sudden expansion (Kays, 1950). According to this assumption, expansion after the vena contracta to the downstream region can be modeled similarly to that in the flow expansion that was mentioned in the previous section. For incompressible flows it can be written as: 2 2 2 ) 2 1 1 2 2 1 ( 1 c C c C c C a d k a c k (5 19) PAGE 79 79 Figure 5 3. Schematic of single phase flow in Contraction where = / 1, is the Venacontraction coefficient. Geiger (1964) suggested an experimental correlation for 5371 0 ) 1 ( 08 2 1 1 a a c C (5 20) In this study to have consistency with the collected data of Toufik (2008), the same simple relation was used to compare our modeling results to the experimental data. Therefore, both modeling and experimental loss coefficients have been calculated based on the = 1 assumption as 1 2 2 1 2 1 a u P c k (5 21) Two Phase Pressure Change across Area Change (Conventional Models) It is common practice that using the same analysis for single phase flows, the pressure drop in sudden expansions and contractions in twophase flows wit hout phase change and for flat velocity profiles, can be expressed as follows (see Abdelall et al. 2005, Kawahara et al. 2002, and Toufik, 2008): PAGE 80 80 ) (1 2 3 2 2 1 3 2 1 2 a a eG P P P (5 22) where, G1, is the mixture mass flux in the smallest channel and subscr ipts refer to Stations 1 and 3 in Figure 52. is the momentum density which was defined by Lahey and Moody (1993) according to the following equation: 1 2 2] ) 1 ( 1 [ G Lx x (5 23) where x is the vapor mass quality and is the void fract ion of the flow. Toufik et al. (2008) assumed that if both phases are incompressible, x and both would remain constant during the flow area expansion and contraction. So Equation ( 5 22) can be simplified to ) 1 (2 2 1 a eG P (5 24) In the next section we will show that in microchannels assuming constant during the expansion or contraction may not be valid since each phase accelerates or decelerates at a different rate because of the different densities, wall shear forces and visc osity. In order to use the above equation, a closer equation needs to correlate vapor mass quality to the void fraction. The quality void relation for one dimensional flows are related to the slip ratio S according to the following: G Lx x S 1 1 (5 25) PAGE 81 81 For homogenous flows S=1 Two phase flows across sudden contractions are considerably more complicated than those across sudden enlargements. In the flow area contraction in two phase flows, it is still unclear whether the charac teristics of the vena contracta in two phase flow are the same as those of single phase flows. However, for twophase flows, in analogy with single phase flows using the vena contracta concept Collier (1972) and Hewitt et al. (1993) suggested the followin g: )] 1 ( ) 1 ( 2 [' 1 2 1 2 3 2 2 1 2 2 2 1 c a c h cC C G P (5 26) where G1 is the mixture mass flux and Cc is the coefficient of contraction which is a function of the area ratio. Chisholm (1983) recommended 1 1 639 0 1 a cC (5 27) where h is the average homogenous flow density between Points 2 and C (the minimum area) which can be found from the average slip ratio between these points G L h ) 1 ( (5 28) and 2 1 2 2 3 2 2 3] ) 1 ( ) 1 ( [ G Lx x (5 29) Also the m omentum density is defined in Equation ( 5 23) and void fraction is defined as x x S xL G L ) 1 ( (5 30) PAGE 82 82 in which, x, is the vapor mass quality. With the same analogy for the pressure drop in the pipe, effective mixi ng caused by sudden contraction may justify the assumption of homogenous flow, and leads to L L cP P (5 31) )] 1 ( ) 1 1 [( 22 2 2 2 1 a c L cC G P (5 32) L G L Lx / ) ( 1 (5 33) However, Schmidt and Friedel (1997) have shown that the vena contracta phenomenon may not occur in twophase flow at all. In the next sections we will introduce an analytical model for slip ratio value in vena contracta location. We will then introduce a new model for pressure drop in contraction for two phase slug flow in microchannels. We will show through the latter model that that the homogenous flow assumption or It is common practice that the converging section of the flow (until the vena contracta) in which deceleration takes place from the vena contracta to the fully developed flow region can be modeled as flow through a sudden expansion (Kays, 1950). According to this assumption, expansion after the vena contracta to the downstream region can be modeled similarly to that i n the flow expansion that was mentioned in the previous section (see Equation ( 5 19)) T he constant void fraction in contraction is not a perfect assumption for flow through the contraction in microchannels. PAGE 83 83 New Analytical Model for Void Fraction in Contra ction As we have shown for pressure drop in the straight microchannel pipe flow, if the flow behaves like the slug flow, the nature of flow is similar to a pulsing flow and the frequency of the pressure pulses in the vicinity of the contraction dep ends on the length and velocity of liquid and gas slugs. In all cases, not only does void fraction not remain constant but the two phase flow regime may also change from that of a big pipe to that of a small pipe. Figure 5 4 shows a schematic of the flow i n a microchannel with a sudden contraction. We are assuming that the liquid slugs are incompressible while gas slugs are compressible. Even though the constant static pressure far upstream of the contraction is not a perfect assumption for this pu lsating type flow, because of the friction, damping, and simplicity we have set the far upstream pressure constant while the downstream pressure was allowed to oscillate. Since we are looking for the pressure difference between the up and downstream locat ions, we expect that this assumption would not have a major effect on the final results since with a simple shifting we can set either the downstream or upstream pressure as the reference pressure Even though obtained results show that all movement s and oscillations have been affected by gas and liquid slugs length in the bigger pipe published experimental work does not reveal much information on these lengths or correlate them to the other properties. For a fixed gas slug length, LG, liquid slug length, LL, is just a function of quality (or void fraction). However, from experimental photographs in the literature, gas slug length varies from 1 to 15 times the pipe diameter (see Figure 5 1) and for this case we have assumed the initial length of the gas slug to be seven times that of the bigger pipe diameter, D3, and for different qualities, the corresponding liquid slug lengths PAGE 84 84 could be calculated. At this point we recommend for that subsequent experimental studies provide detailed information on th e time averaged values of slug lengths under different conditions. The proposed model assumes that the liquid slug hits the facing wall of the contraction at t=0 and after this point, the velocity of the center of mass of the liquid slug and those of the front and tail of this slug as well as the pressure and velocity of the gas slug following the liquid slug would be monitored until the latter passes the contraction and reaches a steady state condition. Figure 5 4 Schematic of two phase slug flow in a microchannel facing a contraction A brief algorithm of the present model is shown in Figure 5 5. The gas slug pressure behind the target liquid slug is characterized by an equation of state. If the ideal gas assumption is assumed to apply on the gas slug, the pressure can be expressed as RT x x P ) ( ) ( (5 34) The liquid slugs center of mass location is determined by ) ( 22 3 3 2 1L A L A m Xf l l c (5 35) where l and ml are liquid slug density and mass, r espectively. Lf and L3 were defined in the schematic in Figure 5 4. The center of mass velocity, uc, can be obtained by PAGE 85 85 differentiating Eq. ( 5 35) with respect to time. Applying the continuity equation to the liquid slug, the front liquid slug velocity in the smaller pipe, uf, can be correlated to the center of mass velocity by this nonlinear relationship ) ( ) (3 1L L u A m uf c l l f (5 36) When liquid slug is going through the contraction, the pressure in the vicinity of the front walls reaches values very close to the stagnation pressure. The actual integration of the numerical results (CFD modeling which is not presented in this chapter ) over this surface reveals that the average pressure on this wall is almost 94 percents of the stagnation pr essure and slightly changes by changing the area ratio. However, for simplicity we have assumed that this pressure is the same as the stagnation pressure. On the other hand nonlinear shear forces act on the liquid slug. For laminar flows, the friction fac tor, f is a sole function of the Reynolds number. For turbulent flow, which is the case here (based on the real Reynolds number), since the diameter is very small, the fully rough wall assumption is appropriate. Thus, using the Colebrooks natural roughne ss function that is independent of the Reynolds number, the friction factor can be expressed as ) / log( 2 14 1 1D f (5 37) Conventional values of turbulent and laminar Reynolds numbers in two phase slug flows may be questionable since t here is no obvious length scale to define the Reynolds number. One can define it based on the diameter of the pipe while because the small length of the slugs (the same order of the diameter) one may prefer to define it based on the average length of the s lugs. Moreover, it seems that if we fix our frame of PAGE 86 86 reference on the moving slugs, because of the short length of the slugs, most or all of the slug length would be in the developing region. This indicates that the length of the slug plays an important role in slug flows. While this seems to be the case, no one seems to have reported on the time averaged values of the slug length. Because of this fact, even for laminar flow, since it is expected that for the most part the slug is located in the developing region in which the velocity gradient and consequently the overall friction factor are greater than what is expected from the parabolic velocity profile for laminar flows, Equation ( 5 37) may present a closer approximation of the real friction factor for s lug flows. As can be seen in the schematic figure and algorithm, since the net forces and moving mechanism are different before and after passing of the liquid slug tail from the contraction area, calculations should be performed separately for the se regions (see the algorithm) As long as the liquid slug tail travels in the bigger pipe, the front velocity and the center of mass velocity will correlate well with Eq. ( 5 36) and the net force acting on the liquid slug can be calculated as foll ows: ) ( ) (3 3 3 1 1 3 f f f aveF L D F L D A P P F (5 38) where, Ff, is the wall shear stress, P is the variable gas slug pressure and Pave is the average pressure acting on the front face of the contraction. They can be calculated from the following equations: f f L fu fu F8 1 (5 39) a i a s aveP P P ) 1 (3 (5 40) PAGE 87 87 In each time step, the liquid slug center of mass velocity can be calculated based on the momentum change and the net forces acting on the slug dt m F u ul i c i c 1 (5 41) where i, i+ 1 represent the current and next time steps, respectively. Then the front velocity, uf would be obtained from Eq. (5 36 ) and the tail velocity can be easily found by u3auf. With knowing the tail velocity, u3(x), the gas slug pre ssure can be calculated by finding the volumetric compression of the gas slug. Therefore, applying the equation of state for an ideal gas, the following correlation can be obtained )) ( /( ) (3 3 3 1u u dt L A RT m Pi f i g i (5 42) where mg is the gas slug mass (using the initially7D3 assumption for the gas slug length mentioned earlier). Since walls periodically wet and dry with the liquid and gas slugs, it is expected that for such a small diameter, the temperature of the gas slug would remain constant; hence, in dev eloping Equation ( 5 42), we have assumed that this process is almost isothermal. This seems to be valid for small compression ratios, but it may not be valid for cases that include phase change between vapor and liquid slugs where the mass of the gas slug would not be conserved anymore. These processes continue until the liquid slug tail passes the bigger diameter pipe and enters the smaller pipe. Then the net force acting on the liquid slug can be calculated from f L jF L D A P P F1 1 1) ( (5 43) wher e LL 1 is the liquid slug length in the smaller pipe. Practically Pj is the gas slug pressure subsequent to the liquid slug; however, for simplicity since the gas slug length PAGE 88 88 is very long in comparison with the diameter of the smaller pipe, it can be assumed that Pj is equal to the downstream pressure P1. So as we mentioned, calculation starts when the liquid slug hits the facing walls of the contraction at t=0 (s). At this stage Lf=0 (Figure 5 4) but the front velocity is uf=u3a and the tail velocity at this infinite small time is equal to u3 (boundary conditions). Because the gas viscosity is much smaller than that of a liquid, the gas slug can travel through in the contraction rapidly while the liquid slug (having t he higher momentum (Order (1000))) hits the front walls of the contraction and produces pressure waves that are directly transferred from the incompressible liquid slug to the gas slug and causes a phenomenon similar to that of a water hammer but with smal ler magnitudes. This is because gas slugs are much more elastic than liquid slugs and this fact reduces the strength of the pressure waves. Water slugs pass through the contraction with more restriction due to the larger friction and the smaller average pr essure difference on both sides of the liquid slug (Eq. ( 5 43)). This causes a small chocking and reduction in the liquid slug velocity. On the other hand, the slugs that follow come with higher velocities and compress the gas slug behind the liquid slug p assing through the contraction. The gas slug pressure, P increases until the liquid slug completely gets inside the smaller pipe. Then, the gas slug that follows, which has already been compressed, shoots the liquid slug in the smaller pipe. Because of th is fact, the shooting speed of the liquid slug exceeds the steady state velocity of the downstream flow, which leads to a decrease in the pressure of the gas slug that follows, while the wall shear force is also affecting the slug motion. PAGE 89 89 Figure 5 5. Su mmary of algorithm used for modeling the average void fraction in the vicinity of the flow area contraction This springdamper type (here spring behavior of gas slug is not linear) of oscillation may continue until the wall shear force damps this o scillation and the liquid slug velocity reaches the steady state condition of the downstream flow. PAGE 90 90 This phenomenon can be seen in more detail in Figure 5 6. In this figure the liquid slugs center of mass velocity, uc, slightly decreases at the initial sta ges d ue to the chocking (adverse front wall forces and excess shear forces) that was mentioned a. Figure 5 6. Gas slug pressure liquid center of mass velocity and liquid front velocity The liquid front velocity in the smaller pipe initially starts from the steady state value in the small pipe and decreases sharply due to the reduction of uc as well as the change in the center of mass position that is dynamically moving. After passing the contraction, both the slugs center of mass velocity and the front velocity become equal PAGE 91 91 to each other (it becomes a normal slug in the smaller pipe in which the center of mass velocity is the same as the front and the tail velocities) and oscil late with the pressure and reach a steady state (because of friction damping) after many oscillations. Since the pressure volume relation and damping forces are not linear, none of these oscillations are sinusoidal in nature, thus, all the main frequencies can be extracted with a Fast Fourier Transform algorithm. To find the gas void fraction at the vena contracta, knowing the location of the vena contracta is important. Although extensive experimental studies have been reported on air water twophase flows in mini and microchannels, there seems to still exist considerable discrepancies, largely due to the difficulties in experimental setup and measurements (Kandlikar, 2002 and He and Kasagi, 2008). To understand the detailed physics of slug flow in microchannels, advanced numerical simulations were performed to help develop a more precise analytical model. Numerical studies on singlephase air water flows, showed that for turbulent water flow the position of the minimum cross sectional area (due to the vena contracta) occurs at 0.25D1 downstream of the sudden contraction (this number for laminar flow is 0.15D1) and the minimum cross sectional diameters are 0.85D1 and 0.92D1 for turbulent and laminar flow, respectively. The conventional definition of the slip ratio is the ratio of the timeaveraged value of the gas phase velocity to the liquid phase velocity. For a small pressure drop it is common practice that the compressibility effects be neglected. However, this assumption is not valid in the vicin ity of the contraction where the pressure gradient is significant and more important since in two phase flows the density ratio between phases is typically in the order of 1000 and all changes take place in considerably shorter distances, thus, the PAGE 92 92 gradient of the gas density could not be neglected. With this analytical model the length of liquid and gas slugs can be determined at any time. Therefore, if compressibility effects of the gas slug were taken into account, the timeaveraged void fraction in the vena contracta may be expressed by the time averaging of the volume of the gas slug. Time averaged void fraction of the gas slug in the vicinity of the vena contracta location (x=0.25D1) can be calculated from 2 2 20 0 0) ( ) (t G t L t Gdt t V dt V dt t V (5 44) where VG and VL are the gas and liquid slug volumes, respectively. Time starts when the gas liquid interface reaches the vena contracta location, while t2 corresponds to the time that the tail of the two liquid and gas slugs that follow pass the vena contracta location. New Pressure Drop Model for Slug Flow through a Contraction in Microchannels Pressure drop across sudden contraction for twophase flow assuming overall long and short slug flow regimes can be determined by the same concept that was used for the singular pipe. Pressure drop through the contraction for a singlephase flow can be found using Equation ( 5 32). So for the flat velocity profile which is a reasonable assumption for slug flow in an portion of time, the pressure drop is due to the gas phase motion, and in (1) portion of time, the pressure drop is due to the liquid phase motion, so the timeaveraged pressure drop can be written as ) 1 ( ) ) 1 1 ( 1 ( 22 2 2 1 2 L G c cC u P (5 45) where Cc was d efined in Eq. ( 5 27) and u1 is velocity of slugs in the smaller pipe. This was defined in Eq. ( 5 4) and is the void fraction in the vicinity of the vena contracta, which can be found from the new analytical model of averaged void fractions presented PAGE 93 93 in t he previous section. It is well known that void fraction values around the vena contracta controls the pressure drop. Different empirical void fraction models can be found in the literature, however each model is valid for a specific regime or experimental condition. Moreover, these models are tuned for straight pipe flows. As mentioned earlier, for the sudden area change it is common to use one of those correlations to find an expression that makes a better fit for the experimental data (see Chalfi and Ghi aasiaan, 2008, Abdelall et al. 2005, Kawahara et al. 2002, Hewitt et al. 1993, and Collier, 1972). So, a good amount of the previous work is practically based on curve fitting of existing void fraction or slip ratio correlations for the straight pipe t o the experimental data associated with the flow in contractions. However, the one dimensional model discussed in the previous section showed that the gas phase accelerates quicker than the liquid phase in a sudden contraction due to the lower density and wall shear forces. These analyses led to introducing a new model for the void fraction in the vicinity of the contraction. In the following, results of the new models for void fraction and pressure drop (Eq. ( 5 45)) will be compared to the conventional mod els and experimental data. Table 5 1 shows the experimental data of Abdelall et al. (2005). In their experiment, they found the pressure drop using linear interpolation of the pressure data up and dow nstream of the contraction (Figure 53). Two pha se regime maps in microchannels are slightly different than those of macrochannels mainly due to the larger surface tension effects and the reduced gravitational forces. Kawahara et al. (2002) studied the flow regime map for a 100 m diameter pipe and com pared it to other available maps for microchannels. PAGE 94 94 The closest map to Abdelall et al. s (2005) test rig is the test case of Tripllet et al. (1999) whose data were obtained for a 1.1 mm diameter tube (close to the average diameter of Abdelall et al. s test rig). Therefore, it is expected that Figure 5 7 would show the overall flow regime map of this case study. It can be seen that all of the data in the bigger pipe (upstream of the contraction) are in the slug or ring slug flow modes, which makes th em more appropriate for our model. The downstream regime, on the other hand, is in multiple zones, which means that all bubbly, ring slug, slug annular, churn, annular and slug flows can be observed and our model for the slug flow may not be appropriate. H owever, since there is a transitional process for regime change and we are just looking to examine the vena contracta location, which is located a short distance downstream of the contraction, it can be expected that the flow regime would not change from t he initial slug regime state to a completely multiple zone regime in this short distance. Also, even if it changes, since we are time averaging void fractions through the contraction process, the side effects of this issue would be damped with this integration. Figure 5 8, shows the relation between the void fraction and homogenous void S=1 S. Armand (1946) Based on the experimental data of a 100 m diameter tube Kawahara et al. (2002) introduced the following equation: 5 0 5 097 0 1 03 0 (5 46) PAGE 95 95 Table 5 1 Two phase flow sudden contraction experimental data Abdelall et al. (2005) (Kg/s) (kg/s) (pa) 0.00115 5.4E 06 5251 1754.4 0.00115 6.6E 06 5353 1756.2 0.00115 8.6E 06 5606 1759.3 0.00115 1.1E 05 6139 1762.4 0.00115 1.5E 05 6235 1768.9 0.00115 2.0E 05 6349 1776.2 0.0017 5.4E 06 9887 2580.6 0.0017 8.1E 06 10826 2584.7 0.0017 1.1E 05 10904 2588.8 0.0017 1.3E 05 11114 2592.7 0.0017 1.6E 05 11470 2596.8 0.0017 1.9E 05 11579 2601.2 0.00232 4.6E 06 17160 3516.1 0.00232 5.8E 06 17185 3517.9 0.00232 7.3E 06 17678 3517.9 0.00232 8.8E 06 18194 3520.2 0.00232 1.2E 05 18280 3526.6 0.00232 1.4E 05 18939 3530.7 0.00232 1.7E 05 19179 3535.2 0.00232 2.0E 05 20378 3539.8 0.00258 4.9E 06 21398 3908.5 0.00258 6.2E 06 22112 3910.5 0.00258 7.8E 06 22892 3912.8 0.00258 9.3E 06 23325 3915.1 0.00258 1.2E 05 23963 3919.5 0.00258 1.5E 05 24565 3923.8 Figure 5 7. Two phase flow regime map for a 1.1 mm diameter tube, Triplett et al. (1999) Bold big circles are Abdelall et al. s (2005) data which were used for the analytical modeling (Hollow big bol d circles are for the larger diameter pipe and solid big bold circles are for the smaller diameter pipe) PAGE 96 96 Figure 5 8. Relationship between homogeneous void fraction and void fraction for different models in the vicinity of the flow area contraction For the annular flow when the minimum entropy generation rule is applied, Zivi (1964) introduced the following expression for the slip ratio 3 / 1) (G Lc S (5 47) where c=0.7 in the modified model. If curve fitting were to apply on the data points to minimize the error in the pressure drop using the presented model in Eq. ( 5 45), the PAGE 97 97 following correlation can be proposed for the average void fraction in the vicinity of the con 5 0 5 097 0 1 47 0 5 0 (5 48) All of the above correlations (except the curve fit) were tuned for the straight pipe flow. However, the new analytical dynamic model for the void fraction (Equations ( 5 34) to (4 5) ) was obtained based on the slug and semi slug flow assumptions in which the void fraction dynamically changes in the transitional area due to the difference in acceleration or deceleration of each phase. Calculated averaged void fractions with the proposed dynamic model are shown in Figure 5 8. As can be seen, for this area ratio and flow conditions, the data correlate well in the large volume fraction region and correlate well with the Armand correlation for the smaller volume fraction region. Moreover, against the other models, the dynamic models data are not on a line or a special curve and the data depend on the flow conditions as well as the geometrical constrains. It is interesting that the results of the new analytical model show outstanding agreement with the experimental curve fit, so it seems that it can describe the physics of this phenomenon quite well. Figure 5 9 illustrates a comparison of different models with the experimental data for the pressure drop across the sudden contraction in microchannels. Before taking a look at the details, it can be seen that all of six figures show the robustness of Equation (5 45) for slug flows PAGE 98 98 (a) (b) (c) (d) (e) (f) Figure 5 9. Comparison of six different void fraction models a) New analytical model (Equations ( 5 34) to ( 5 45)) b) Armand correlation c) Kawahara et al. s correlation d) Homogeneous model e) Pure curve fit of experimental data f) Zivis correlation PAGE 99 99 If we compare the maximum error in these figures with the latest conventional models (Equation ( 5 26)), that wer e applied by Chalfi and Ghiaasiaan (2008) and Abdelall et al. (2005) who reported 500% error for the homogenous model and 100% error for the slip flow models, the robustness of the proposed model for pressure drop should be clear. These figures show th at the presented model (Equation ( 5 45)) is robust enough that the maximum error would not exceed than 310% (for the homogenous model) vena contracta location. In Figure 5 9.d, it can be seen that the homogenous model has the worst accuracy. The error in the Kawahara et al. s correlation in Figure 5 9.b is also significant. Chalfi and Ghiaasiaan (2008) and Abdelall et al. (2005) reported 100% error for the pressure drop usi however, in Figure 5 9.f it can be seen that Equation ( 5 45) could reduce the maximum error to 50%. Armands correlation in Figure 5 9.b which was used by Chalfi and Ghiaasiaan (2008) and Abdelall et al. (2005) and was their most accurate model which showed 30% error, with the presented model shows 20% error. With a curve fit and tuning based on the experimental data defined by Equation ( 5 48), the error in the pres ented model can be reduced to 8% which is shown in Figure 5 9.e. As was mentioned earlier, all of the above models are either based on the empirical model or the homogenous model and more importantly they were tuned for PAGE 100 100 the straight pipe flow. But the new purely analytical model for the void fraction could have a good agreement with experimental data (see Figure 5 9.a). After curve fitting the data (Figure 5 9.e) this analytical model has the second most accurate prediction and reduced error to 20% More importantly, since this is an analytical model, all changes in the flow properties and different geometrical constrains can also be monitored and studied. Furthermore, it can provide better understanding of the physics taking place in slug flows fac ing the contraction that were harder to study before using the more conventional models. Figure 5 10. Nondimentional pressure drop versus two phase Reynolds number for the experimental data of Abdelall et al. (2005) Pressure drop in nondimensi onal form versus the two phase flow Reynolds number can be seen in Figure 5 10. Nondimensionalization is obtained based on the new correlation introduced in Equation ( 5 45) in which Equation ( 5 48) is used to define PAGE 101 101 G So, the nondimensional pressure drop is defined as ) 1 ( 2 1 ) ) 1 ( (2 2 exp a L GG P P (5 49) While the twophase Reynolds number is defined a s 1 1 2)) )( 1 ( ) ( ( Re D uG G L L (5 50) It can be seen that all nondimensional values of the pressure drop are in the range between 0.945 and 1.16 (close to unity). Also the data are showing a weak function of the twophase Reynolds number which indicates that this nondimensional parameter for the pressure drop is the main nondimensional parameter. To have more precise prediction of the pressure drop, some curves can be introduced to distinguish the data based on their void fraction values. It can b e seen that for any given two phase Reynolds number there are two possible regions for the nondimensional pressure drop. This may come from the fact that for any given twophase Reynolds number in reality there are two different gas liquid slug lengths. Th e smaller length of liquid slug corresponds to the larger pressure drop and vice versa. So it seems that with additional experimental data about the slug length, the pressure drop in slug flows in microchannels can be more precisely predicted to even have better accuracy than the current accuracy (8% error). So to achieve this goal many experiments should be conducted on the slug lengths. PAGE 102 102 Conclusions In this chapter a new analytical model for void fraction calculation was developed. The model reveal ed some important phenomena that are occurring in two phase slug flows through a sudden area contraction. Against the previous models which used slip ratio correlations of straight pipes for the flow in the contraction, the new model directly tackled this problem and provided accurate predictions of the void fraction. Moreover, based on the new model, we were able to find the value of the real void fraction in the vicinity of the vena contracta, which helped to develop another correlation for the pressure drop in twophase slug flows through the area contraction. Results showed excellent progress in pressure drop prediction. This new model could reduce the minimum (60%) and maximum (500%) errors of the more conventional models to a minimum of 8% and a maxim um of 310%. PAGE 103 103 CHAPTER 6 NUMERICAL SIMULATION OF THERMOFLUID CHARACTERISTICS OF T WO PHASE SLUG FLOW IN M ICROCHANNELS Outlook A fundamental study of heat transfer characteristics of twophase slug flow in microchannels is carried out employing the Volumeof Fluid (VOF) method. Despite of the fact that numerical simulations of twophase flows in microchannels have been attempted by many investigators, most efforts seem to have failed in correctly capturing the flow physics, especially thos e pertaining to the slug flow regime characteristics. The presence of a thin liquid film in the order of 10 m around the bubble is a contributing factor to the above difficulty. Typically, liquid films have a significant effect on the flow field and heat transfer characteristics. In the simulations reported in this chapter the film is successfully captured and a very high local convective heat transfer coefficient is observed in the film region. A strong coupling between the conductive heat transfer in th e solid wall and the convective heat transfer in the flow field is observed and characterized. Results showed that unsteady heat transfer through the solid wall in the axial direction is comparable to that in the radial direction. Results also showed that a fully developed condition could be achieved fairly quickly compared to singlephase flows. The fully developed condition is defined based on the Peclet number (Pe) and a dimensionless length of the liquid slug. Local and timeaveraged Nusselt numbers for slug flows are reported for the first time. It was found that significant improvements in the heat transfer coefficient could be achieved by short slugs where the Nu number was found to be 610% higher than in single phase flows. The study revealed new fi ndings related to slug flow heat transfer in microchannels with constant wall heat flux. PAGE 104 104 Introduction As we mentioned in Chapter 2, Walsh et al. (2009) seem to be the only group of researchers who had performed experimental studies on heat transfer of slug flows in microchannels. One of the objectives of this chapter is to attempt simulating the aforementioned experiment employing the Volumeof Fluid (VOF) method. In the chapter heat transfer of the slug flow will be compared to that in the conventional p lug and Poiseuille flows where a constant heat flux boundary condition is applied. Figure 6 1 shows a schematic of such flows. Figure 6 1 Schematic of thermal boundary layer development for plug flow (left) and Poiseuille fl ow (right) Using a simple energy balance the mean temperature, Tm, can be defined as: G G L L s ave w mcp m cp m x p q T T. 0 ( 6 1) where To is the inlet temperature, ave wq is the average heat flux, ps is the perimeter of the microchannel, x is the axial distance from the edge of the heater, and Lm and Gm are the water and liquid mass flow rates, respectively. The dimensionless position d ownstream of the heated section is defined by the reciprocal of the inverse Graetz number according to the following equation : PAGE 105 105 1 *Pr Re / Gr DPe x D x xD ( 6 2) where D is the diameter of the microchannel, Pe is the Peclet number, Pr, is the Prandtl number and Gr is known as the Graetz number. For a fully developed laminar flow in a circular pipe with a constant heat flux boundary condition, Graetz et al. (1883) reported an analytical solution for the temper ature field. For the fluid with Prandtl number greater than one, Muzychka and Yovanovich (2004) introduced a correlation for the Nusselt number which is valid in both the thermal entry and fully developed regions as follows: 5 / 1 5 3 / 1 5 1] ) 302 1 ( 36 4 [ x Nu ( 6 3) For inviscid flows or very low Prandtl number fluids, a flat velocity profile ( p lug flow) is a reasonable assumption. For such flows, Muzychka and Yovanovich (2002) introduced an expression for the Nussel t number that is valid for both the thermal entry and fully developed regions as follows: 2 / 1 2 2 / 1 2 1 ,] ) 886 0 ( 96 7 [ x NuPlug ( 6 4) In what follows, heat transfer characteristics of slug flows will be compared to t he above correlations in addition to relevant experimental data. Problem Formulation and Solution Procedure Governing E quations: The Volumeof Fluid (VOF) method is employed here to simulate gas liquid twophase slug flows in microchannels. This method wa s developed by Hirt and Nichols (1981) In it, two or more phases are not allowed to interpenetrate each other. The PAGE 106 106 summation of the volume fractions of all phases in each control volume is set to unity. The VOF method is probably one of the most convenient methods relative to other twophase flow modeling methods. The method has reasonable accuracy, is less computationally intensive (compared to the Eulerian model), is relatively simple to use, and can solve highly complex free surface flows. In this chapter both the liquid and gas phases are assumed as incompressible. The method is structured to solve a single momentum equation throughout the domain, while the resulting velocity field is shared among the phases. The momentum equation, shown below, is dependent on the volume fractions of all phases through the properties and ( 6 5) w here and in the above equation can be calculated from the following equations: ( 6 6) and (6 7) where and are the volume fractions of the liquid and gas phases, respectively. The term in Equation ( 6 5) represents surface tension and wall adhesion effects on the flow field. The surface tension is due to the strong intermolecular attract ive forces among molecules in a fluid. These forces contribute to SF TF g v v p v v v t )] ( .[ ) .( ) ( SFF G G L L G G L L G G G L L L L G SFF PAGE 107 107 minimizing the surface area of the slug. The wall adhesion is due to the stronger attractive forces among liquid molecules relative to their attraction to the wall. This causes the fluid to create some contact angle with the wall. In this chapter the Continuum Surface Force (CSF) model proposed by Brackbill et al. (1992) is used to model surface tension. The surface curvature is computed from the local gradients in the surface normal at the interface. Thus, the source term, in the momentum equation can be specified as follows (see Br ackbill et al. 1992): ] ) ( 5 0 [ G L G G L L SFn F ( 6 8) w here n is the surface normal vector and is the surface curvature Ln ( 6 9) )] ( ) [( 1 ) ( n n n n n n (6 10) A wall adhesion angle in conjunction with the surface tension model is also included in this simulation. The model is again taken from the work of Brackbill et al. (1992) The contact angle that the fluid makes with the wall boundary is used to adjust the surface normal in cells near the wall. This so called dynamic boundary condition results in the adjustment of the curvature of the surface near the wall. If is the contact angle at the wall, then the surface normal at the live cell next to the wall is sin cos w wt n n (6 11) w here wn and wt are the unit vectors normal and tangential to the wall, respectively. Some researchers believe that the wall adhesion properties, such as wall contact angle, play an important role in slug flow simulation. However, this may not be entirely true SFF PAGE 108 108 (Gupta et al. 2009). The role of the wall contact angle is only significant when both of the phases are in contact with the solid wall or the liquid film around the gas bubble is so thin th at the van der Waal forces can act across the film. This may be the case for T or Y type junctions where both phases remain in contact with the wall but once a slug flow is formed with a liquid film at the wall, wall adhesion plays no role in determining t he bubble shape. Moreover, the modeling of wall adhesion requires inclusion of the contact angle hysteresis effect to obtain meaningful results (Gupta et al. 2009). The tracking of the interfaces between the phases is achieved through solving the continuit y equation for the volume fraction of the second phase according to the following equation (see Brackbill et al. 1992): (6 12) The term is the mass transfer from the gas phase to the liquid phase and is the m ass transfer from the liquid phase to the gas phase. T he mass transfer between gas and liquid phases has been neglected. The volume fraction equation, Equation ( 6 12), will not be solved for the primary phase; T he primary phase volume fraction will be comp uted based on the following constraint: (6 13) Also a single energy equation will be solved through the domain and the results will be shared among the phases according to the following equation: (6 1 4) ) ( ) .( ) (. LG GL L L L L Lm m v t GLm. LGm. 1 G L )] ( .[ )) ( .( ) ( T k p E v E teff PAGE 109 109 where the energy term E in the above equation can be obtained by the following mass averaging method: (6 15) where the energy Ei for each phase is based on the specific heat of that phase and the shared temperature among the phases. In Equation ( 6 14) the effective thermal conductivity, keff is shared by the phases Test Cases Examined It seems that the only experimental heat transfer study of two phase slug flows in microchannel is reported by Walsh et al. (2009). We aim to numerically simulate their experimental results. The test consisted of a 6 cm microchannel with 1.5 mm internal diameter and an external diameter of 2 mm, giving a solid wall thickness of 0.25 mm. The solid stainless steel wall is accounted for in the simulation in order to be able to model the effects of the axial and the radial conduction heat transfer on the fluids heat transfer characteristics. Table 6 1 Property v alues e mployed Material Density (kg/m3) Specific heat (J/kg K) Thermal conductivity (W/m K) Water 998 4182 0.6 Air 1.22 1006.4 0.0242 Stainless steel 8030 502.2 16.27 At the inlet, two parallel streams (water and air) enter the channel with an inlet void fraction equal to the homogenous void fraction. Properties of the materials used in this study are shown in Table 6 1. n i i i n i i i iE E1 1 PAGE 110 110 Table 6 2. Test c ases e xamined Test Case Pipe diameter D (mm) Liquid flow rate (ml/min) Air flow rate (ml/min) Average wall heat flux (w/m 2 ) slug length/pipe diameter (L s /D) Capillary number ) (G L lJ J Ca 1 1.5 8 16 8,850 0.9 0.0031 2 1.5 8 16 8,850 1.7 0.0031 By changing the inlet cross sectional area for water and air streams, two slug lengths could be obtained. These have been compared to the experimental r esults of Walsh et al. (2009). A brief summary of the cases is given in Table 6 2. Model Verification and Validation Hydrodynamic characteristics of slug flows were extensively validated by Mehdizadeh et al. (2009b) and Mehdizadeh et al. (2009c). T hey showed that the film thickness, bubble shape and velocity, void fraction and pressure drop could be accurately simulated using the mesh adaption methods. They performed their simulations with 96% less mesh and with a superior accuracy than the only suc cessful study available in the literature (see Gupta et al. 2009). Figure 6 2 shows a comparison of the experimental data obtained by Walsh et al (2009) with the present simulations (Case 2 of Table 6 2 ). As can be seen, results agree well with t he experimental observations, where the wall is seen to experience oscillatory temperatures. Walsh et al. (2009) observed these oscillations when the wave length was of the order of the slug length. Oscillations were clear for the long slugs while for the short slugs Walsh et al (2009) could not clearly see the oscillatory behavior except in the vicinity of the thermal entry region. We assume that since they were measuring temperature from outside the pipe (using an infrared camera), it was hard for them t o capture the oscillations for the short slugs. In those cases, the PAGE 111 111 frequency of temperature oscillation is high, with the resulting effect of a high wall thermal storage and the dampening of the highfrequency oscillations. The circulation in the slugs pumps the cold liquid from the axis of the slug toward the wall. This circulation enhances the heat transfer and causes the temperature to oscillate near the wall. Figure 6 2 Comparison of simulated wall temperature wit h experimental data The characteristic wave length of the oscillation, 2Ls+D is the distance that has to be traveled by a fluid particle when making a single complete internal circulation. Three characteristic time scales can be defined for the heat trans fer of slug flows. The characteristic time for a single complete circulation can be expressed as follows: s s ciru D L 2 (6 16) where us is the slug velocity. The characteristic time for conductive heat tr ansfer inside the solid wall is diff,w=h2s, where, h, s is the wall thermal diffusivity. A similar time scale can be defined for conduction heat transfer for PAGE 112 1 12 the liquid slug, diff,L=D2l. Barid and Mohseni (2008) defined the relative effect o f the internal circulation of the slug as: L D D Pe Ncir L diff2, (6 17) w here Pe=usl is the Peclet number of the liquid. They concluded that a greater value of N provides a larger heat transfer. Thes e parameters will be employed in the following sections of this chapter Results and Discussion Figure 6 2 shows the temperature variation along a microchannel with a length of 6 centimeter. To make sure that a hydrodynamically developed flow was achieved before the heating section, heating started at the middle of the microchannel (3 cm). Heat is applied to the microchannel by a heater whose thickness is 0.25 mm. This configuration is identical to the experimental test rig of Walsh et al (2009) (see T est C ase 1). Figure 6 3 indicates that at the beginning of the thermal entry region, the wall temperature rises rapidly where in the downstream region the flow is thermally fully developed. The radial conduction heat transfer across the wall causes differences between the inner and outer wall temperatures to be noticeable. The small shifting of the temperature difference between the inner and outer walls is due to the axial wall conduction heat transfer. In other words, the outer wall feels the existence of the cold liquid slug with a small time lag compared to the inner wall. The scale of this time lag is diff,w as defined earlier. Figure 6 3 indicates that there are large differences among the gas and liquid temperatures along the axis of the microchannel. Si nce heat can easily PAGE 113 113 penetrate the gas slug through the thin liquid film, the gas slug is always significantly hotter than the liquid slug. Figure 6 3 Temperature rise along the microchannel ( T est Case 1) Figure 6 4 shows velocity vectors and temperature contours for two adjacent liquid and gas slugs in the stationary frame of reference. The length of the liquid slug is used to non dimensionalize the axial position. This figure also shows variations of the dimensi onless heat flux and temperature for the slugs in question. As can be seen, the instantaneous slope of the temperature rise, ( TwTm), changes, but the timeaveraged slope is zero (as it should be) in the fully developed regime (see Figure 6 3). Figure 6 4 shows that the minimum heat transfer occurs at the end of the film region in which a small circulation exists. Generally, axial velocity magnitudes are small in the vicinity of the circulation region compared to the far stream velocity. Circulation occurs due to the reduction of pressure which causes a sharper interface curvature to PAGE 114 114 develop and causes a thinner liquid film to exist. As the film gets thinner, the bubble diameter gets larger. Thus, in order to preserve a constant volumetric flow rate across t he channel, the liquid velocity at the end of film region would need to change direction and causes the circulation discussed above. The low velocity at the end of film region produces a weak heat transfer mechanism in which heat transfer is less than 15% of the average heat transfer (see Figure 6 4 at x/Ls Figure 6 4. Temperature contours, velocity vectors and dimensionless wall heat flux for a train of a liquid and gas slug at the stationary frame of reference (thick line in the first figure indicates gas liquid interface) a b c d PAGE 115 115 Against the expectation of previous researchers who largely ignored the heat transfer in the film region, Figure 6 4 shows that the heat transfer in this region has a contribution to about 25% of the overall heat transfer. This is significant and should not be ignored. During the period in which the gas slug is passing, the wall stores a significant amount of heat due to the small wall to fluid heat transfer. The excess stored heat is released while liquid slug passes. The non uniformity of the heat released from the wall enhances heat transfer in the liquid slug. This result in having more than 150% enhancement in heat transfer in the middle of the liquid slug (see Figure 6 4 at x/Ls In Figure 6 4 other heat transfer peaks are observed in front of the gas slug (Region a). Since the velocity of the gas slug is slightly bigger than the velocity of the liquid slug, the liquid has to be accelerated in Region a thus pumping the cold water into th e liquid film region (in a frame of reference attached to the slugs). A combination of flow acceleration and low flow temperature creates very large peaks of heat transfer. Fluctuations of the heat transfer in Region a and in the film region (Region b ) are mainly caused by strong effects of the surface tension force which creates highfrequency pressure waves. Similar fluctuations were observed in the pressure and velocity fields by Mehdizadeh et al (2009b). One may thus conclude that the fluctuations of the velocity locally enhances mixing and heat transfer similar to what happens in turbulent flow s The nature of velocity fluctuations close to the interface is different than those in traditional turbulent flow. In slug flow, instability of the interfa ce is PAGE 116 116 the source of fluctuations. This means that fluctuations are imposed by the vibrating boundaries even in the low Reynolds number regime. According to the contours of temperature in Figure 6 4, the coldest region is located in the liquid slug right a t the back of the gas slug. It is interesting to note that the maximum heat transfer in the gas slug occurs at the back (tail) of the gas slug too. In order to compare the heat transfer characteristics of the flow with those in the Poiseuille and p lug fl ows it is necessary to introduce the homogeneous void fraction, G L GJ J J (6 18) where JLand JG are the superficial velocities of the liquid and gas phases, respectively. Since the liquid phase is the main contributor toward heat transfer and only has contacts with the wall in (1) portion of time, heat transfer must be normalized by (1 ) in order to have a fair comparison of the heat transfer characteristics of slug flow relative to those in the p lug and Poiseuille flows. As mentioned earlier, heat transfer in slug flow s is inherently unsteady. Therefore two different Nusselt numbers are defined here, an average one and a local one ; The average Nusselt number, Nuave, can be defined as: ) (, m x w ave w aveT T k D q Nu (6 19) where ave wq is the average heat flux provided by the heater. Thi s is easy to measure experimentally. PAGE 117 117 Figure 6 5. Comparison of the averaged and local Nusselt numbers of a slug flow in a microchannel Also x wT, is the wall temperature at any given time and location. Equation ( 6 19) is the only practical way to describe the Nusselt number in the experiments. However, as seen earlier in Figure 6 4, the instantaneous heat flux near the wall has significant fluctuations because of the ability of the solid wall to store and release heat. Thus the real Nusselt number that the wall experiences can be defined by the local Nusselt number, Nux, defined as follows: ) (, m x w x w xT T k D q Nu (6 20) where x wq is the i nstantaneous wall heat flux at any given position. It is obvious that the time averaged value of x wq will be equivalent to ave wq PAGE 118 118 Figure 6 5 displays comparisons of the calculated normal ized averaged Nusselt number, aveNu and the local Nusselt number, xNu with those in the Poiseuille and plug flows. Experimental averaged Nusselt numbers obtained by Walsh et al. (2009) for a sl ug flow (SL/D=1.6) are also shown in Figure 6 5. As can be seen, the Nusselt number decreases rapidly in the thermal entry region. Interestingly enough, the local Nusselt number, xNu can decrease to much lower values than th ose in the Poiseuille flow in the thermal entry region. Figure 6 5 shows that for small values of the dimensionless position, ) Re( slugx plug flow is associated with a thinner thermal boundary layer compared to slug flow. The thicker thermal boundary layer of slug flow causes the averaged Nusselt numbers to have lower values than those in the plug flow. The time scale needed for the wall to feel circulation in the slug flow is proportional to the relative effect of the inte rnal circulation, N which was described earlier. Once the wall feels the circulation of the slugs (x*4 for this length of the slug), the averaged Nusselt number, aveNu will be superior than the Nusselt number for the Plug flow, 1 PlugNu With three examined slug lengths, Wals h et al. (2009) introduced an empirical correlation for the fully developed Nusselt number for slug flows in microchannels as follows: 2 / 1 ,) / ( 5 31S Poiseuille dev slug devL D Nu Nu (6 21) The above correlation shows that the enhancement of heat transfer is inversely proportional to the square root of the slug lengthto diameter ratio. Figure 6 5 indicates PAGE 119 119 that the fully developed Nusselt number can be as high as 27. This value is 610% larger th an that in the traditional Poiseuille flow. One of the significant findings of the research upon which this chapter partly reports is the effect of the solid wall on the overall heat transfer pattern. Figure 6 5 clearly shows that the local Nusselt number, xNu can exceed a value of 100 and drop to a value of 3. This wide range of oscillations is attributed to wall heat storage and release capability. This finding shows a strong dependency of the local Nusselt number on the material and the thickness of the solid wall. In contrast, the average Nusselt number is independent of the wall material and thickness. The effect of the solid wall on the local Nusselt number can be characterized by four dimensionless par ameters. The Biot number (Bi=h(t)Lc/kw), a dimensionless length (L*wt/L*2) and the dimensionless heat generation (g*= g wt/(kw not practical to investigate the effect of solid wall on the tem perature field of slug flows because of the complexity of the problem, huge computational cost and the need to account for the numerous free parameters that characterize the problem. Therefore, conduction in the solid wall will be studied separately by imp osing a time dependent boundary condition for the convective heat transfer coefficient. The latter can be obtained from the results of the simulation. Equation ( 6 22) shows the governing energy equation in cylindrical coordinates. The axial conduction in t he solid wall is neglected for purposes of simplifying the model a bit. dt dT k g dr dT r dr T d1 12 2 (6 22) PAGE 120 120 In the above equation T denotes the temperature, g generation rate and thermal diffusivity of the wall, respectively. The outside of the heater is assumed to be well insulated as Walsh et al. (2009) reported in their study. The boundary conditions for the above equation can be written as: 0 Ro rr T k (6 23) ) )( (m R r w R rT T t h r T k (6 24) where h(t) is the time dependent convective heat transfer coefficient on the interior wall of the microchannel. Since the wall experiences the same trend when slugs are passing, a periodic condition can be introduced for the temperature field in the solid wall as follows: ) ( ) ( t r T t r T (6 25) where, is the time needed for two adjacent liq uid and gas slugs to pass from any specific location. The above equation can be analytically solved by Greens theorem in integral form. A finite difference code was developed to calculate the temperature field inside of the solid heater. The timedependen t convective heat transfer coefficient can be found by interpolating a curve on the local heat transfer coefficient obtained from numerical simulations for Test Case 1. Figure 6 6 (right) shows the convective heat transfer coefficient of slug flow al ong the microchannel. The bold black line with square symbols shows a sample variation of this coefficient over one period of time, By knowing the slug velocities, the time dependent convective heat transfer coefficient h(t*), or al ternatively the time dependent PAGE 121 121 Nusselt number, Nu(t*), can be introduced by a polynomial curve fit. The error of curve fitting in Equation ( 6 26) was less than 10%. 1 t if 1) Nu(t ) Nu(t 1 t 0 if 1.2238 19.232t 194.85t 936.54t 2114.9t 2131.8t 777.91t ) (* * 2 3 4 5 6 aveNu t Nu (6 26) where t* is defined as /*t t The above equation is valid for one period, where t< (=(Lg+LL)/ v and where v is the slug velocity). This is repeated in subsequent periods. Figure 66 Wall temperature variation of the solid heater for several times, left. Convective heat transfer coefficient in different locations, right Figure 6 6 (left) shows dimensional variation of wall temperature in several nondimensional times, t*=t/ and positions, r*=r/D. It can be seen that when the gas slug is passing, the heat transfer coefficient decreases which results in a small Biot number (Bi). In such a case, heat penetration happens significantly faster than the heat removal by convection. This makes for almost a flat temperature profile inside of the wall, as can be seen in the time range 0.1 PAGE 122 122 average temperature becomes smaller. Figure 6 6 clearly shows the effect of heat storage and heat release by the solid wall on the local convective heat transfer coefficient. If the steady state solution of Equation ( 6 22) is used for nondimensi onalization, the temperature profiles would be smoother. Therefore, the dimensionless temperature T* may be defined as: ) ( ) ( ) ( ) (, *t T t r T r T t r T Tss w ss (6 27) where Tss is the steady state solution for different values of the convective h eat transfer coefficient, while Tw,ss is the steady state inside wall temperature for different values of the convective heat transfer h(t) The solution of the steady state version of Equation ( 6 22) can be expressed as: 2 1 24 c Lnr c k r g Ts ss (6 28) where the constants are: ) 2 4 2 2 ( 22 2 2 2 2 1 s o s o m s ok LnR R k R Rh R h R g T c and k gR c (6 29) Figure 6 7 shows the dimensionless temperature profiles, T*, in the solid wall for various dimensionless times, t*(=t/ ). The figur e clearly shows that for small Biot numbers (Bi<0.05), the dimensionless temperature profiles are uniform and their values approach unity. PAGE 123 123 Figure 67 Dimensionless temperature in the solid wall for various dimensionless times The observed uniform tempe rature profile is a sign of faster heat penetration compared to the rate of heat removal caused by the convective heat transfer on the inner wall. During the time in which the Biot numbers are small, the wall stores heat. The wall will release heat when th e Biot numbers are larger than the average Biot num b er (Bi>Biaverage), and depending on the slope of the convective heat transfer coefficient and the amount of the stored heat, the values of T* may become larger or smaller than one. The above findi ng can have many practical applications especially in places where heat is generated periodically (such as in CPU or Microprocessors). By controlling the slug lengths and velocity, spots of high local Nusselt numbers (as was observed earlier in Figure 6 5) can be matched with the spots of intense heat release. Conclusions In this study, slug flow in a microchannel with a constant heat flux boundary has been simulated using the volume of fluid (VOF) method. While the liquid film around the PAGE 124 124 bubble was not successfully captured in many previous numerical studies, mainly due to employing course meshes to minimize computational time, in this study the liquid film was successfully captured by a new strategy of mesh generation. Heat transfer characteristics of the slug flow were compared to those in the Poiseuille and plug flows. It was found that the combination of gas slug and liquid film contribute as much as 20% to the overall heat transfer. A strong coupling between the conductive heat transfer in the s olid wall and the convective heat transfer in the flow field was observed. Local and time averaged Nusselt numbers for slug flows were reported for the first time. Significant enhancement of heat transfer was observed in short slugs where the Nusselt numbe r was found to be 610% higher than that in Poiseuille flow. One of the unique findings of this study is the considerable oscillation of the local Nusselt number where the latter was found to vary from 3 to 100 depending on the location and length of the sl ugs. Results of this study can be applied to microscale chemical and mechanical processes. PAGE 125 125 CHAPTER 7 APPLICATIONS OF TWO PHASE SLUG FLOW IN MICROCHANNELS TO ESTIMATING INTERNAL LEAKAGE IN ROTARY VANE TURBOMACH INERY Outlook Understanding the physic s of two phase flows in micro scale leakage paths in any devices is one of the applications of the study reported in the previous chapters. Rotary machines have played an important role for many years in refrigeration and air compression applications because of their inherent simplicity and reliability. They are also very attractive machines since as positive displacement devices; they are more suitable for low flow rates (low specific speeds). In this chapter the thermodynamic and fluid mechanic characte ristics of a rotary vane air motor are analyzed. The optimum geometrical and operational characteristics of the machine are presented. Experiments are conducted to understand the working principles and operational constraints of the machine. This study als o helps formulate design procedures that can be utilized to modify air motors into optimized expanders for singlephase flow applications. The model has been used to evaluate geometrical parameters such as the optimum intake and exhaust port locations, their spreads and the geometric volume ratio, as well as evaluating performance parameters such as the work produced and the mechanical, isentropic and total efficiencies of the machine. This chapter demonstrates that a large internal leakage that inherently exists in rotary vane expanders undermines its functionality. We have used models developed in the previous chapters to model the internal leakages in a rotary vane expander. PAGE 126 126 It is anticipated in a follow up study that the model developed will be the basis for an expander design tool that uses twophase working fluids in relevant industrial applications. General Background Refrigeration Cycle The purpose of this and the following sections is first to explain a fundamental concept of refrigeration cycles and second to explain the feasibility of COP incensement by installing a rotary vane expander. Refrigeration is the process in which heat is removed from an enclosed space, or from a substance Thermodynamic Overview on Refrigeration Cycle The vapor compre ssion refrigeration cycle is used in most household refrigerators as well as in many large industrial refrigeration and HVAC systems. Figure 7 1 provides a schematic diagram illustrating the various components of a typical refrigeration system. Figure 7 1 A schematic of a typical refrigeration system PAGE 127 127 The thermodynamic processes of the ideal cycle are shown in Figure 7 2 on a T s diagram. In this cycle, a circulating refrigerant (working fluid) enters the compressor as a vapor at point 1. From point 1 to point 2, the vapor is compressed isentropically and exits the compressor as superheated vapor. From point 2 to point 3 and on to point 4, the superheated vapor travels throughout the condenser which first cools the superheated vapor and then condenses the saturated vapor into saturated liquid by removing additional heat. Since the pressure drop in condenser is negligible, ideally the condensation is an isobaric process. Between points 4 and 5, the saturated liquid refrigerant goes through the expansion valve (also called a throttle valve) or capillary tube where its pressure suddenly drops, causing flash evaporation, typically, for less than half of the liquid. Figure 7 2 Temperature vs. Entropy diagram for a typical refrigeration cycle That results in a liquidvapor mixture at a lower pressure and temperature as shown at point 5. The cold mixture of liquid and vapor then travels through the evaporator and is completely vaporized by getting heat from the lower temperature PAGE 128 128 space (from the space being refrigerated). The heated vapor returns to the compressor inlet at point 1 to complete the thermodynamic cycle. In the actual vapor compression refrigeration cycle there are many effects like frictional pressure drop in the system and slight thermodynamic irreversibility in the compression process that was not considered in the ideal cycle. The coefficient of performance for cooling, of refrigeration cycle is defined as: 1 2 5 1h h h h W Q COPcold ( 7 1) w here coldQ is the absorbed heat from the space to be cooled and W is the work done by the compressor. Therefore the COP is the ratio of the heat removed from the cold space to the input work. Potential of Increasing Coefficient of Performance The higher the COP rating of a refrigeration unit, the higher its energy efficiency. There are many ways to increase the COP of a system for a set of low and high temper atures. Any irreversibility in the cycle can reduce the COP. One of the main irreversibility sources is the flash evaporation process in the throttle valve. In the throttling process high pressure liquid refrigerant becomes a low pressure liquidvapor mixt ure while it doesnt do any work. In the ideal cycle, isenthalpic process (process 4 to 5 in Figure 7 2 ) can be replaced by an isentropic process in which flow expands isentropically that is shown by dashed line 4 5 in Figure 7 2 The ideal COP for a refr igeration system employing ideal rotary vane expander can be defined as PAGE 129 129 ) ( ) (' 5 4 1 2 5 1 exph h h h h h W W Q COPander c cold exander (7 2) The performance improvement of a refrigeration system after installation of a rotary vane expander is extensively studied in the PhD dissertation of Mahmoud (2008). Figure 7 3 Various components of a rotary vane expander machine Figure 7 3 shows various components of a rotary vane machine. The geometry is based on two eccentric cylinders and sliding vanes positioned in the slots of the rotor This study reveals that the increase of system COP by installing a rotary vane expander is possible, but low efficiency is a strong barrier to its functionality. It is found that the losses in this device are mainly caused by large leakage through inter nal leakage paths. A large internal leakage in rotary vane expander motivated us to start an extensive fundamental study on twophase flows in microchannels which was reported in the previous chapters Now after all those extensive studies we are trying to employ the knowledge we obtained from those studies to estimation the internal leakage in a rotary vane machine. PAGE 130 130 Microscale Leakage Paths in Rotary Vane Expander Figure 7 4 shows a schematic of two phase flow through a microscale leakage path exist at the tip of a vane. Figure 7 4 Schematic of two phase flow passing through a microscale leakage path at the tip of a vane Understanding of the physics of twophase flow in leakage paths is one (out of thousands) example of the applications of the rep orted study. In this chapter a physical description of the combined phenomena of leakage, friction and irreversibility due to the processes that occur within a rotary vane motor are formulated. Using these results, the potential of improving the design of the device with respect to these losses is discussed. This study is a contribution to ongoing research investigating the possibility of replacing conventional capillary tubes or throttle valves with rotary vane expanders in refrigeration systems. Expansi on in a conventional refrigeration system occurs by means of an isenthalpic process The potential to use a rotary vane expander to carry out expansion by means of an isentropic process will consequently increase the coefficient of performance of the cycle. Before doing so, modifications to a commercial air motor are necessary. PAGE 131 131 Analytical Model (a) ( b ) (c) Fig ure 7 5 The tested rotary vane air motor (a and b) and the s chematic of the device Geometr ical Analysis The rotary vane air motor used to collect the experimental data is shown in Figure 7 5 ( a ) and 7 5 (b ) The schematic device is also shown in Figure 7 5 (c ). The PAGE 132 132 commercial device used for the purpose of this study is a Gast Corporation air motor model NL32NCC1. The geometry is based on two eccentric cylinders and sliding vanes positioned in the slots of the rotor. When the rotor rot ates, the inner and outer cylinders and the side vanes create a chamber that can expand or compress according to its position. Inlet flow through the admission arc is admitted into the rotor slot cavity beneath the vane and forces the vanes out. The fluid then enters the cavity through a pressurized slot represented by a dotted line in Figure 7 5 (c ). If the center of a polar coordinate system is set at the center of the rotor, an equation of the outer circle, the stator, can be written as: 2 2) cos ( sin d R d ro ( 7 3 ) In which d is the eccentricity, Ro is the stators radius and is the angle shown in Figure 7 5 (c ). If the number of vanes, n, is known then the area of eac h chamber can be defined as d R r An i 2 2 21 1) ( 2 1 ) ( (7 4) In which Ri is the rotor radius, 1 is the angle of a particular vane from the origin and is the angle corresponding to the thic kness of a vane. Since the majority of a vane is located far from the center, it can be assumed that does not change during rotation. An average value of is found to be 0.157 Radian. This integration can be evalu ated analytically or numerically. The cross sectional area of one chamber and the variation of the chamber volume as a function of can be estimated as shown in Figure PAGE 133 133 7 6 The focus of the current study is for an air motor in which no expansion occurs. In this case, a chamber that is disconnected from the inlet is immediately connected to the exhaust port (Figure 7 5 ). Since the pressure differences between the chambers of an air motor are much greater than that in an expander, a higher pressure gradient is the leakage potential. The exhaust ports location is selected so that the exhaust process begins when the volume between adjacent vanes (cavity) is maximum. The angle of the maximum volume is a strong function of the number of vanes and is independent of the eccentricity. The optimum angular location for the inlet and outlet ports is obtained by differentiating Equation (7 2 ) with respect to and setting it equal to zero. 0 ] cos sin ) 2 ( cos ) 2 [sin( 2 2 cos )) 2 ( 2 cos(2 2 2 2 k n k n n (7 5 ) Whe re oR d k For a four vane air motor, the angles 2250 and 45o correspond to the maximum and minimum chamber volume positions, respectively. The exhaust port arc has a spread of 135o that starts from 225o (location of maximum volume) to the minimum volume. Ideally, this will eliminate compression losses. If the vane thickness is neglected the intake port should start from 45o (minimum volume) to 45o corresponding to the maximum volume. Spread angles of 90o and 135o for the admission and exhaust arcs allow the air motor to intake and discharge air freely. Robertson and Wolgemuth (1978) showed that the chamber pressure at the inlet and the exhaust ports corresponds to the supply and exhaust pressures, respectively. PAGE 134 134 Fig ure 7 6 The calculated v ariation of cross sectional area of a chamber during one revolution ( m2) Vane Friction Friction forces are determined by utilizing a free body diagram and calculating the forces acting on a vane. Since vanes can freely slide in their slot s the clearance between the end plates and vane sides is small and the frictional losses due to the vane sides are neglected. An enlarged schematic of a vane and the acting frictional and pressure forces, weight and centrifugal forces are presented in Figure 7 7 The Coriolis acceleration is neglected. After solving the force and torque balance equations for these 12 forces, N Nb and Nt are determined. The normal force, acting on the tip of the vane is written as ] 3 ) sin 2 ( 2 ) 5 ) ( 4 / 2 ( [ 1 1 Wx o P c g r m c LW c x L W c n hW c i P c N (7 6 ) In th e above equation, x can be found as a function of PAGE 135 135 Figure 7 7 Schematic of forces acting on a vane (not to scale ) The constants in this equation are L x h L h c ) 2 / ( 2 ) 1 (2 1 ) (2L x h L x c ) ) 2 / ( 2 (3 L x h x L c ) ) ( (4 L x h x L c ) (5 L x h L c w here W is the width of chamber, is vane stator coefficient of friction and L is the length of the vane. The parameter n in Equation ( 7 6 ) is introduced since at higher rotational speeds the large friction force ( N ) at the tip of the vanes quickly erodes the graphite vanes leaving the tips edges dull. Leakage past the tip of the vanes cannot be estimated. If there is leakage, then a significant pressure gradient exists across the tip of the vane. If this is the case a correction factor for the vertical pressure force must be introduced to account for the average pressure force at the tip. PAGE 136 136 Figure 7 8 shows the variation of normal forces N, Nb and Nt during one rotation at 4.2 bar inlet pressure. Fig ure 7 8 Variation of normal forces acting on a vane as a function of (deg) This figure shows that the normal force at the tip of the vane, N is approximately constant and a significant ju mp occurs when pressure loading begins. It should be noted that the centrifugal forces generated are substantially larger than the weight force since N is nearly constant during periods of noload. During the inlet, the side friction, Fft and Ffb (Figure 7 7 ), help decrease N and consequently decrease the frictional force at the tip of the vane. This decreases the possibility of leakage past the vane tip. It should be examined whether lighter vanes with a larger coefficient of friction can be utilized to ac hieve smaller frictional loss during the period of no load. This may lead to a more efficient device. PAGE 137 137 Torque a nd Power Torque is produced as the result of the forces acting on vanes. The torque produced by a vane is defined as ) ( L x R N R N Ti t i b v (7 7 ) Figure 7 9 shows the torque produced by a single vane during one revolution. It is seen that due to friction there are negative torques generated during a large portion of one revolution. When the vane is exposed to the high pressure region, both vane friction and torque increase. Significant pressure forces on the vane surface generate significant torque that is maximum at 270o. This angular location corresponds to the largest protruded area of the vane, x Vane f riction is not the only source of frictional loss. In practice, after each run, when the device is opened, uniform rubbing is noticed on both end plates. For sealing purposes the side clearances are so tight that rubbing is inevitable. The frictional force s increase with increased speed due to larger centrifugal forces. It is assumed that the sides frictional forces are not a function of speed since the tolerances do not change with speed. The side frictional forces are solely a function of the torque that is spent to tighten the screws that hold the end plates and stator cylinder together. Before each test run, an identical amount of torque is used to close the device. There is no guarantee, however, that the side friction is identical in each run. They are also a function of time, pressure and other unknown parameters. For example, at the higher pressures, the strain causes the side end plates to be deformed. This consequently decreases the sidewall frictional forces. However, the only thing that is PAGE 138 138 certain during each run is that the torque lost due to side friction is constant. Ideally this practical constant, k should be a function of the coefficient of friction. Figure 7 9 Single vane torque during a revolution In reality, however, the rotor and side end plates slightly rub against each other. The standard friction law ( N ) is hence not valid and hence this constant must be found by empirical means. The total torque and power due to the contribution o f each vane and side fricti onal loss is expressed by i vR k d T n Torque 2 02 ( 7 8 ) Torque RPM Power 60 2 ( 7 9 ) PAGE 139 139 Figure 7 10 Measured and the calculated t orque and power for 2.8 bar inlet pressure It can be seen that the torque decreases at higher speeds since centrifugal forces impose a larger normal force and hence a larger friction force. The analytical estimate of power shown in Figure 7 10 is in good agreement with the experi mental data. Ideally, power should increase linearly with speed but frictional losses actually reduce it. This leads to a maximum value for power after which the power decreases with increased speed. The maximum error of the analytical model is approximately 10%. In addition, the speed corresponding to the maximum power is efficiently predicted by this model. PAGE 140 140 To understand the frictional losses in better detail, Figure 7 10 shows the two sources of losses separately. The difference between ideal and no v ane friction in Figure 7 10 is representative of the side frictional loss. Moreover, the difference between the no vane friction and analytical curves represents the vane frictional loss. It seems that both of these losses have the same order of magni tude. Ideal Volumetric Flow Rate The ideal volumetric flow rate is defined as the sweeping volumes per revolution. The outlet volumetric flow rate is defined as a o cslot c slot c o i iV V V V P P RPM n Q ) ( ) ( ) 60 ( (7 10) w here c is the cut off angle, o is the discharge angle and a is the arc of the exhaust port. The parameter n is the number of vanes. cV is the volume of cavity at the corresponding cut off angle. slotV is the vol ume of the small cavity underneath the vanes (a function of x ) in addition to the constant volume of the pressurized slot (Figure 7 5 ). Since the volumetric flow rate is not conserved and the outlet volumetric flow rate should be determined, the first term in the bracket should be multiplied by a density ratio which is proportional to the pressure ratio for an ideal gas at constant temperature. It is obvious that the ideal flow rate should equal zero at zero rotational speed and increase linearly with speed However, experimental data shows significant leakage during the stationary (zero RPM) and highspeed modes of operation. The experimental volumetric flow rate also increases linearly with speed but with a lower slope than that of the analytical model. The difference between the experimental data and analytical PAGE 141 141 solution represents the amount of leakage through the device. It can be seen that the leakage decreases at higher speeds. Leakage The dimensions of the leakage paths in the stationary mode can be accurately measured. There are four possible leakage paths. There exists a leakage path due to clearances between the vanes and the slots, vane sides and end plates, rotor side and end plates and between the stator and rotor cylinders (see Figure 7 5 ). Li ke all turbomachines there is no simple way to measure these clearances during operation. As an example, at 4.2 bar, the inlet pressure force is approximately 1035N ( area ojected P Pr ) acting on the rotors shaft to the left hand side (see Figure 7 5 ). However these side effects of strain on the size of the top clearance or on the axial gap between the rotor and end plates cannot be predicted. During operation, the location of the vanes impacts the leakage flow rate. This effect can only be predicted during the stationary mode of operation and is neglected The greatest difference between the stationary and dynamic modes of operation is the orientation of the vanes. I t is assumed that the vanes are always blocking the path between the vanes and the sl ots. Figure 7 7 shows that this assumption is valid since a significant pressure force on one side of the vane and/or the frictional force at the tip of the vane rotate the vane in the slot and seal this path. The remaining leakage paths can be modeled as a Couette type flow. A schematic of a leakage path is shown in Figure 7 12 For all leakage paths, it should also be checked if the flow is laminar, turbulent and choked. PAGE 142 142 Since there is a significant pressure gradient in each path, the constant density as sumption is not valid. Even though the pressure gradient may not be linear in each path, a linear assumption does not produce erroneous results since the lengths of the various leakage paths are very small. To find the velocity profile in each path, the co ntinuity equation for an ideal gas is solved. The first law formulated for a control volume is solved assuming that the process through the leakage paths is isenthalpic. Hence the temperature is assumed not to vary significantly. The flow is also assumed q uasi developed in each section and hence the y direction component of velocity remains zero (see Figure 712) Fig ure 7 11 Measured and calculated (ideal) v olumetric flow rate at the outlet at Pi=2.8 bar The governing equations of the flow in the leak age paths are: PAGE 143 143 iP x l p P (7 11 ) RT P (7 12) 0 ) ( x u (7 13 ) In Equation ( 7 11 ), l is the leng th of the leakage path. Solving these equations simultaneously, the velocity profile is expressed as ) ( ) ( y g P x l P P y x ui i (7 14 ) Fig ure 7 12 Schematic of leakage paths The x momentum equation is a nonlinear equati on that cannot be solved analytically. To find g(y), however, the convective term in the x momentum equation is neglected Physically it has a small effect on flow field as verified by a scale analysis. The velocity profile is hence PAGE 144 144 y a U ay y l P P x l P P y x ui i 0 2) ( 2 1 ) ( (7 15 ) w here a is the gap between the plates and U0 is the velocity of the moving wall. For any gap Uo can be found as a function of the radial position multiplied by the angular velocity ( 60 / 20RPM r U ). The volumetric flow ra te can then be found by integrating the velocity profile as : 312 60 a l P a RPM r P x l P P W A U Qi i ave l (7 16) Since the rotor is symmetric, it can be assumed that if one vane at the bottom has a velocity of Uo another vane is located at the other side (exactly 180o apart) has a velocity of Uo Because of this, the influence of rotational speed is assumed negligible. The side leakage is found to be a strong function of the pressure ratio and is also independent of the rotational speed. The largest source of leakag e corresponds to the minimum clearance between the rotor and stator cylinders. Since the flow at that point can be modeled as that in a convergingdiverging nozzle, the volumetric flow rate can be expressed as o i i k k LT P A k R k Q/ ) / )] 1 /( 2 [ ) / ( () 1 ( 2 ) 1 ( (7 1 7 ) w here A is the cross sectional area at the throat. The rotor stator gap size is difficult to measure during operation. It is evaluated during stationary operation and is estimated to be 0.14 mm. Table 7 1 shows the variation of QL for different pressures Table. 7 1. Leakage flow rate from the clearance between the rotor stator cylinders for different pressures Pressure (bar) 1.4 2.8 4.2 5.6 QL (L/sec) 1.71 3.42 5.14 6.85 PAGE 145 145 The total leakage can be seen in Figure 7 13 which shows that the analytical model developed can predict the leakage and volumetric flow rate reasonably well. At high speeds, however, the actual leakage decreases. This does not agree well with the model. One of the reasons for this is the fact that practically, leak age is not a steady state process. In general, a vane is located between two highand low pressure sources. The flow must accelerate in this leakage path. Because of momentum, this acceleration needs time. At higher frequencies (higher speeds) the flow has a shorter time to accelerate. This leads to a smaller average velocity in the leakage path. Consequently the leakage flow rate decreases. Fig ure 7 13 The Measured and the modeled volumetric flow rate and leakage (Pi=2.8 bar) This conclusion can als o be shown by a scale analysis. For the unsteady Navier Stokes (N.S) equations, the scale of the average velocity in the leakage path can be PAGE 146 146 expressed as ) /( l RPM P Uave At higher speeds, the momentum of the flow does not allow it to completely acc elerate. At low speeds, in an ideal device, we do not expect there to be any flow. In reality a significant amount of flow exists. The contribution of this research lends itself to determining the weaknesses, performance losses and the physics of leakage w ith single phase fluids as a first step in modifying the device. The experimental data have been fit and non dimensionalized with respect to the pressure ratio. This is best represented by ) ( ) (2 1 exp o i n o iP P c RPM P P c Q (7 18) w here Qexp is th e experimental flow rate (L/sec), 3 110 07 1 c 36 12 c and 65 0 n Fig ure 7 14 The calculated l eakage ratio for different inlet pressure s PAGE 147 147 The second term on the right hand side of Equation ( 7 18 ) represents the l eakage during the stationary mode [ Equation ( 7 17 )] The first term represents the combination of the polytropic process for free expansion and a density ratio that is proportional to the pressure ratio. That density ratio factors in when converting the ma ss flow rate to a volumetric flow rate. The leakage ratio is defined as the ratio of actual volumetric flow rate to the ideal flow rate (analytical without leakage): i RQ Q Lexp (7 19 ) Figure 7 14 shows the variation of the leakage ratio as a function of RPM for various inlet pressures. It is obvious that at zero RPM, at which the ideal flow rate is zero, all the curves asymptote to infinity. As the speed increases, the curves asymptote to 1 or deviate a small distance fr om unity due to the different slopes of the Qexp and Qi. Since the portion of the leakage to the total flow rate is smaller at higher flow rates, the analytical flow rate approaches the actual flow rate. It can also be seen that as the pressure increases, the curves tend to be closer to each other and the effect of leakage becomes smaller. Efficiencies We can define three different efficiencies for this device. friction noi actual i mw m w m (7 20 ) s w i m actual w i m s (7 21 ) PAGE 148 148 s w l m i m actual w i m ) ( (7 22 ) The mechanical efficiency, m reflects the effects of mechanical frictional losses in the device. The second law efficiency, s repres ents how close this device operates relative to an isentropic machine. The overall efficiency represents the actual (total) efficiency of the device by which leakage and friction are taken into account. Figure 7 15 shows the behavi or of the different efficiencies at a 2.8 bar inlet pressure. Since at zero speed there is start up friction the mechanical efficiency is less than 100% and decreases as the frictional losses increase with higher rotational speeds. If this device works ideally, i.e. no leakage, the isentropic efficiency will decrease as the rotational speed increases. If we compare it with the mechanical efficiency curve, we can conclude that at low speeds, frictional losses are small and the process is not ideal. In an ai rmotor the pressure suddenly drops from the inlet pressure to the ambient pressure and creates significant irreversibilities. If the device is modified to expand the working fluid, we may achieve a process that resembles an isentropic process much more cl osely. Due to significant leakage, however, this goal is not easily achieved. Some method by which the leakage paths may be reduced without dramatically increasing the friction is needed. The distance of the mechanical efficiency curve with the 100% effi ciency line represents the frictional losses. The distance between the isentropic efficiency and mechanical efficiency curves represents free expansion losses (irreversibility due to free PAGE 149 149 expansion). At lower speeds, the irreversibility of free expansion i s much more important than frictional losses and vice versa. In Figure 7 15 the total efficiency, with leakage considered, is represented as a solid line. At zero speed when there is no work, significant leakage flow escapes the device without doing any work. This leads to zero efficiency. The maximum efficiency is about 30% and remains almost constant for a wide range of speeds (1200 to 2300 RPM). Fig ure 7 15 The calculated e fficienc ies at Pi=2.8 bar Even though the maximum power occurs at 2500 RPM the maximum efficiency occurs at 1600 RPM and decreases slightly at higher speeds. There exist significant potential modifications to seal the leakage paths in order to maximize the isentropic efficiency. For example, at a speed of 1500 RPM, there are PAGE 150 150 35 %, 12% and 27% potential for efficiency improvement due to frictional losses, free expansion loss and leakage losses, respectively. However, these parameters also depend on each other and any change in one parameter greatly affects the others. To reduce fr ictional losses it may be possible to use lighter materials for the vanes with a smaller solid solid frictional coefficient. The same can be done for the end plates. On the other hand modifications to reduce leakage increase friction. Irreversibilities due to free expansion may be overcome by allowing the flow to expand in the chamber. This modification is easily achieved by implementing an earlier cut off angle (Fig ure 7 5 ) that allows the flow to expand and produce more work. Ideally, the isentropic effic iency should increase and match the mechanical efficiency. This however, cannot be done until leakage is reduced. Small leakage flows from one chamber to the next (high pressure to low pressure) during the expansion process also affect expansion. The most complicated and important modification would be to reduce leakage in such a way that does not increase friction. This is quite difficult, since with any reduction of leakage, friction inherently increases. It is shown that efficiency is more sensitive to f riction than leakage at higher speeds. In conclusion, modifications made to this device are a challenging optimization problem that has many known and unknown variables. Two Phase Internal Leakage Estimation Internal leakage of two phase working fluid in the studied rotary vane turbomachine can be estimated by the correlations developed and presented in the previous chapters. The geometry of the single phase leakage passing across the side of the vanes was shown in Figure 712. For two phase working fl uid, the geometry is the same but PAGE 151 151 flow most likely will behave as the slug flow. If a rotary vane turbomachine replace with a throttle valve in a refrigeration cycle, a high pressure single phase working fluid will enter to the leakage path and due to the pressure drop a portion of the liquid will evaporate and forms a slug flow. For the converging section of the flow going to the leakage path pressure drop can be simply calculate by Equations (5 19 ), ( 5 20 ) and ( 5 21 ). Pressure drop of twophase slug flow entering the micro scale leakage path can be calculated by one of the correlations presented in Table 4 4 Since the leakage paths are not circular the hydrodynamic diameter is used for modeling. For the limiting case were the clearance is much smaller t han the length of the vane, hydraulic diameter can be estimated as two times of the clearance (DH = 2 a ). Table 7 2. Leakage flow rate (mean) from the clearance between the side of the vanes and end plates Pressure (bar) 1.4 2.8 4.2 5.6 Estimated velocit y (m/s) 4.41 9.39 12.53 15.04 QL, (L/sec .m ) 0.88 1.87 2.50 3.00 Table 7 2 shows the calculated side leakage per unit length of vane for the several inlet pressures. If the leakage from the top clearance of rotor and stator adds up to the above flow rate values, the total two phase int ernal leakage in rotary vane turbomachine can be estimated. Conclusion In this chapter the thermodynamic and fluid mechanic characteristics of a rotary vane air motor were analyzed. The optimum geometrical and operational characteristics of the machine are presented. The model has been used to evaluate geometrical PAGE 152 152 parameters such as the optimum intake and exhaust port locations, their spreads and the geometric volume ratio, as well as performance parameters such as the work produced and the mechanical, i sentropic and total efficiencies of the machine. Results show that there are three main issues that affect the performance of this device; frictional losses, leakage, and irreversibility due to free expansion. Results show quantitatively that at low speeds leakage dramatically suppresses the efficiency while at the high speeds frictional loss reduces the performance. PAGE 153 153 CHAPTER 8 RECOMMENDATIONS AND FINAL CONCLUSIONS In this dissertation many concepts and physical investigations dealing with the two phase slug flow s in microchannels were described. Below is a list of the contributions of this study: The effect of gas liquid interface curvature on the calculated pressur e drop was investigated. It is found that the surface tension effect which produces t he gas liquid interface curvature decreases the calculated pressure drop values especially for the short slugs A simple correlation capable of predicting the Moody fricti on factor for slug flows has been introduced. Guidelines to set the refinement thresholds for dynamic mesh adaption were developed. The negative effects of using a large elemental aspect ratio have been described. An advance mesh generation method was dev eloped and explained to be able to employ the minimum number of elements to capture liquid film successfully. The former was successfully captured employing 96% less mesh compared to previous numerical investigations. A novel model capable of predicting tw o phase slug flow void fraction in vicinity of a flow area change was introduced. An analytical model was developed for the pressure drop of two phase slug flow through a sudden area change. Results showed excellent progress in pressure drop prediction. This new model could reduce the minimum (60%) and maximum (500%) errors of the more conventional models to a minimum of 8% and a maximum of 310%. PAGE 154 154 For the first time, this study reveals that the liquid film region contribute s as much as 20% to the overall heat transfer. A strong coupling between the conductive heat transfer in the solid wall and the convective heat transfe r in the flow field was found. Local and timeaveraged Nusselt numbers for slug flows were reported for the first time. Significant enhancement of heat transfer was observed in short slugs where the Nusselt number was found to be 610% higher than that in Poiseuille flow. One of the unique findings of this study is the existence of the considerable oscillation of the local Nusselt number where the local Nusselt number was found to vary from 3 to 100 depending on the location and the length of the slugs. Results of the study on a rotary vane turbomachine showed that there are three main issues that affect the performance of this device; fric tional losses, leakage, and irreversibility due to free expansion. Results show quantitatively that at low speeds leakage dramatically suppresses the efficiency while at the high speeds frictional loss reduces the performance. One of the main challenges o f this study was the cost of the numerical simulation. As was mentioned earlier the computational time of some cases was more than 3 months. It seems that employing parallel processing might decrease the needed computational time. Another difficulty we had was the lack of the data from experiment. For example slug lengths were not measured in most of the experiments while it has a significant effect on the overall heat transfer and pressure drop. Also the local Nusselt number which was reported in Chapter 6 has never been measured experimentally before. PAGE 155 155 Results of this study can be applied in many microscale chemical and mechanical processes. PAGE 156 156 LIST OF REFERENCES Abdelall, F.F., Hahn, G., Ghiaasiaan, S.M., Abdel Khalik, S.I., Jeter, S.S., Yoda, M., and Sadowski, D.L., (2005) Pressure Drop Caused by Abrupt Flow Area Changes in Small Channels, Experimental Thermal and Fluid Science Vol 29, No. 4, 425 434. Angeli, P., and Gavriilidis, A., (2008) Hydrodynamics of Taylor Flow In Small Channels, IMECE2008, Part C: Mechanical Engineering Science Vol. 222, 737 751. Armand, A.A. and Treschev, G.G., (1946) The Resistance during the Movement of a Two Phase System in Horizontal Pipes, Izv. Vses. Teploteck. Inst. Vol 1, 16 23 (AERELib/ Trans 828). ASHRAE, (2001) ASHRAE Handbook: Fundamentals (Chapter 2), the American Society of Heating, Refrigerating and Air Conditioning Engineers, Inc., Atlanta, GA. Aussillous, P., and Quere, D., (2000) Quick Deposition of a Fluid on The Wall of a Tube, Physics of Fluids Vol 12 (10), 2367 2371. Baird, E. and Mohseni, K., (2008) Digitized Heat Transfer: A New Paradigm for Thermal Management of Compact Micro Systems, IEEE Transaction on Components and Packing Technologies Vol 31, No. 1. Baker, O., (1954) Simultaneous Flow of Oil and Gas, J. Oil Gas. Vol 53, 184 95. Black, K., Bruno, A.B., and Martel, J., (2007) A Preliminary PIV Analytical Investigation of Wall Shear in Micro Channel Slug Flow, IMECE2007, Seattle, Washington, U SA. Brackbill, J.U., Kothe, D.B., and Zemach, C., (1992) A Continuum Method for Modeling Surface Tension, J. of Computational Physics Vol 100, No. 2, 335354. Brauner, N., (1990) The Relation between Two Phase Flow under Reduced Gravity and Earth Ex periments, Int. Communications Heat Mass Transfer, Vol. 17, 271 282. Bretherton, F.P., (1961) The Motion of Long Bubbles in Tubes, J. Fluid Mechanics Vol 10, 166188. Chalfi, T.Y. and Ghiaasiaan, S.M., (2008) Pressure Drop Caused by Flow A rea Changes in Capillaries Under Low Flow Conditions, Int. J. Multiphase Flow Vol 34, 2 12. Chisholm, D., (1983) TwoPhase Flow in Pipelines and Heat Exchangers, Pitman Press, Bath, England. PAGE 157 157 Collier, J.G., (1972) Convective Boiling and Condensati on, McGraw Hill Book Co, UK. Damianides, C.A. and Westwater, J.W., (1988) Two Phase Flow Patterns in a Compact Heat Exchanger and Small Tubes, in: Proceedings of Second UK National Conference on Heat Transfer, Glasgow, September 1416. Mechanical Engin eering Publications, London, 12571268. Darhuber, A.A., and Troian, S.M., (2005) Principles of Micro Fluidic Actuation by Modulation of Surface Stressesm, Annual Review of Fluid Mechanics Vol 37, 425455. Darvish, M., and Moukalled, F., (2006) Conv ective Schemes for Capturing Interfaces of Free Surface Flows on Unstructured Grids, Numerical Heat Transfer, part B, Vol. 49, 1942 De Schepper, Sandra, C,K., Heynderickx, Geraldine J. Marin, and Guy B. (2008) CFD Modeling of All Gas Liquid and Vapor Liquid Flow Regimes Predicted by the Baker Chart, Chemical Engineering Journal. Vol 138, 349357. in Two Phase Flow: A Comparison of Existing Correlations for Pressure Lo ss and Holdup, AIChE J. Vol 10, No. 1, 3843. Eckard, S.E, (1975) Multi Vane Expander as Prime Mover in Low Temperature Solar or Waste Heat Applications, Record of the Tenth IECEC Published by IEEE, Newark, Delaware, August 1822. Fairbrother, F., and Stubbs, A.E., (1935) The BubbleTube Method of Measurement, Journal of Chemical Society Vol 1, 527 529. Fluent 6.3.26 Users' Guide, 2006. Fluent Inc., Lebanon, NH, 2006. Fukano, T., and Kariyasaki, A., (1993) Characteristics of Gas Liquid Two Phase Flow in A Capillary, Nucl. Engng. Des. 141,5968. Fukagata, K., Kasagi, N., and Ua arayaporn, P., (2007) Numerical Simulation of Gas Liquid Two Phase Flow and Convective Heat Transfer in a Micro Tube, Int. J. Heat Fluid Flow Vol 28, 7282. Galbiati, L. and Andreini, P., (1994) Flow Pattern Transition for Horizontal Air Water Flow in Capillary Tubes. A Microgravity "Equivalent System'' Simulation, Int. Commun. Heat Mass Transfer, Vol. 21, 461468. Geiger, G.E., (1964) Sudden Contraction Losses in Single and Two Phase Flow, PhD thesis, University of Pittsburgh, Pittsburgh, PA. PAGE 158 158 Ghiaasiaan, S.M. and Abdel Khalik, S.I., (2000) TwoPhase Flow in Microchannels, Adv. Heat Transfer Vol 34, 145254. Graetz, L., (1883) Ueber die Wrmeleitungsfhigkeit von Flssigkeiten, Annalen der Physik und Chemie, Vol. 18 79. Granick, S., Zhu, Y., and Lee, H., (2003) Slippery Questions about Complex Fluids Flowing Past Solids, Nature Materials Vol 2, 221227. Gupta, R. Fletcher, D. F., and Haynes, B. S., (2009) On the CFD Modeling of Taylor Flow in Microchannels, J. Chemical Engineering Science, Vol. 64 29412950 Hayashi, S., Kasagi, N., and Suzuki, Y., (2007) The Effects of Inlet Flow Conditions On Gas Liquid Two Phase Flow in A Micro Tube, Proceeding of ASME JSME Thermal Engineering Summer Heat Transfer Conference, Vancouver, British Columbia, CANADA, July 812. He, Q., Fukagata, K., and Kasagi, N., (2007) Numerical Simulation of Gas Liquid Two Phase Flow and Heat Transfer with Dry Out in A Micro Tube, Proceedings of 6th International Conference on Multiphase Flow, Leipzig, Germany. He, Q. and Kasagi, N., (2008) PhaseField Simulation of Small Capillary Number TwoPhase Flow in a Microtube, Fluid Dynamics Research Vol 40, No. 78, July August, 497509. Hewitt, G.F., Shires, G.L., and Bott, T.R., (1993) Process heat transfer CRC Press, Ann Arbor, Michigan. Hirt, C.W., and Nichols, B.D., (1981) Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries, J. Computational physics. Vol. 39, 201225. Issa, R.I., (1986) Solution of Implicitly Discretized Fluid Flow Equations by Operator Splitting, J. Computational physic, Vol 62, 4065. Kandlikar, S.G., (2002) Fun damental Issues Related to Flow Boiling in Minichannels and Microchannels, Experimental Thermal and Fluid Science, Vol 26, 389407. Kashid, M.N., Platte, F., Agar, D.W., and Turek, S., (2007) Computational Modeling of Slug Flow in a Capillary Microrea ctor, J. of Computational and Applied Mathematics, Vol. 203(2), 487497. Kashid, M.N., Agar, D.W., and Turek, S., (2007) CFD Modeling of Mass Transfer with and without Chemical Reaction in the LiquidLiquid Slug Flow Microreactor, J. of Chemic al Engineering Science, Vol. 62, 51025109. PAGE 159 159 Kawahara, A., Chung, P.M. Y., Kawaji, M., (2002) Investigation of TwoPhase Flow Pattern, Void Fraction and Pressure Drop in a Microchannel, Int. J. Multiphase Flow, Vol 28, 14111435. Kays, W.M., (1950) Loss Coefficient for Abrupt Changes in Flow Cross Section with Low Reynolds Number Flow in Single and Multiple Tube System, Trans. ASME Vol 72, 1067 1074. Kreutzer, M.T., Kapteijn, F., Moulijn, J.A, Kleijn, C.R., and Heiszwolf, J.J., (2005) I nertial and Interfacial Effects on Pressure Drop of Taylor Flow in Capillaries, A.I.Ch.E. Journal, Vol 51 (9), 2428 2440. Kumar, V., Vashisth, S., Hoarau, Y., and Nigam, K.D.P., (2007) Slug Flow in Curved Microreactors: Hydrodynamic Study, Che mical Engineering Science, Vol .62 (24), 7494 7504. Kuru, W.C., Sangalli, M., Uphold, D.D., and McCready, M.J., (1995) Linear Stability of Stratified Channel Flow, International Journal of Multiphase Flow Vol 21, No. 5, 733753. Lahe y Jr., T.R. and Moody, F.J., (1993) The Thermal Hydraulic of Boiling Water Nuclear Reactors, 2nd ed., American Nuclear Society, La Grange Park, Illinois. Langhaar, L., and Lafayette, H., (1942) Steady in the Transition Length of a Straight Tube J. Applied Mechanics. Vol 64, A 55. Li, L., Zhao, Y., Guo, B., Shu, P., Shen, J., and He, S., (2003) Wrap of Cylinder and Its Effect on Main Features of Rotary Vane Compressor for Automobile Air Conditioning System, International Journal of Refriger ation Vol 26, 566574. Mahmoud, A. M., (2008) Analytical and Experimental Investigation of Rotary Vane TwoPhase Expanders in Vapor Compression Refrigeration Systems, PhD Thesis, University of Florida, Gainesville, FL. Marlene, E., (2000) Healing wi th Aromatherapy p. 9, McGraw Hill, ISBN 0658003828 Mehdizadeh, A., Sherif, S.A., and Lear, W.E., (2009a) "A Combined Analytical Numerical Model for a Two Phase Flow through a Sudden Area Change in Microchannels," 47th AIAA Aerospace Sciences Meeting and Exhibit Orlando World Center Marriott, Orlando, Florida, 58 January. Mehdizadeh, A., Sherif, S. A., and Lear, W.E., (2009b) CFD Modeling of TwoPhase Gas Liquid Slug Flow Using VOF Method in Microchannels, A Proceedings of Fluids Engineering Division Summer FEDSM2009, August 26, Vail, Colorado, USA, FEDSM2009 78438. PAGE 160 160 Mehdizadeh, Sherif, S. A., and Lear, W.E., (2009c) Numerical Simulation of Two Phase Slug Flows in Microchannels, A Proceedings of ASME Summer Heat Transfer Conference, HT2009, July 19 23, San Francisco, California, USA, HT200988126. Muzaferija, S.M., Sames, P., and Schell in, T., (1998) A Two Fluid Navier Stokes Solver to Simulate Water Entry, In Proc 22nd Symposium on Naval Hydrodynamics Washington, DC, 277289. Muzychka, Y.S., and Yovanovich, M.M., (2004) Forced Convection Heat Transfer in the Combined Entry Region o f Non Circular Ducts, ASME Journal of Heat Transfer Vol .126,5461. Muzychka, Y.S., and Yovanovich, M.M., (2002) Laminar Flow Friction and Heat Transfer in NonCircular Ducts and Channels Part II Thermal Problem, Compact Heat Exchangers A Festschrift on the 60th Birthday of Ramesh K. Shah, Grenoble, France August 24, 2002, 131139. Myers, R., (2003) The Basics of Chemistry Greenwood Publi shing Group p. 14, ISBN 0313316643 Schmidt, J. and Friedel, L. (1997) Two Phase Pressure Drop across Sudden Contraction in Duct Areas, Int. J. Multiphase Flow Vol 23, 283299. Serizawa, A., Feng, Z., and Kawara, Z., (2002) TwoPhase Flow in Microchannels, Experimental Thermal and Fluid Science Vol 26(67), 703714. Triplett, K.A., Ghiaasiaan, S.M., Abdel Khalik, S.I., and Sadowski, D.L., (1999) Gas Liqu id Two Phase Flow in Microchannels Phase Flow Pattern, Int. J. Multiphase Flow, Vol 25, 377394. Peterson, C.R. and McGahan, W.A., (1972) Thermodynamic and Aerodynamic Analysis Methods for Oil Flooded Sliding Vane Compressors, in: Proceedi ng of the Purdue Compressor Technology Conference, July 25 27, 18. Qian, D., and Lawal, A., (2006) Numerical Study on Gas and Liquid Slugs for Taylor Flow in A T Junction Microchannel, Chemical Engineering Science, Vol. 61 (23), 7609 7625. Pitman, V., (2004) Aromatherapy: A Practical Approach, Nelson Thornes p. xi, I SBN 0748773460 Robertson, G.F., and Wolgemuth, C.H., (1975) Analysis and Test Apparatus for a Vane Expander Using Steam, Record of the Tenth IECEC Published by IEEE, Newark, Delaware, August 1822. PAGE 161 161 Robertson, G.F and Wolgemuth, C.H., (1978) Exper imental and Analytical Study of Friction, Leakage, and Heat Transfer in a Vane Expander, Society of Automotive Engineers (SAE) Paper No.789521, 1430 1441. Serizawa, A., Feng, Z., and Kawara, Z., (2002) "Two Phase Flow in Microchannels," Experimental Thermal and Fluid Science Vol 26, No. 67, 703714. Taylor, G.I., (1960) Deposition of A Viscous Fluid on the Wall of A Tube, Journal of Fluid Mechanics Vol 10, 161 165. Tretheway, D.C., and Meinhart, C.D., (2002) Apparent Fluid Slip at Hydrophobic Microchannel Walls, Physics of Fluids Vol 14, No. 3, 912. Triplett, S.M.G.K.A., Abdel Khalik, S.I., LeMouel, A., and McCoud, B.N., (1999), Gas Liquid Two Drop, Int. J. Multiphase Flow Vol 121, B 263276. Walsh, A.P., and Walsh, J.E., (2009) Laminar Slug Flow Heat Transfer Characteristics Wi th Constant Heat Flux Boundary, A Proceedings of ASME Summer Heat Transfer Conference HT2009 July 1923, San Francisco, California, USA, HT200988428. Youngs, D. L., (1982) Time Dependent Multi Material Flow with Large Fluid Distortion, Numeri cal Methods for Fluid Dynamics K. W. Morton and M. J. Baines (ed.), Academic Press, London, 273285. Yu, Z., Hemminger, O., and Fan, L.S., (2007) Experiment and Lattice Boltzmann Simulation of TwoPhase Gas Liquid Flows in Microchannels, Chemical Engin eering Science. Vol. 62 (24), 7172 7183. Zhao, T.S. and Bi, Q.C., (2001) Co Current Air Water Two Phase Flow Pattern in Vertical Triangular Microchannels, Int. J. Multiphase Flow Vol 21, 765782. Zivi, S.M., (1964) Estimation of Steady State Steam V oidFraction By Means of the Principle of Minimum Entropy Production, J. Heat Transfer Vol 68, 247252. PAGE 162 162 BIOGRAPHICAL SKETCH Ayyoub M. Mehdizadeh was born in 1981, in Tehran, Iran. Ayyoubs excellence in academic background received his Bachelor of S cience from Guilan University in m echanical e ngineering. He was admitted to Sharif University of Technology where he received his Master of Science in 2006. During a presentation in an A SME conference Ayyoub meet Professor Sherif and was offered to peruse his Ph D work at University of Florida under his supervision. Ayyoubs multiple publications and team leader ship at the University of Floridas Industrial Assessment Center and practical experience in the fields of manufacturing and energy management has enabled him to propose > 5.1 million dollars saving s in cost for 21 audited big industri al facilities Ayyoub has also been active in a plethora of professional activities while in graduate school. He is a member of American society of Mechanical Engineer ing (ASME) and a member of American Institute of Aeronautics and Astronautics (AIAA). Ayyoub will start his career as a postdoc toral fellow and lab manager in the Two Phase F low L aboratory under the supervision of Professor Klausner at the M echanical and A erospace Engineering D epartment of the University of Florida. 