UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 1 FORCED CONVECTION FLOW BOILING A ND TWOPHASE FLOW PHENOMENA IN A MICROCHANNEL By YUN WHAN NA A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2008 PAGE 2 2 2008 Yun Whan Na PAGE 3 3 To my parents and my wife PAGE 4 4 ACKNOWLEDGMENTS This dissertation could not been accomplished if not for the continuous advice, assistance, and support of m y advisor, professors, my parents and my wife. First, I would like to thank to my advisor Dr. Jacob N. Chung who supports, initiate, and guide the res earch during my Ph.D. study. His support, advice, and encouragement have made it possible for me to finish the study. Many thanks are expressed to Dr. James Klausner, Dr. William Lear, Jr., Dr. S. A. Sherif, and Dr. ChangWon Park for serving as member s of dissertation committee. Their precious advice and recommendations has considerably help improve this study. A special thanks goes to my colleagues in my research group who have given their valuable time and experience. I would like to express my appr eciation to my parents, brot her and his wife, and nephews for their patience and encouragement. Most of all I would like to express a great debt of gratitude to my wife, Ji Hye Choi, for her love, patience, and support with me over the past few years. I dedicate this dissert ation to my wife. PAGE 5 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS ............................................................................................................... 4LIST OF TABLES ...........................................................................................................................8LIST OF FIGURES .........................................................................................................................9NOMENCLATURE .................................................................................................................. ....13ABSTRACT ...................................................................................................................... .............18CHAPTER 1 INTRODUCTION .................................................................................................................. 201.1 Research Background .......................................................................................................221.2 Research Objectives ..........................................................................................................231.3 Organization of the Dissertation ....................................................................................... 242 BACKGROUND AND LITER AT URE REVIEW ................................................................ 282.1 TwoPhase Flow Pattern ................................................................................................... 282.1.1 Flow Patterns in Horizontal Tubes ......................................................................... 282.1.2 Flow Patterns in Minicha nnels and Microchannels ................................................ 302.1.3 Flow Regime Map ..................................................................................................342.2 Pressure Drop during Flow Boiling in Minichannels a nd Microchannels ....................... 392.2.1 SinglePhase Pressure Drop ...................................................................................402.2.2 TwoPhase Pressure Drop ......................................................................................422.3 TwoPhase Flow Boiling in Microchannels ..................................................................... 462.3.1 PhaseChange Phenomena in Meniscus ................................................................. 472.3.2 Simple Model for TwoPhase Flow Boiling .......................................................... 502.3.3 Numerical Simulation of TwoPhase Boiling ........................................................ 523 SEMIANALYTICAL EVAPORATION MODEL ............................................................... 623.1 Assumptions .....................................................................................................................623.2 Mathematical Modeling .................................................................................................... 633.2.1 Mass and Energy Balance ......................................................................................633.2.2 Liquid Pressure Gradient ........................................................................................643.2.3 Vapor Pressure Gradient .........................................................................................653.2.4 Interfacial Shear Stress ...........................................................................................663.2.5 Liquid Temperature Variation ................................................................................673.2.6 Pressure Difference at th e LiquidVapor Interface.................................................683.2.7 Liquid Film Thickness ............................................................................................683.2.8 Interfacial Temperature .......................................................................................... 69 PAGE 6 6 3.2.9 Mean Velocities of Li quid and Vapor Phases ........................................................ 703.3 Solution Procedures ..........................................................................................................713.4 Results and Discussion .................................................................................................... .743.4.1 Constant Heat Flux .................................................................................................753.4.2 Constant Wall Temperature ....................................................................................793.4.3 Pressure Drop .........................................................................................................824 NUMERICAL SIMULATION OF FORCED CONVECTIVE BOILING IN A MICROCHANNEL ..............................................................................................................1014.1 Background and Literature Review ................................................................................1014.1.1 Onset of Nucleate Boiling (ONB) ........................................................................ 1034.1.2 Bubble Departure Diameter .................................................................................. 1054.1.3 Bubble Departure Frequency ................................................................................ 1064.2 Forced Convection Flow Boiling Numerical Simulation Method .................................. 1084.2.1 Governing Equations ............................................................................................1084.2.2 Numerical Method ................................................................................................ 1094.2.2.1 Finite volume method (FVM) .................................................................... 1094.2.2.2 Volume of Fluid (VOF) .............................................................................. 1134.3 Problem Formation and Description ............................................................................... 1184.3.1 TwoDimensional Model ...................................................................................... 1184.3.2 Numerical Infrastructure ...................................................................................... 1184.4 Results and Discussion ................................................................................................... 1224.4.1 Onset of Nucleate Boiling .................................................................................... 1224.4.2 Effects of Channel Height ....................................................................................1244.4.3 Bubble Coalescence and Departure ......................................................................1254.4.4 TwoPhase Flow Patterns .....................................................................................1264.4.5 Heat Transfer and Flow Asso ciated with Elongated Bubbles ..............................1274.4.6 Heat Transfer ........................................................................................................1275 TWOPHASE FLOW INSTABILITIES IN A TRAPEZOI DAL MICROCHANNEL ....... 1475.1 TwoPhase Flow Instabilities .........................................................................................1485.1.1 Static and Dynamic Flow Instabilities .................................................................. 1485.1.1.1 Static instabilities .......................................................................................1495.1.1.2 Dynamic instabilities ..................................................................................1505.1.2 Microchannel Forced Convecti on Flow Boiling Instabilities .............................. 1515.1.3 Research Objectives .............................................................................................1565.2 Problem Descriptions ......................................................................................................1575.3 Computational Model .....................................................................................................1585.4 Results and Discussion ................................................................................................... 1605.4.1 Flow Boiling Patterns ...........................................................................................1615.4.2 Wall Temperature Fluctuations ............................................................................ 1625.4.3 Pressure Fluctuations ............................................................................................1645.4.4 Stability Map ........................................................................................................ 1666 CONCLUSIONS AND RECOMMENDATIONS ............................................................... 185 PAGE 7 7 6.1 Summary and Conclusions .............................................................................................1856.2 Recommendations for Future Research ..........................................................................187APPENDIX A DERIVATION OF FILM THICKNESS AND E QUILIBRIUM NONEVAPORATING FILM THICKNESS ..............................................................................................................188A.1 Liquid Film Thickness ................................................................................................... 188A.2 Equilibrium NonEvaporating Liquid Film ................................................................... 192B DERIVATION OF LIQUID AND INTERFACI AL TEMPERATURE ..............................194B.1 Liquid Temperature Profile ............................................................................................ 194B.2 Interfacial Temperature .................................................................................................. 197C NUSSELT NUMBER AND HEAT TRANSFE R COE FFICIENT CORRELATION ........200C.1 Nusselt Number ............................................................................................................ ..200C.2 Heat Transfer Coefficient Correlation ...........................................................................200LIST OF REFERENCES .............................................................................................................201BIOGRAPHICAL SKETCH .......................................................................................................215 PAGE 8 8 LIST OF TABLES Table page 21 Classification of channel dimensions .................................................................................5922 Empirical parameter C in LockhartMartinelli correlation ............................................... 5923 Summary of twophase fricti onal pressure drop correlations ............................................ 6041 Dimensions of the microchannel used in computational simulation ............................... 14542 Thermophysical properties of working fluid and heat sinks used in computational simulation .................................................................................................................... .....14543 Grid density for each case used in computational simulation .......................................... 14544 Effects of channel he ight on nucleate boiling .................................................................. 14651 Summary of investigation on flow instabil ities in minichannels and microchannels ...... 18052 Dimensions of the microchannel used in computational simulation ............................... 18453 Grid density used in computational simulation ...............................................................18454 Summary of the amplitude and oscillati on of inlet pressure under different heat fluxes and mass fluxes ..................................................................................................... 18455 Summary of stable and un stable flow boiling regime under different heat fluxes and mass fluxes ................................................................................................................... ....184 PAGE 9 9 LIST OF FIGURES Figure page 11 Microchannel heat sinks. ...................................................................................................2612 Heat transfer coefficient variations for natural convection, singlephase and twophase forced convection for different coolants. .................................................................2613 Vapor reverse flow occurred in parallel microchannels. ................................................... 2721 Schematic representations of flow regi mes observed in horizontal, cocurrent gasliquid flow. .................................................................................................................. .......5522 Flow patterns in small para llel rectangular channels. ........................................................ 5523 Twophase flow patterns and transiti ons of flow boiling in microchannels. ..................... 5624 New diabatic coalescing bubble map for eva porating flow in circular microchannels. ....5725 Flow regions through a liquid thin film in a microchannel. ..............................................5726 Interline junction of vapor, adsorb ed evaporating liquid thin film and nonevaporating adsorbed liquid film. ...................................................................................... 5831 Physical model for evaporation phenomena through the liquidvapor interface in a microchannel. ................................................................................................................. ....8532 Geometry and coordinate system used in a numerical simulation. ....................................8533 Microchannel heat sink unit cell used in a numerical simulation. ..................................... 8534 Solution procedure for ob taining unknown variables. .......................................................8635 Effects of heat flux on the liquid film thickness. ...............................................................8736 Evaporative mass flux variations at different heat fluxes. .................................................8737 Liquid pressure grad ient variations at di fferent heat fluxes. .............................................. 8838 Effects of heat flux on pressure difference along the axial direction. ............................... 8839 Mean velocity variations of the li quid phase at different heat fluxes. ...............................89310 Mean velocity variations of the va por phase at different heat fluxes. ...............................90311 Interfacial shear stress variati ons at different heat fluxes. ................................................. 90312 Effects of heat flux on interf acial temperature variations. .................................................91 PAGE 10 10 313 Comparison of channel wall and heat sink bottom tem perature to experimental results investigated by Qu and Mudawar (2004). ..............................................................91314 Local heat transfer coefficient va riations for different heat fluxes. ................................... 92315 Average heat transfer coefficient va riations with various heat fluxes. .............................. 92316 Heat transfer coefficient variations with different thermodynamic equilibrium quality. ...................................................................................................................... .........93317 Liquid film thickness variations for different wall temperatures. ...................................... 93318 Effects of wall temperat ure on evaporative mass flux. ...................................................... 94319 Liquid pressure gradie nt variations for differe nt wall temperatures. ................................. 94320 Effects of wall temperature on pressure difference along the axial direction. ................... 95321 Mean velocity profiles of the liquid phase at different wall temperatures. ....................... 96322 Mean velocity profiles of the vapor phase at different wall temperatures. ........................ 96323 Influences of wall temperatur e on interfacial shear stress. ................................................97324 Interfacial temperature profiles for different wall temperatures. ....................................... 97325 Effects of wall temperature on lo cal heat transfer coefficient. .......................................... 98326 Average heat transfer coefficient variations for different wall temperatures. ................... 98327 Microchannel heat sinks operating at subcooled condition. .............................................. 99328 Comparison of pressure drops of the cu rrent model with the experimental results observations by Qu and Mudawar (2004) for different inlet temperatures. .................... 10041 Ebullition mechanism of nucleation and bubble growth a nd departure at an active cavity site. .................................................................................................................. ......13042 Bubble release frequency with the waiting time and the growth time. ............................ 13043 Control volume used to disc retize the governing equations. ........................................... 13044 Twophase cell with an im mersed phase interface. .........................................................13145 Measurement of the contact angle near a heated wall. .................................................... 13146 Liquidvapor interface calculation us ing a volume fraction in each cell. ........................ 13147 Calculation of volume flux th rough the computational face. .......................................... 132 PAGE 11 11 48 Physical model of a microchannel us ed in a com putational simulation. ......................... 13249 Grid system and boundary conditions used to simulate twodimensional twophase flow boiling model. .......................................................................................................... 133410 Computational procedure in cluding mass transfer and phase change calculations. ........ 133411 Wall temperature profiles at incipient boiling for different channel heights. .................. 134412 Effects of mass flow rate on wall temper ature profile during in cipient boiling at Hch = 500 m. .........................................................................................................................135413 Effects of heat flux on wall temperat ure profile during incipient boiling for Hch = 300 m. ...................................................................................................................................135414 Wall temperature fluctuations while bubbles generate everywhere in a microchannel for Hch = 30 m. ...............................................................................................................136415 Onset of nucleate boiling on a heated su rface at different mass flow rates for Hch = 1mm. .......................................................................................................................... ......136416 Bubble departure procedure in a microchannel with Hch = 700 m. ...............................137417 Bubble departure procedure in a microchannel with Hch = 500 m. ...............................138418 Bubble departure procedure in a microchannel with Hch = 300 m. ...............................139419 Bubble departure diameter variations for different channel heights at same heat flux and mass flux. ................................................................................................................ ..139420 Contour plots of vapor volume frac tion at various simulation times for Hch = 1 mm. .... 140421 Contour plots of vapor volume frac tion at various simulation times for Hch = 30 m. ...140422 Wall temperature and velocity profile underneath a bubble with local dryout near the downstream region for Hch = 30 m. .........................................................................141423 Wall temperature and velocity profile underneath a slug bubble near the downstream region for Hch = 30 m. ....................................................................................................142424 Contour and vector plots of va por bubbles in a microchannel with Hch = 300 m. ........143425 Heat transfer coefficient variations w ith quality for different channel heights. .............. 143426 Heat transfer coefficient variations with quality for different heat fluxes. ...................... 14451 Ledinegg instability based on pressu re difference and mass flow rate. ........................... 168 PAGE 12 12 52 Physical diagram of a trapezoidal microchannel. ............................................................ 16853 Symmetric flow channel geometry us ed in a computational simulation. ........................ 16954 Grid system of a threedimensional CFD model. ............................................................ 17055 Flow patterns in the middle of a tr apezoidal microchannel for Case I. ...........................17056 Flow patterns in the middle of a tr apezoidal microchannel for Case II. ..........................17157 Wall temperature fluctuations with small amplitude and shortperiod oscillation and flow patterns for Case I. ...................................................................................................17258 Wall temperature fluctuations with large amplitude and shortperiod oscillation and flow patterns for Case II. .................................................................................................17459 Wall temperature fluctuations with large amplitude and shortperiod oscillation for Case III. ..................................................................................................................... .......175510 Inlet and outlet pressure fluctuations with small amplitude and shortperiod oscillation for Case I. .......................................................................................................176511 Inlet pressure fluctuations with large amplitude and shortperiod oscillation for Case II. ........................................................................................................................... ...........177512 Inlet pressure fluctuations with large amplitude and shortperiod oscillation for Case III............................................................................................................................ ..........178513 Stability map in subcooling number versus phase change number based on Chang and Pan (2007). ................................................................................................................179 PAGE 13 13 NOMENCLATURE A Hamaker constant [J] A Area [m2]; Constant coefficient a, b Constant coefficient Bo Boiling number C Empirical parameter; Constant coefficient cp Specific heat [J/kgK] Ca Capillary number Co Confinement number; Convection number D Diameter of the flow cha nnel [m]; Bubble diameter [ m] Eo Etovos number f Fanning factor or frictio n factor; Frequency [Hz] FTD Dimensionless parameter for wavyannular and wavyintermittent transitions Fr Froude number G Mass flux [kg/m2s] g Gravitational acceleration constant, [m/s2] H Channel height [mm] hlv Latent heat of vaporization [J/kg] hc Heat transfer coefficient [W/m2K] j Superficial velocity [m/s] Ja Jakob number Ja* Modified Jakob number k Heat conductivity [W/mK] K Curvature of the meniscus [m1] PAGE 14 14 KTD Dimensionless parameter for stratified flowtowavy flow transition L Length [mm] M Liquid molecular weight [kg/mol] M Mass flow rate per unit length [kg/ms] m Mass flow rate [kg/s] m Mass flux [kg/m2s] N Total number; Dimensionless parameter Npch Phase change number Nsub Subcooling number Nu Nusselt number Pw Heated perimeter [m] p Pressure [Pa] Q Heat transfer rate [W] q Heat flux [W/m2] R Channel radius [m] r Radial coordinate [m]; Cavity size [m] g R Gas constant [J/kgK] u R Universal gas constant [J/molK] Re Reynolds number S Slip ratio T Temperature [K] T Bulk temperature [K] TTD Dimensionless parameter for bubbly fl owtointermittent flow transition t Time [s] PAGE 15 15 u Velocity [m/s] u Mean velocity [m/s] V Molar volume [m3/mol] v Blowing velocity [m/s] W Channel width [mm] We Weber number x Coordinate [m]; Vapor quality X Martinelli parameter Y Coordinate [m] Z Coordinate [m] Greek letters Accommodation coefficient; Volume fraction Momentum flux coefficient Variable difference; Firstorder derivative of Film thickness [m] 0 Nonevaporating film thickness [m] Firstorder derivative of Contact angle []; Inclination angle of a channel [] Laplace constant; Mean free path [m] Dynamic viscosity [kg/ms2] Kinematic viscosity [m2/s] Density [kg/m3] Surface tension [N/m] Shear stress [N/m2] Twophase multiplier PAGE 16 16 Subscripts 2 Twophase mixture a Acceleration b Boiling position bot Bottom of microchannel heat sinks c Capillary pressure; Critical radius; Crosssectional area ch Flow channel conf Confinement CBD Convective boiling dominant d Disjoining pressure; Bubble departure e Thermodynamic equilibrium eff Effective heat flux imposed on microchannel heat sinks evap Evaporative fr Frictional g Growth time; Gas phase i Liquidvapor interface in Inlet h Hydraulic HS Heat sinks l Liquid lo Liquid only NBD Nucleate boiling dominant p Phase p; Plenum pch Phase change ONB Onset of nucleate boiling PAGE 17 17 out Outlet q Phase q s Shear force; Heated surface sat Saturation sub Subcooling t Total th Threshold top Top v Vapor vo Vapor only w Wall of a heated flow channel; Width; Waiting time z Axial direction PAGE 18 18 Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy FORCED CONVECTION FLOW BOILING A ND TWOPHASE FLOW PHENOMENA IN A MICROCHANNEL By Yun Whan Na December 2008 Chair: Jacob N. Chung Major: Mechanical Engineering The present study was performed to numeri cally analyze the evaporation phenomena through the liquidvapor interface and to investigat e bubble dynamics and heat transfer behavior during forced convective flow boiling in a microchannel. Flow instabilities of twophase flow boiling in a microchannel were studied as well. The main objective of this research is to investigate the fundamental mechanisms of twophase flow boiling in a microchannel and provide predictive tools to design thermal manage ment systems, for example, microchannel heat sinks. The numerical results obtained from this study were qualitativ ely and quantitatively compared with experimental re sults in the open literature. Physical and mathematical models, accoun ting for evaporating phenomena through the liquidvapor interface in a microc hannel at constant heat flux and constant wall temperature, have been developed, respectively. The heat transf er mechanism is affected by the dominant heat conduction through the thin liquid film and vaporization at the liquidvapor interface. The thickness of the liquid film and the pressure of the liquid and vapor phases were simultaneously solved by the governing differential equations. Th e developed semianalytical evaporation model that takes into account of the interfacial phenomena and surface tension effects was used to obtain solutions numerically usi ng the fourthorder RungeKutta method. The effects of heat flux PAGE 19 19 and wall temperature on the liquid film were evaluated. The obtained pressure drops in a microchannel were qualitatively consistent with the experimental results of Qu and Mudawar (2004). Forced convective flow boiling in a single microchannel with different channel heights was studied through a numerical simulation to in vestigate bubble dynamics, flow patterns, and heat transfer. The momentum and energy equations were solved using the finite volume method while the liquidvapor interface of a bubble is captured using the VOF (Volume of Fluid) technique. The effects of different constant heat fluxes and different channel heights on the boiling mechanisms were investigated. The eff ects of liquid velocity on the bubble departure diameter were analyzed. The obtai ned results showed that the wa ll superheats at the position of nucleate boiling are relatively inde pendent of the mass flow rates at the same channel height. The obtained results, however, showed that the heat flux at the on set of nucleate boiling strongly depends on the channel height. With a decrease of the channel height and an increase of the liquid velocity at the channel in let, the departure diameter of a bubble was smaller. The periodic flow patterns, such as the bubbly flow, elongated slug flow, and churn flow were observed in the microchannel. Flow instabilities of twophase fl ow boiling in a trapezoidal microchannel using a threedimensional model were investigated. Fluc tuation behaviors of flow boiling parameters such as wall temperature and inlet pressure ca used by periodic flow patterns were studied at different heat fluxes and mass fl uxes. The numerical results showed large amplitude and short period oscillations for wall temperature and inlet pr essure fluctuations. Stable and unstable flow boiling regime with short period oscillations were investigated. Those flow boiling regimes were not listed in stable and unstable boiling regime map proposed by Wang et al. (2007). PAGE 20 20 CHAPTER 1 INTRODUCTION As the techn ologies of electronic systems have developed in the last two decades, the size of the devices becomes smaller and some of them require dissipating high heat fluxes (up to 300 W/cm2). The recently hottest issues of thermal manage ment in electronic systems are to maintain high device performance and reli ability. The cooling device for highpowerdensity electronic systems needs the very large heat dissipation fr om small areas, and mini/microdimensions in recent years. It is necessary to develop new methods for distributing and removing heat from highpowerdensity of electronic devices. The surface area to volume ratio increases as the diameter of a channel from macr oscale to microscale (up to 10 m) decreases. Thus, with very high heat transfer coefficients microchannels having a hydraulic diameter less than the Laplace constant, lvg are well suited for electronic cooling systems requiring for removing higher heat fluxes. Figure 11 shows a schematic of microchannel heat sinks for dissipating a high heat flux. Microchannel heat sinks comp rise numerous parallel microchannels. The microchannel heat sinks dissipate the heat generated from a heat source by heat conduction and remove it by forced convection of bulk liqui ds flowing through the microchannels. Microchannel cooling systems using singlephase flow have been extensively studied in the past years (Tuckerman and Pease, 1981; Philips, 1990; Peng et al., 1994; Fedorov and Viskanta, 2000; Wu and Cheng, 2003a and 2003b; Li et al., 2004; Lee and Garimella, 2006). Tuckerman and Pease first proposed electronic c ooling systems using microchannel heat sinks having a width of 50 m and a depth of 302 m and separated by fins with a width of 50 m. Fluid flow and heat transfer model for microcha nnel heat sinks using si nglephase flow were investigated by Philips (1990). He poi nted out that microchannel heat sinks increase significantly thermal performance. It was found that thermal resistance as low as 0.1 C/(W/cm2) and PAGE 21 21 pumping power requirement less than 10 W/cm2 can be obtained when using microchannel heat sinks as well. The detailed st udies on fluid flow and thermal transport phenomena of singlephase microchannels can be found in review pa pers (Sobhan and Garimella, 2001; Garimella and Sobhan, 2003; Agostin i et al., 2007). Due to higher heat removal capabilities, on th e other hand, the theoretical and experimental studies of twophase flow boiling in microcha nnels have been investigated by a number of researchers in late years as a prime application for thermal management of next generation highpowerdensity commercial and defense elec tronic devices. Garimella and Sobhan (2003) reviewed the recent results of theo retical and experimental studies on fluid flow and heat transfer in microchannels and microtubes. Thome (2006) provided the comprehensive review of boiling and twophase flow including boiling phenomena and flow instabilities in microchannels. However, transport phenomena including flow instabilities of twophase flow boiling in microchannels have not been su fficiently understood due to differe nt characteristics comparable to flow boiling in macroscales. Therefore, prope r understanding and predic tive tools of fluid flow and heat transfer in microchannels are essential for their desi gn and safe operation. This research works focused on three topics for twophase flow boiling in a microchannel. First, a semianalytical model of elongated slug bubble flow has b een developed to provide the fundamental understanding of transport phenomena of twophase flow boiling in a microchannel. Second, fluid flow near a vapor bubble, onset of nucleate boiling, and flow boiling in the microchannel have been investig ated to analyze b ubble dynamics under various channel heights and flow boiling conditions. Finally, twophase fl ow instabilities occurred during a flow boiling process in a microchannel have been numerically studied to an alyze stable and unstable flow boiling regimes at different heat fluxes and mass fluxes. PAGE 22 22 1.1 Research Background W ith the rapid development of microelectro mechanical system (MEMS) technology in the past decades, microchannels are widely used in miniature heat exchangers, nuclear reactors, cooling devices in highpower electronic syst ems, microchemical reactors, biotechnology systems, and heat rejection systems of a spacecr aft. The cooling method using singleor twophase flows is extensively used for the therma l management of highpowerdensity electronic systems to dissipate higher heat fluxes in recent years. The cooling syst em using singlephase flow has significantly progressed in recent years (Zhang et al., 2005; Gernert et al., 2005). It utilizes the sensible heat exchange to rem ove heat leading to in creasing the bulk liquid temperature in microchannels. The singlephase cooling method has several disadvantages as follows: To reduce the increased bulk liquid temper ature in and surface temperature of microchannel heat sinks requires more mass flow rates. Higher pumping power may be needed due to the increased mass flow rate. Surface temperature drastically increases al ong the streamwise direction due to heat dissipation from a heated surface. The high surface temperature seriously affects device materials. It causes the performance and reliab ility of electronic devi ces to deteriorate. Pressure drops may largely increase due to the frictional force caused by the liquid velocity increased to maintain appropriate cooling effects. The cooling system using twophase flow boilin g dissipates higher heat by the latent heat of vaporization rather than the sensible heat exchange. It has several advantages as follows: Mass flow rate can be reduced because liquids in microchannels are mostly consumed by the latent heat of vaporization rather than sensible heat. It causes pumping power to decrease. Due to phase change from liquid to vapor while boiling initiates, uniform surface temperature can be obtained. It may provide much better performance and reliability of electronic devices. PAGE 23 23 Heat transfer coefficients of twophase flow boiling are much higher than those of singlephase flow as shown in Figure 12. Heat is efficiently removed from local hot spots due to the higher heat transfer coefficients as well. Although the cooling technique usin g the latent heat of vapori zation rather th an using the sensible heat is efficiently applicable to micr oelectronic devices requiri ng higher heat removal, the cooling method using twophase flow boiling ha s disadvantages such as large pressure drops, twophase flow instabilities causing severe pr essure drop and wall temperature fluctuations causing mechanical damage, and reverse flow ca used by rapid vapor bubble generation, growth and expansion in microchannels as shown in Figure 13. The heat transfer mechanism of twophase fl ow boiling in microchannels, such as flow regimes, flow patterns, pressure drops, flow boili ng heat transfer, and flow instabilities has been extensively studied by a number of researchers in the recent year s (Qu et al., 2004; Revellin et al., 2006; Wang et al., 2007; Huh et al., 2007). Transport phenomena of flow boiling have not been fully understood due to complex behavior of vapor bubbles when boiling in microchannels initiates. No optimum methodology has been still developed to explain the transport phenomena of twophase flow boiling in microchannels Thus, alternative predictive methods for understanding the physical phenomena of twophase flow boiling in microchannel are required to provide a better design t ool and safe operation. 1.2 Research Objectives The physical phenom ena of twophase flow boiling in microchannels should be entirely understood to be applied to microscale devices. Fu rthermore, the predictive tools for flow boiling in a microchannel are quite limited. To contribute the prediction of twophase flow boiling in a microchannel, two main mechanisms for twopha se flow boiling are numerically investigated: heat transfer through the liquidvapor interface and bubble dynami cs such as nucleate boiling, PAGE 24 24 bubble growth and departure. Flow instabilities leading to mechanical failure and an operating problem are numerically studied under certain operating conditions as well. Through this research works, the following objectives are achieved: Understanding the boiling transport phenomena of fluid flow and heat transfer in a microchannel. Developing the predictive tool of pressure drop and heat transfer for twophase flow boiling. Investigating bubble dynamics such as bubble ge neration, departure, and onset of nucleate boiling on a heated surface to provide better understanding of twophase flow boiling in a microchannel. Identifying stable and unst able boiling regimes in a trapezoidal microchannel under various heat fluxes and mass fluxes. 1.3 Organization of the Dissertation Current study focuses on the num erical simulatio n of twophase flow boiling to investigate heat transfer phenomena at the liquidvapor inte rface, bubble dynamics, and flow instabilities in a microchannel. This disserta tion is organized as follows: Chapter 1 provides an introduction and research objectives to the current work. Chapter 2 introduces the background of twopha se flow boiling and heat transfer in minichannels and microchannels. Previous experimental and numerical results for twophase flow are briefly reviewed and qualitatively assessed. Chapter 3 explains the mathematical model of heat transfer mechanism through the liquidvapor interface in detail. The governing equati on, boundary conditions, and solution procedure are described. The numerical results of constant wall temperature and constant heat flux are given as well. The pressure drops predicted by numerical results and experimental results are qualitatively and quantitatively evaluated. PAGE 25 25 Chapter 4 describes the computational simu lation method of twophase flow boiling to investigate bubble dynamics and heat transfer in a microchannel. Twodimensional numerical simulations are carried out using a finite vo lume technique for solving the momentum and energy equations while the liquidvapor moving interface is tracked by the volume of fluid (VOF) method. The effect of channel heights on vapor bubble and periodic flow boiling are discussed. Chapter 5 presents flow instabilities of twophase flow boiling in a trapezoidal microchannel under several operating conditions. Threedimensional computational simulations are carried out using transient tw ophase flow boiling model incor porated with mass transfer and phase change. Wall temperature and inlet pressure oscillations are desc ribed. To validate the numerical results, the results obt ained from the simulation are compared with the experimental data in open literature. Stab le and unstable boiling regi mes are identified as well. Chapter 6 summarizes the results obtained from the current studies and provides some recommendations for future studies. Appendix A describes derivations of film th ickness and equilibrium nonevaporating film thickness. Appendix B gives more detailed li quid temperature and in terface temperature. Appendix C describes Nusselt number for a threeside heated rectangular channel and heat transfer coefficient correlation of flow boiling in microchannels. PAGE 26 26 Figure 11. Microcha nnel heat sinks. Figure 12. Heat transfer coe fficient variations for natural convection, singlephase and twophase forced convection for different coolants (Mudawar, 2000). PAGE 27 27 Incoming subcooled liuiqd Liquid Vapor Microchannel wall Inlet plenum Figure 13. Vapor reverse flow o ccurred in parallel microchannels. PAGE 28 28 CHAPTER 2 BACKGROUND AND LITERATURE REVIEW Fluid flow and heat transport phenom ena in mi crochannels are significa ntly different from those in the macrochannels. First, flow boiling with energy transfer a nd momentum exchange through the liquidvapor interface is important in the microscale. Second, the surface tension force plays a significant role in microchannels due to the relative mu ch smaller gravitational force and thus the gravitational force in microcha nnels may be neglected. This effect dominantly affects the liquid thin film and controls the heat transfer mechanisms. The flow pattern such as bubblyslug flow observed in the microchannels is mainly affected by th e surface tension. Flow regimes and flow patterns are of important char acteristics of twophase flow boiling. They are essential to the development of reliable predictive tools for transport phenomena such as pressure drop, heat transfer, and flow st ability in microchannel heat si nks using the phasechange of liquids. Although numerous researches for twophase fl ow boiling in microchannels have been studied during the past several decades, the fundamental understanding of fluid flow and heat transfer characteristics in the microscale has been very limited. In this chapter, the previous works including twophase flow boiling, twoph ase flow pattern and flow regime maps, interfacial phenomena, pressure drop model, a nd numerical simulation for flow boiling in macroand microscale are reviewed. The current study in microchannel flow boiling is highlighted as well. 2.1 TwoPhase Flow Pattern 2.1.1 Flow Patterns in Horizontal Tubes The understanding of twophase flow patterns in horizontal tubes prov ides the inform ation of pressure drop and heat transfer phenomena. Dryout can be occurred because most liquids may PAGE 29 29 change into vapor bubbles while boiling takes place in macrocha nnels. It leads to drastically increasing surface temperature and then causes heat transfer to significantl y decrease. Therefore, the condition at which dryout or cr itical heat flux (CHF) occurs can be investigated by studying the flow patterns in flow boiling systems. Baker (1954) investigated early flow pattern maps for oilgas flow in macrochannels. The generalized flow pattern maps for gasliquid flow in horizontal pipes were investigated by Mandhane et al. (1974). They proposed various flow patterns: dispersed flow, bubble and elongated bubble flow, stratified flow, wave flow, slug flow, and annular and annular mist flow. Taitel and Dukler (1976) presented theoretically fl ow pattern maps with transition criteria for gasliquid flow in horizontal tubes and these flow pattern maps were widely used for horizontal flows. Empirical based flow pattern maps were investigated by Weis man et al. (1979) and Spedding and Spence (1993). Flow pa tterns in inclined tubes were observed by Taitel and Dukler (1976), Barnea et al. (1980), Weisman and Ka ng (1981), and Mukherjee and Brill (1985). The adiabatic flow patterns in circular tubes were comprehensively investigated by Hewitt (2000). The typical flow patterns for twophase flow in a horizontal round tube are illustrated in Figure 21. The characteristics of twophase flow pattern are de fined as follows (Carey, 1992): Bubbly flow: The vapor bubbles te nd to flow the upper portion of the horizontal tube due to their buoyancy effect at the very low quality. Plug flow: As vapor quality increases, the small bubbles merge each other and then produce large bubbles. The bubbles fl ow in the top of the tube. Stratified flow: The liquid flowing in the bottom of the tube is separated from the vapor in the top of the horizontal tube by a smooth surface at the low flow rates and higher qualities. Wavy flow: At the increased quality and fl ow rate, the interface becomes wavy and it becomes the wavy flow regime. Slug flow: As the vapor velocities increase, the flow regime transit from the wavy flow to the slug flow. The vapor slugs tend to move in the top of the tube due to the buoyancy. PAGE 30 30 Annular flow: Annular flow is observed at high vapor velocities. The film thickness on the bottom of the tube becomes thick, while the film on the top of the tube tends to be thin due to the buoyancy effect. 2.1.2 Flow Patterns in Minichannels and Microchannels Twophase flow characteristics of m inichanne ls and microchannels have been known to be significantly different from thos e in macrochannels. The gravitati onal force is more dominant than the surface tension force in macrochannels, while the surface tension force plays important role in minichannels and microchannels. The f undamental study of twophase flow patterns in minichannels and microchannels are still in vestigated due to the complex morphological configurations with different geom etries and operating conditions. Several investigators proposed the microc hannel criterion given by Laplace constant L or the Etovos number Eo L lvg (21) where g, l and v are the surface tension, gravitati onal acceleration, densities of the liquid and vapor at the satura ted pressure, respectively. 2 lvgD Eo (22) with a channel diameter, D. Suo and Griffith (1964) provide d the following criterion for a channel diameter, D, Confined number3.3LD (23) Brauner and MoalemMaron (1992) f ound the criterion as follows: 22 Eo (24) PAGE 31 31 Kandlikar and Grande (2003) presented the cl assification of channel dimensions based on the Knudsen number, mhKnD where the mean free path 2m R T as listed in Table 21. Flow patterns in isothermal gasliquid for small diameter tubes were experimentally studied by Fukano et al. (1989). They observed four flow patte rns: bubbly, plug, slug, and annular flow in the horizontal capillary tubes having a diameter of 1.0 mm, 2.4 mm, and 4.9 mm. Wambsganss et al. (1991) investigated compre hensively adiabatic twophase flow patterns and transition of airwater mixtures in a small horizontal rectangular channel with 19.05 mm 3.18 mm. Their results showed that the flow pattern maps in a macrochannel are not applicable to the small diameter channel. There is no stratified flow in those channels because of the surface tension is dominant in comparison with th e gravitational force. Cornwell and Kew (1992) observed experimentally three fl ow patterns such as isolated bubble, confined bubble, and annularslug flow in small parall el channels with R113 as show n in Figure 22. They showed that the flow patterns are strongl y related with heat transfer coefficient. Mertz et al. (1996) studied the flow boiling heat tr ansfer with water and R141b in rectangular channels with channel widths of 1 mm, 2 mm, a nd 3 mm and aspect ratios of up to 3. They investigated the presence of nucleate boiling, confined bubble flow for medium hat fluxes and annular flow for high heat fluxes. Their results s howed that the vapor generated in the channels blocks the liquid flow, causing a reverse flow to occur in the channels. Kuznetsov and Shamirzaev (1999) investigated flow patterns during flow boiling of R318C in an annul ar channel with 0.9 mm gap. They observed four flow pattern s: small bubbles, confined bubbles (Taylors bubbles ), cell flow (annularslug flow) and annular flow. In their observation, the isolated bubble was called as the small bubble region. They pointed out that the capi llary forces strongly a ffect the flow pattern PAGE 32 32 and heat transfer in the horizon tal annular channel. Coleman a nd Garimella (1999) investigated flow patterns and flow transitions for airwater flow in round a nd rectangular tubes with small diameter ranging from 5.5 to 1.3 mm. They observed bubble, dispersed, elongated bubble, slug, stratified, wavy, annularwavy, and annular flow patterns. They pointed out that the channel diameter and surface tension effects play a signifi cant role in determining the flow patterns and flow transitions in the case of tubes with smaller diameter. Furthermore, their results showed that these effects suppress the stratified flow regime in small cha nnels. They proposed that the aspect ratio, hydraulic diameter, and surface tension ma y be important factors to identify the flow transitions. Triplett et al (1999a) investigated flow pa tterns for airwater mixtures in circular microchannels with 1.1 and 1.45 mm inner diameters and in microc hannels with semitriangular crosssections with hydraulic di ameters of 1.09 and 1.49 mm. Five flow patterns were observed in the experiments: bubbly, churn, slug, and sl ugannular and annular flow except the stratified flow. Flow patterns for airwater flow and steamwater flow in mi crochannels were visualized by Feng and Serizawa (1999 and 2000) and Serizawa et al. (2002) Their flow patterns were dispersed bubbly flow, gas slug fl ow, liquid ring flow, li quid lump flow, annul ar flow, frothy or wispy annular flow, rivulet flow, and liquid droplets flow. The liquid ring flow transited from the slug flow at high gas velo cities and the liquid lump flow deve loped from the li quid ring flow at the increased gas flow rate in liquid ring flow. Th ey noticed that the liqui d lump flow is similar to the wavy flow in a horizontal macrosize tube. Kawahara et al (2002) investigated experimentally nitrogen gaswater flow patterns in a circular mi crochannel with a diameter of 100 m. Different flow patterns at a given flow condition were obs erved: liquid slug, gas core with a smooththin liquid film, gas core with a smooththick liquid film, gas core with a ring PAGE 33 33 shaped liquid film, and gas core with a deform ed interface. Bubbly and churn flow were not observed in the experiments. They also investig ated the ringshaped li quid films first observed by Feng and Serizawa (1999 and 2000). The results indi cated that the ringshaped liquid film and smooththin liquid film are affected by the surface tension rather than the gravitational effect. Kandlikar (2002a) comprehensively reviewed tw ophase flow patterns in minichannels and microchannels. It was noted that the confined flow is similar to the plug flow pattern observed in conventional channels with a larger diameter (> 3 mm). Kawaji and Chung (2004) investigated adiabatic nitrogen gaswater flow in microchannels with hydraulic diameters of 100 m and 96 m. Slug flow was only obser ved in the microchannels. The effect of the gravity effect on the tran sition from symmetry flow to asymmetry flow patterns for condensation flow in miniand microtubes was studied by Li and Wang (2003). They obtained the critical and thre shold values of tube diameter, cD and thD in terms of the Laplace constant as follows: 0.224cLD (25) 1.75th LD (26) For cDD the gravitational effect in the channels ca n be negligible and possible flow patterns are annular, lengthened bubble, and bubble flow. For cDD the surface tension effect is smaller than the gravitational effect. Thus, the fl ow regime may similar to those in macrosize tubes. When ct hDDD the flow patterns are asymmetrical annular and asymmetric elongated bubble flow because both the surface tension and gravitational effects ar e dominant in this criterion. Xu et al. (2005) studied microscale boiling heat transfer using acetone as the working fluid in ten parallel sili con microchannels with hydraulic diameter of 155.4 m. Four flow patterns were observed: liquid plug/vapor slug flow, bubble slug entrained in a liquid plug, PAGE 34 34 paired or triplet bubbles entrained in a liquid plug, and transient a nnular flow. It was investigated that all microchannels repeat the transient fl ow patterns in milliseconds. Revellin et al. (2006) studied diabatic twophase flow patterns for R134a in microchannels with an internal diameter of 0.5 mm as shown in Figure 23. They observed four flow patterns: bu bbly, slug, semiannular, and annular flow. For higher mass fluxes, the an nular flow was quickly investigated due to coalescence among small bubbles. 2.1.3 Flow Regime Map Flow regime maps represent graphically fl ow regimes at different conditions using superficial gas and liquid velocities defined by v vGx j (27) 1l lGx j (28) where tcGmA is the mass flux, x is the vapor quality, tm is the total mass flow rate, and cA is the crosssectional area. Baker (1954), Mandhane et al (1974), and Barnea et al. (1980) proposed experimentally flow regime maps for gasliquid fl ow in horizontal or slightly in clined round tubes. Theoretical models to predict flow pattern transition boundari es for adiabatic twopha se flow in horizontal macrochannels were developed by Taitel and Dukl er (1976) and Weisman et al. (1979). Flow regime map represented by Taite l and Dukler (1976) has b een more commonly used in macrochannels. The horizontal axis of the flow regime map is the Martinelli parameter X as defined by 1/2 l vdpdz X dpdz (29) PAGE 35 35 where ldpdz and vdpdz are the frictional pressure gr adients for the liquid and vapor phases flowing along a channel as 2 221l l l f Gx dp dz D (210) 222v v v f Gx dp dzD (211) with the friction coefficients fo r the liquid and vapor phases are l f and v f respectively. The two vertical positions of the map are the dimensionless parameters, TDK for the stratified flowtowavy flow transition, TDF for the wavyannular and wavyintermittent transitions and TDT for the bubbly flowtointermittent flow tr ansition, respectively, defined by 1/2 2cosvvl TD llvjj K vg (212) 1/2 2cosvv TD lvj F Dg (213) 1/2cosl TD lvdpdz T g (214) where lv is the kinematic viscosity of the liquid and is the inclination angle of the channel to the horizontal. A lthough the flow regime map repr esented by Taitel and Dukler (1976) is widely used to predic t the flow patterns for horizontal flow in macrochannels, it can hardly analyze the flow pattern s in minichannels and microchannels. The reason is that the gravitational effect is dominant in macrochannels a nd has considerably in fluence on the flow patterns. The channel orientation has a significantly effect on th e twophase flow map as well. As aforementioned, the surface te nsion in miniand microchannels, however, is a foremost PAGE 36 36 factor rather than the gravitati onal effect. The channel orientati on and the gravitational effect have no longer an important effect on the twophase flow map in the miniand microchannels due to their smaller diameter. Instead, both cha nnel size and geometry of microchannels affects flow pattern transitions. Co rnwell and Kew (1992) proposed the Confinement Number to represent the restriction to the general flow regimes caused by channel size effects as follows: lv hg Co D (215) They proposed that the Confinement Number can be used to find the transition from isolated to confined bubble regimes and the first confined boiling occurs when 0.5Co. According to Suo and Griffith (1964), buoyancy effects may be negligible for 3.3Co Twophase flow boiling in microchannels has typically a large value of Co, which points out the predominance of the surface tension force over the gravitational force. Flow regime maps for airwater mixtures in small round and rect angular tubes with diameters ranging from 5.5 to 1.3 mm were de veloped by Coleman and Garimella (1999). The liquid and gas superficial veloci ties ranged from 0.01 to 10.0 m/s and from 0.1 to 100 m/s. The flow regime maps showed that flow regime tr ansitions are dominantly affected by the tube diameter. The transition from a plug and slug regime to a disper sed or bubbly regime found at higher superficial liquid velocitie s with the decreased tube diameter. It may be dominant surface effects. At higher superficial vapor velocities, the transition to the annular flow regimes, such as wavyannular and annular flow was observed. At higher and nearly consta nt superficial vapor velocities, the transition from wa vyannular flow to pure annular flow was investigated. Triplett et al (1999b) constructed a fl ow regime map for gasliquid twophase flow in circular microchannels with 1.1 and 1.45 mm inner diamet er and semitriangular microchannels with PAGE 37 37 hydraulic diameters of 1.09 and 1.49 mm. The range of superficial velocities of liquid and gas phases were 0.02 8 m/s and 0.02 20 m/s, respec tively. The comparison of their flow regime with relevant flow transition models and correlations showed poor agreement. Tabatabai and Faghri (2001) proposed a new tw ophase flow regime map to emphasize the surface tension effects in horizontal miniature an d micro tubes. Based on force balance including shear, buoyancy and surface tension forces, a tran sition boundary was proposed. At the flow regime map with transition boundari es, the horizontal axis represen ts the ratio of pressure drop due to surface tension forces to that due to shear forces, which is the sum of the liquid and gas phases defined by surface tension surface tension shear ,ratio lsgsdp dp dz dz p dpdpdp dzdzdz (216) The vertical axis denotes the ratio of superficial gas to superfic ial liquid velocities. The proposed boundaries account for the surface tensiondominat ed regime and sheardominant regime; the surface tensiondominated regime represents bubb ly, plug, and slug flow and the sheardominant regime refers to stratified, wavy, annular, and misty flow. They pointed out that the transition boundary is valid for small diameter tubes of less than about 4 mm. Zhao and Bi (2001) identified flow regimes of cocurrent airwater twophase flow in vertic al triangular channels with hydraulic diameters of 2.886, 1.443, and 0.866 mm. The superf icial velocities of liquid and gas phases ranged from 0.08 to 10 m/s and from 0.1 to 100 m/s, re spectively. The flow regimes indicated that the flow pattern s encountered in the conventional vertical channels, such as dispersed bubbly, slug, churn, and annular flow are observed in the channels with hydraulic diameters of 2.886 and 1.443 mm. The disperse d bubbly flow was not found in the channels PAGE 38 38 having a hydraulic diameter of 0.866 mm. Instead, capillary bubbly flow was observed in the smaller channels. It was found that the transi tion boundary from slug to churn flow and from churn to annular flow in the flow regime maps shifts to the right with the decreased channel diameters. Flow pattern maps for airwate r twophase flow in circular microchannels with a diameter of 20 m were investigated by Serizaw a et al. (2002). They identif ied four flow regimes: bubbly flow, slug flow, liquidring flow, and liquid lump flow except liqui d droplet flow. Kawahara et al. (2002) presented a detailed flow pattern map for a 100 m circular tube with nitrogen gaswater. The superficial velocities of liquid and gas phases were in 0.02 4 m/s and 0.1 60 m/s, respectively. Bubbly and churn flow patterns were not found in their results. They pointed out that the absence of bubbly and c hurn flow may result from the low liquid Reynolds and laminar nature of liquid flow in the microchannels. Chung et al. (2004) developed flow pattern maps for adiabatic gasliquid twophase flow in a square channel with 96hD m and a circular channel with 100hD m. Based on the observed flow patterns (slugring, ringslug, multiple, and semiannular flow), they pointed out that th e transition boundaries were dependent on the microchannel shape. Qu et al. (2004) represented a flow pattern map for twophase flow in rectangular microchannels with 0.4062.032 mm. They showed significant discrepancy between their map and the present microchannel data. Stratified, wa vy, and dispersed flows were never observed in their study. They noted that the discrepancy betw een flow patterns and tr ansition boundaries is due to the significant sensitivity of gasliquid fl ow pattern to working fluid, channel size and geometry for minichannels and microchannel flow s at low superficial velocities. Revellin and Thome (2007) proposed a new diabatic flow pattern map for boiling heat tr ansfer of refrigerant PAGE 39 39 R134a and R245fa in microchannels with diameters of 0.509 and 0.790 mm as shown in Figure 24. They proposed that the transitions are governed by the rate of coalescence. The reason is that small bubbles fast coalesce to form elongated slug ones at higher mass flux. The new diabatic flow pattern map based on a mass flux a nd vapor quality for evaporating flow in the microchannels is characterized by four zones: isolated bubble regime, coalescing bubble regime, annular, and postdryout zone. They noted that the flow pattern map can be used to determine the reasonable operating range of microevaporators. A general flow pattern map in minichannels and microchannels ha s not been entirely constructed due to simultaneous ly considering several dimens ionless parameters such as Re and Rellh vvh lv lv j Dj D (217) 22 and llh vvh lv j Dj D We We (218) 22 and lv lv hh j j Fr Fr gDgD (219) ll j Ca (220) v l j S j (221) where Weber numbers lWeand vWe are the ratios of inertia to su rface tension, Froude numbers lFrand vFrinertia to gravity, Capillary number Caliquid viscous force to surface, and S(slip ratio) vapor inertia to liquid inertia. 2.2 Pressure Drop during Flow Boiling in Minichannels and Microchannels Pressure drop in channel flow is the mo st important parameter to design microscale devices. The overall pressure drop for twophase flow boiling is the sum of the singlephase and PAGE 40 40 twophase flow. The pressure drop in microcha nnels significantly deviates from that in macrochannels because of the channel size, r oughness, and geometry of microchannels. The accurate predictive model to estimate pressure drop, especially a twophase region, is very limited in those channels. The previous studies for pressure drop model including singphase and twophase regions are summarized in this section. 2.2.1 SinglePhase Pressure Drop The singlephase pressure drop over a channel le ngth is affected by the friction effects of the flow: 22huL pf D (222) where the friction factor f for a circular channel is defined by 16 Ref (223) The significant discrepancies between the expe rimental results obtained in microchannels and the predictions of the convent ional theory were reported by se veral researchers (Mala and Li, 1999; Qu et al., 2000; Jiang et al ., 2001), while others have re ported good agreement (Xu et al, 2000; Judy et al., 2002). Mala and Li (1999) studied experimentally si nglephase pressure drop in microtubes with diameters ranging from 50 to 254 m. Their experimental results showed that the experimental pressure drop for small Reynolds number is a pproximately equal to th at predicted by the conventional flow theory; however, as the Reynolds number increases, measured pressure drop is significantly higher than predicted by the Poiseuille flow theory. In addition, they showed that the friction factor is higher than that obtained in the conventi onal theory. It was found that the reason of the large discrepancy is due to the effects of surf ace roughness of the microchannels. PAGE 41 41 Qu et al. (2000) reported pressu re gradients and flow frictions of water in trapezoidal silicon microchannels having hydraulic diameters of 51 169 m. The pressure gradients and friction factors obtained from their experiments in micr ochannels were higher than those estimated by the conventional laminar flow theory. It was found that the surface roughness of the microchannels may affect the higher pressure gradients and fricti on factors. Jiang et al. (2001) investigated experimentally pre ssure drop in microchannels with 0.2 mm wide and 0.6 mm deep. They found that the transition from laminar to turbulent flow in the microchannels occurs at Re600 compared to the conventional channels. Xu et al. (2000) investigated experimentally flow friction of liquid flow in microchannels with hydraulic diameters in the range of 30 344 m at Reynolds number ranging from 20 to 4000. Flow characteristics flow in microchannels agreed with conventional behaviors predicted by NavierStokes equations, although very small channel dimensions caused micro effects. Judy et al. (2002) studied pressure dr op for liquid flow in round and r ectangular microchannels in the range of 15 150 m. They did not find discrepancies be tween the experimental results and macroscale Stokes flow theory for any channel cr osssection, diameter, and material. Wu and Cheng (2003a) investigated fric tion factors of laminar flow of deionized water in smooth trapezoidal silicon microchannels with hydrau lic diameters in the range of 25.0 291.0 m. The experimental results showed that the friction factor is significantly affected by the crosssectional aspect ratio. It was found th at transition from laminar to turbulent flow occur atRe15002000 in the microchannels. They emphasized that the NavierStokes equations are still valid for the laminar flow of deionized wa ter in the silic on microchannels having hydraulic diameter as small as 25. 9 m. PAGE 42 42 Bahrami et al. (2006) studied the influen ce of surface roughness on the laminar, fully developed, incompressible flow in microtubes and proposed a new model considered with roughness effects. They showed that the increased pressure drop is influenced by the effect of surface roughness in microtubes. Steinke and Ka ndlikar (2006) comprehensively reviewed the available literature on singlephase liquid friction factors in microchannels and conducted experiments on square silicon microchannels with 200 m sides and 10 mm in length. They pointed out that the Stokes and Poiseuille flow theories can apply for singlephase flow in microchannel flows. 2.2.2 TwoPhase Pressure Drop The twophase pressure drop consists of tw o components: accelerat ional and frictional pressure drops. The acceleration pressure dr op is caused by the flow acceleration during the phasechange process as (Loc khart and Martinelli, 1949): 22 22 2 in11 11out in out in a voutlout vlinxx xx pG (224) where the void fraction can be eval uated by Zivis correlation (1964) 1 2/31 1v lx x (225) The twophase frictional pressure drop in macr ochannels is calculated by the homogeneous flow model or the separated flow model associated with the twophase multiplier. The twophase multiplier is widely used for predicting the twophase frictional pressure drop in channel flow. The parameter can be determined by Lockhart and Martinelli correlation (1949), MartinelliNelson correlation (1948), or Fr iedels correlation (1979) and is defined by (Carey, 1992) PAGE 43 43 1/2 fr l ldpdz dpdz (226) 1/2 fr lo lodpdz dpdz (227) 1/2 fr v vdpdz dpdz (228) The twophase multiplier of the liquid and vapor phases is given by (Chisholm, 1983) 21 1lC X X (229) 21vCXX (230) where the empirical parameter C depends on the flow region of each phase, laminar or turbulent region as listed in Table 22 and the Martinelli parameter X is estimated by Equation (29). The pressure drop of flow boiling in m icrocha nnels has been studied by many investigators in the past years. Table 23 summarized twophase pressure drop corr elations in the open literatures. Bergles and Dormer (1969) investigated ea rly the pressure drop for subcooled flow boiling of water in horizontal sm all tubes with diameters ranging from 1.57 to 5.03 mm. Bowers and Mudawar (1994a, 1994b, and 1994c) investigated flow boiling pressure drop of R113 in minichannel (2.5 mm i.d.) and microchannel (510 m i.d.) heat sinks. They evaluated the pressure drop in the twophase region using hom ogeneous flow model and the predicted pressure drop was compared very well to the experimental data. Mishima and Hibi ki (1996) studied the frictional pressure drop for airwater mixtures in capillary vertic al tubes with diameters ranging from 1 to 4 mm. They developed a new correlation for the empirical parameter C related with an inner diameter to improve the accuracy of the frictional pressure drop as follows: PAGE 44 44 319211hDCe (231) Triplett et al. (1999b) experimentally investig ated the twophase frictional pressure drop for gasliquid in circular microchannels with 1.1 and 1.45 mm diameter and in microchannels with semitriangular crossse ctions with 1.09 and 1.49hD mm. They showed that for bubbly and slug flow patterns, the twophase friction factor based on ho mogeneous model predicts very well their experimental data with good accuracy For annular flow, however, the homogeneous model and other widely used correlations significantly overpredict the frictional pressure drop. Tran et al. (2000) studied the twophase pressure drop of three refrigerants (R134a, R12, and R113) in round tubes with 2.46 and 2.92 mm inne r diameter and a rectangular channel with 4.06 mm wide and 1.7 mm deep during flow boiling process. They proposed a new frictional pressure drop correlation associated with the effects of surface tension and channel size. The deviation between the predicted values and experimental data was within 20% Lee and Lee (2001a and 2001b) studied flow boiling pressure drop for wate r and airwater mixtures in horizontal small rectangular channels with gaps of 0.4 4 mm an d for R113 in single rectangular channels of 20 mm wide and 0.4, 1 and 2 mm high Several previous correlations developed for both macroand mini/microchannels were evaluated for the fricti onal pressure drop. Large discrepancies except for a modified LockhartMartinelli correlation ta king into account of mass fluxes and channel gaps were investigated by Lee and Lee (2001a ). The modified LockhartMartinelli correlation for the frictional pressure drop was within accuracy of 20% Zhao and Bi (2001) investigated twophase pressure drop of airwater twophase flow in vertical miniature triangular channels with hD = 0.866, 1.443, and 2.886 mm. They po inted out that the LockhartMartinelli correlation with the newly propose d friction factor correlation for singlephase can better predict the twophase pressure drop. Flow boiling pressure drop of water in a single horizontal tube with PAGE 45 45 a diameter of 2.98 mm was evaluated by Yu et al. (2002). They develope d a modified LockhartMartinelli parameter and the Chisholm correlation to predict the twophase pressure drop. The new correlation showed better agreement with the experimental results. Kawahara et al. (2002) studied twophase pressure drop of deionized water and nitrogen gas in a microchannel having a diameter of 100 m. They showed that the twophase fric tion multiplier data were overpredicted by the homogeneous flow model; however, the sepa rated flow model proposed by Lockhart and Martinelli (1949) correlated well within 10% Qu and Mudawar (2003a) developed predictive tools of the twophase pressure drop for water in microchannel heat sinks with 21 parallel rectangular channels (231713 m). They investigated the accuracy of twophase pressure drop using generalized twophase pr essure drop correlations. A ne w correlation accounting for the effects of both channel size and coolant mass velo city was developed. Their results showed that the new correlation was more accurate than prev ious correlations for predicting pressure drop. Qu and Mudawar (2004) compared three widely used twophase pressure drop models and correlations, which are the homogeneous flow m odel, MartinelliNelson correlation model and LockhartMartinelli correlation, to experimental data with their previous geometry. They pointed out that the homogeneous flow model shows poor pr edictions of the experimental data for water; however, the LockhartMartinelli correlation represents the good agreement with the experimental data. Lee and Mudawar (2005a) measured twophase pressure drop for R134a in a microchannel heat sink with 231 m wide and 713 m deep. They showed th at the predictions of the homogeneous flow model, separated flow mode l and correlations are not compared very well to the experimental results. They developed a new correlation which modified the twophase pressure drop multiplier with the effects of su rface tension and liquid viscosity and the proposed correlation showed good agreement with the experimental results. Hwang and Kim (2006) PAGE 46 46 investigated twophase flow pre ssure drop of R134a in microtubes with 0.244, 0.430, and 0.792 mm i.d. They proposed a new correlation includin g the effects of the tube diameter, surface tension effect, and Reynolds numbe r to predict the twophase flow pressure drop in microtubes. They reported that the accuracy of the new correlation is within 8.1% However, the predictive tools of twophase pressure drop to increase accu racy compared with experimental data have been still developed by a number of investigators. 2.3 TwoPhase Flow Boiling in Microchannels Boiling is a phase change process through the liquidvapor interface when vapor bubbles are formed either on a heated surface or in a supe rheated liquid layer near th e heated surface. It is utilized in various energy conversion and heat exchange systems and in cooling system of highpower density electronic devices. The liquidvapo r flow configurations in flow boiling change due to the rapid generation of vapor bubbles alon g the flow direction. In consequence, flow boiling through heated channels effici ently enhances heat transfer. In the modeling of boiling nucleation and bubble growth, the fluid mechanics and heat transfer associated with a liquid thin film in the form of a meniscus and a contact line are closely related to the boiling process in microchannels. The liquid thin film plays an important role in nucleate boiling heat transfer. This liquid f ilm underneath a vapor bubble is called the microlayer. Professor P. Wayner and his coworke rs have reported exte nsive experimental and modeling works in the thin film hydr odynamics, and heat and mass transfer. In general, interfacial phenomena have a signif icant effect on heat transfer characteristics in microchannel heat sinks requiring higher heat fluxes. The interfacial phenomena within the evaporating extended meniscus have been studied by a number of investigators. PAGE 47 47 2.3.1 PhaseChange Phenomena in Meniscus The properties of small liquid systems deviate fr om those of a bulk liquid due to the effects of longrange intermolecular forces. Intermolecula r and surface forces result from the electronic structure of atoms and molecules. These forces cause the adhesion of one substance to another, the cohesion in bulk liquids, the free energy associated with interfaces, and the liquidvapor phase distribution in a closed system. The disj oining pressure is caused by the intermolecular forces between two apolar molecules or atoms. The disjoining pressure and its role in preventing liquid film dryout were first proposed by Derjaguin et al. (1957). At an e quilibrium nonevaporating film re gion, liquid flow is prevented from evaporating by interfacial resistance to mass transfer. They suggested that thin film transport, beyond this film region, is capable of significantly increasing the evaporation rate from capillaries. Derjaguin et al (1965) developed a model incorporating an adsorption isotherm of the disjoining pressure as well as mass transport in capillary tubes. They showed that the mass transport in the liquid thin film near the interlin e dramatically increases the evaporating rate over that estimated solely by diffusion theory. Potash and Wayner (1972) experimentally stud ied the transport processes occurring in an evaporating twodimensional meniscus and adsorbed thin film formed on a superheated flat glass plate immersed in a pool of saturated carbon tetr achloride with Derjagui ns disjoining pressure. They proposed the extended meniscus which consists of three distinct regi ons as shown in Figure 25. The intrinsic meniscus region is domina ted by capillary forces; the equilibrium nonevaporating liquid film region is mainly governed by the disjoining pressure due to intermolecular interactions between the wall and the bulk liquid. In the adsorbed film region, superheated liquid molecules are prevented from evaporating by extremely strong interaction forces between the liquid molecules and the soli d substrate. The evaporating thin film region PAGE 48 48 exists between the intrinsic meniscus and the adsorb ed film region. In this region, the effects of capillarity and the disjoining pre ssure are equally important for fl uid flow and heat transfer. The results showed that fluid flow results from a ch ange in the external meniscus profile and the presence of an adsorbed superh eated film results in a smooth tr ansition between the evaporating and nonevaporating portions of the extended meni scus. Miller (1973) inves tigated the stability of liquidvapor interfaces moving as a result of phase change or mass transfer. His result showed that instability during vaporiz ation of a liquid film may o ccur only when pure liquids of relatively low volatility are vaporized under vacuum, or that the liquid has very low surface tension. Wayner et al (1976) first developed a mathema tical model consisting of an adsorbed nonevaporating layer and an evaporating meniscus w ith a constant surface te mperature as shown in Figure 26. The interline region se parates the two parts and the cont rolling forces are due to the disjoining pressure and capillary pressure. Also, a dispersion based model of the evaporating meniscus near the exit of a capillary tube wa s developed by Wayner (1979). He found that the viscous flow in the liquid thin film near the in terline dominantly affected the entire meniscus profile. The results showed that increasing the in terline heat flux reduces the capillary pressure by increasing the apparent contact angle. Renk and Wayner (1979) investig ated an evaporating ethanol meniscus and proposed a mathematical model of capillary flow heat transfer. They showed that the fluid flow resulting from a change in the meniscus resupplies the liquid evaporated in an evaporating meniscus. In additi on, the reported results showed that the fluid velocity, local heat flux, and local temperatur e differences between liquidvapor interface and surface reach a maximum very near the inter line. Moosman and Homsy (1980) developed a mathematical model to account for the transition phenomena in the horizontal extended meniscus PAGE 49 49 using perturbation theory. It was found that th e maximum heat flux occurs in the transition region with small thermal resistance. Over the years the model proposed by Wayner et al (1976) has been improved to include effects of gravity, different li quidsubstrate combinations and detailed meniscus and film structures, which were verified by extensive experimental data (Sujanani and Wayner, 1992; DasGupta et al., 1993 and 1994; Wayner, 1997). The m odel has also been modified to address the motion of a meniscus and a moving contact line (Wayner, 1994). Wayner (1999) comprehensively reviewed on the modeling through the liquidvapor interface and its application to various physical systems. Mirzamoghadam and Catton (1988) developed a physical model of the evaporating meniscus on an inclined, partially submerged copper plate using an integral approach similar to boundary layer analysis. Swanson and Herdt (199 2) proposed a mathematical model of the evaporating meniscus in a capi llary tube combining the full threedimensional YoungLaplace equation, Marangoni convection, Londonvan der Waals dispersion forces, and nonequilibrium interface conditions. They found that the dimensi onless superheat has no apparent effect on the meniscus profile; the dispersion coefficient, how ever, significantly changes in the meniscus profile. Also, the results showed that the local in terfacial mass flux and total mass transfer rate increase drastically with the dispersion coefficient increased. Begg et al. (1999) proposed a steadystate mathematical model of condensation in a miniature tube to predict the liquid film prof ile along the condenser and the length of the twophase flow region. The model incorporated th e liquidvapor frictional interaction and surface tension gradient. They pointed out that pres sure drop in the vapor phase is insignificant compared to the pressure drop in the liquid phase. PAGE 50 50 2.3.2 Simple Model for TwoPhase Flow Boiling Many investigations of micro heat pipes have been conducted to predict heat transfer during evaporation from the extended meniscus based on Wayners model. Stephan and Busse (1992) developed a mathematical model for the radial heat transfer of grooved heat pipe evaporator walls. They pointed out that th e assumption of an in terface temperature Ti equal to the saturation temperature Tsat of a vapor phase may lead to a large overprediction of the radial heat transfer coefficient. Longtin et al. (1994) developed a onedimensional model of the evaporation and adiabatic section to yield pressure, velocity, a nd film thickness along the length of a micro heat pipe. The proposed model included interfacial and vapor shear stress terms. They found that the pressure drop in th e liquid phase is several times gr eater than that of the vapor phase. Peles and Haber (2000) proposed a simple onedimensional model of twophase flow boiling and heat transfer in a single triangular microchanne l. The model was developed by assuming the constant vapor pressure along the channel and the constant vapor temperature. Sartre et al. (2000) studied the effect of interfacial phenomena on evaporative heat transfer in micro heat pipes. They developed three coupled models, solving the mi croregion equations, the twodimensional wall heat conduction problem a nd the longitudinal capillary twophase flow. The results showed that the major part of the tota l heat input in the evaporator section is in the microregion. Both the apparent contact angle an d the heat transfer rate in the microregion increased with increasing wall superheat. Qu and Ma (2002) investigated the polarity effects of working fluids in a capillary with constant wall temperature. They showed that for strong polar working fl uids, the length of the evaporating interfacial region is much longer compared with those of other apolar working fluids. However, their model did not include the thermocapillary effect though nonisothermal interfacial conditions. Park et al. (2003) developed a mathematical model to describe the PAGE 51 51 transport phenomena for a thin film region on a microchannel with a nonpolar liquid under isothermal conditions. They applied the slip velo city conditions to an evaporating microchannel filled with water under very simplified conditions Thermocapillary effects were not considered in their model. They pointed out that the gradient in the liquid pressure had a significant influence on the thin film pr ofile. In addition, it was found th at the length and the maximum thickness of the thin film decrease exponentially as the heat flux increases. Wee et al. (2005) developed a mathematical model to predict th e micro/nanoscale fluid flow and heat/mass transfer phenomena in an ev aporating extended meniscus unde r nonisothermal interfacial conditions. The model considered the polarity eff ect of working fluids, a slip boundary condition on a heated wall, and thermocapillary stresses at the liquidvapor interface. The results showed that for a polar liquid, the transi tion region of the evaporating meni scus is longer than that of a nonpolar liquid. The strong polar at traction with the solid wall le d to decreasing the evaporative heat transfer flux. The results used the slip boundar y condition increased ev aporative heat and mass flux and dropped the liquid pressure grad ients and viscous drag at the wall. Jacobi and Thome (2002) developed a simple heat transfer model accounting for thinfilm evaporation of elongated bubble flows in microchannels. They pointed out that thinfilm evaporation is dominant heat tr ansfer mechanism. The proposed model predicted heat transfer coefficient in microchannels well. However, the appropriate nucleation radius and initial thin film thickness should be known befo re predict heat transfer using their model. Thome et al. (2004) proposed a threezone flow boiling model in microchannels. The model analyzed the temporal variation of the heat transfer coeffi cient in a liquid slug, evaporating elongated bubble and vapor slug. They noted that heat transfer in the thin film of the elongated slug bubble is PAGE 52 52 dominant and that of the vapor slug can be ne gligible. In addition, it was found that bubble frequency strongly affects on heat tr ansfer behavior in microchannels. Although extensive research of evaporating nonpolar thin film s has been performed due to the simplicity of modeling van de r Waals intermolecular forces, co mplete understanding of heat transfer phenomena through the liquidvapor interface has not been achieved yet. 2.3.3 Numerical Simulation of TwoPhase Boiling Numerical simulation of twophase boiling in cluding an analytical model has been extensively studied by several res earchers in the past decades. A variety of numerical techniques for direct numerical simulation of these problem s include the volumeoffl uid (VOF), levelset, phase field, front tracking met hods, and sharpinterface method. Numerical simulation of bubble growth in nucle ate boiling using numerical grid generation technique was developed by Lee and Nydahl (1989). The model in cluded microlayer evaporation but assumed negligible free convection, and a specified bubble and micr olayer shape. They concluded that microlayer evap oration contributes nearly 90 % of the energy for the bubble growth at 8.5 C and 1 atmospheric pressure. Mei et al. (1995) studied numeri cally the vapor bubble growth rate and associated thermal fields in heterogeneous pool boi ling. They constructed four dimensionless parameters to characterize the vapor bubble grow th rate and the temperature distribution of the solid and the liquid. They assumed the liquid microlayer between a vapor bubble and a solid heating surface as a simple wedge shape. Son and Dhir (1997) numerically simulated sa turated film boiling on a horizontal surface using a numerical grid generation method. The results showed that the heat transfer coefficients decreased with increase in wall superheat and at higher pressures, the growth of the interface PAGE 53 53 slowed down. Son and Dhir (1998) used the le velset method to numerically simulate film boiling near critical pressures. Welch (1998) carried out a local simulation of axissymmetric vapor bubble growth. The numerical model used an interface tracking met hod in conjunction with a finite volume method on a moving unstructured mesh. Son et al ( 1999) numerically simula ted a twodimensional growth and departure model of single vapor bubbles during nucleat e pool boiling. They used the levelSet method to capture the liquidvapor interface. A Numerical model of fluid fl ow with phase change was developed by Welch and Wilson (2000). The liquidvapor interface was tracke d by the volume of fluid (VOF) scheme incorporated with mass transfer and surface tension. They showed that the numerical simulation using the VOF technique provides accurate calculation of horizon tal film boiling. Ye et al. (2001) developed a direct numerical simulation method to investigate single bubble deformation and phase change using the sharpinterface method to capture the liquidvapor interface with large density ratio. They po inted out that the sharpinterface approaches successfully simulate bubbly dynamics for density ratio of 1600 or higher. Mukherjee and Dhir (2004) studied threedimensional cases for merger and departure of multiple bubbles in nucleate pool boiling using th e levelset technique. They showed that merger of multiple bubbles significantly increases the ov erall wall heat transfer. The reason resulted from trapping of liquid layer between the bubbl e bases during merger an d caused by liquid flow toward the wall during contracti on after merger. Mukherjee and Kandlikar (2005) performed the numerically simulation of the vapor bubble grow th during flow boiling through a microchannel. The levelset method based on the SIMPLER tec hnique was used to capture the liquidvapor interface. They showed that the bubble growth ra te increases with the incoming liquid superheat PAGE 54 54 while this rate decreases with the Reynolds number. The effect of gravity on the bubble growth rate was observed to be negligible for the case with Tin = 107 C. It was found that the bubble growth obtained from the numeri cal simulation is similar to the experimental data. Threedimensional numerical simulation of single bubble with different contact angles during nucleate pool boiling is conducted by Mukherjee and Kand likar (2007). The SIMPLER method was used to calculate the incompressible NavierStokes equations and the interface location was solved by the level set scheme. It was found that the numerical results are in good agreement with experimental observations. Geisler (2007) simulated natural convection boiling in vertical channel with the channel gap of 7, 0.5, and 0.3 mm using the VOF me thod. The numerical simulation using the VOF method successfully calculated natural c onvection boiling without phase change. Only a few investigations for numerical simu lation of flow boiling in a microchannel have been reported in the literature. Furthermore, there is no threedimensi onal numerical simulation dealing with periodic flow boiling in a microchannel in open literature. PAGE 55 55 Figure 21. Schematic represen tations of flow regimes observe d in horizontal, cocurrent gasliquid flow (Carey, 1992). Figure 22. Flow patterns in small parallel rectan gular channels (Cornwell and Kew, 1992). PAGE 56 56 A B C D E F G Figure 23. Twophase flow pa tterns and transitions of flow boiling in microchannels. A) Bubbly flow, B) Bubbly/slug flow, C) Sl ug flow, D) Slug/semiannular flow, E) Semiannular flow, F) Wavy a nnular flow and G) Smooth a nnular flow (Revellin et al., 2006). PAGE 57 57 Figure 24. New diabatic coalescing bubble map for evaporati ng flow in circular microchannels (Revellin and Thome, 2007). Equilibrium nonevaporating thin film (Adsorbed thin film) Evaporating thin film Intrinsic Q Transverse liquid flow Figure 25. Flow region s through a liquid thin film in a microchannel. PAGE 58 58 0 x Figure 26. Interline junction of vapor, ad sorbed evaporating liquid thin film and nonevaporating adsorbed liquid film (Wayner et al., 1976). PAGE 59 59 Table 21. Classificatio n of channel dimensions (Kandlikar and Grande, 2003) Channel Hydraulic diameters Conventional Channels 3hDmm Minichannels 2003hmDmm Microchannels 10200hmDm Transitional Channels 0.110hmDm Transitional Microchannels 110hmDm Transitional Nanochannels 0.11hmDm Molecular Nanochannels 0.1hDm Table 22. Empirical parameter C in LockhartMartinelli correlation (Chilsholm, 1983) Liquid phase Vapor phase C laminar laminar 5 turbulent laminar 10 laminar turbulent 12 turbulent turbulent 20 PAGE 60 60 Table 23. Summary of twophase frictional pr essure drop correlations Author Frictional pressure drop LockhartMartinelli (1949) ,2 2 2 2 2, 0 ,21eoutx le frle eouthlL fGx pd x xD 2 21 1lC X X ; constant C from Table 22 Friedel (1979) ,2 2 2 2, 0 ,2eoutx lo frloe hleoutfGL pd x Dx 2 23 1 0.0450.0353.24lo A A A FrWe; 2 2 11lvo ee vlo f Axx f ; 0.224 0.78 21eeAxx 0.910.19 0.7 31lvv vllA ; 2 2 hG Fr gD ; 2 hGD We ; 11ee vlxx Chisholm (1983) ,2 2 2 2, 0 ,2eoutx lo frloe hleoutfGL pd x Dx 0.5 220.5111ex loeeBxx ; 0.50.5 lv vl Collier and Thome (1994) 2 22 2,2 1 2eout l fr hl vlfGL x p D ; 20.003 f Mishima and Hibiki (1996) ,2 2 2 2 2, 0 ,21eoutx le frle eouthlL fGx pdx xD 2 21 1l vvvvC X X ; 0.319211hDCe;[]hDmm Trans et al. (2000) ,2 2 2 2, 0 ,2eoutx lo frloe hleoutfGL pd x Dx 0.875 220 8 7 51 7 5111lo confe e eCNxxx ; 2vo lodpdz dpdz ; lv conf hg N D PAGE 61 61 Table 23. Continued Author Frictional pressure drop Lee and Lee (2001b) ,2 2 2 2 2, 0 ,21eoutx le frle eouthlL fGx pdx xD 2 21 1lC X X ; 20.7266.18510ReloC for laminar liquid phase turbulent vapor phase Yu et al. (2002) ,2 2 2 2 2, 0 ,21eoutx le fr le eout hlL fGx pd x xD 21.9lX for laminar liquid phase and turbulent vapor phase 0.5 0.1 0.5Re 1 18.65 Revve le lx X x Qu and Mudawar (2003a) ,2 2 2 2 2, 0 ,21eoutx le fr le eout hlL fGx pd x xD ; 2 21 1lC X X ; 0.3192110.004180.0613hDCeG for laminar liquid phase and laminar vapor phase 0.5 0.50.51lev velx X x Lee and Mudawar (2005a) 2 21 1lC X X ; 0.0470.602.16Revv loloCW e ; 0.250.231.45Revt loloCW e 2Re;hlh lo lo lGDvGD We Hwang and Kim (2006) ,2 2 2 2 2, 0 ,21eoutx le frle eouthlL fGx pdx xD 1 23 0ReC C C loconfCCXN ; Cfrom Table 5 in Hwang and Kim (2006) lv conf hg N D PAGE 62 62 CHAPTER 3 SEMIANALYTICAL EVAPORATION MODEL A physical and m athematical model of an annul ar liquidvapor flow pattern in a circular microchannel has been developed. The schematic of the physical model is given in Figure 31. In this model, an annulus flow pattern is assumed which has been observed in most microchannel twophase flows. The liquidvapor twophase flow is separated by a sharp interface. The conservation of mass, momentum and energy e quations with the axisymmetric cylindrical coordinates as shown in Figure 32 were developed and then so lved for the liquid and vapor phases separately. Then the liquidvapor interface was traced by maintaining continuity of mass flux, heat flux and shear stress between the two phases. The vaporiz ation mass flux at the interface was determined by the kinetic theory of gases. In a microchannel, the pressure jump condition due to the surface tension is an importa nt effect; it, therefore, is included for the interface force balance. 3.1 Assumptions The basic assumptions adopted in th e analysis are listed as follows: An axisymmetric condition is used in the cylindrical coordinates. Continuum is assumed, and therefore the nos lip condition is applied at the tube wall. Liquidvapor interface is sharp and smooth. The radial pressure gradient is negligible. The inertial terms can be neglected for liquid films due to lower mass flow rate. Surface tension effects are included in the liquid thinfilm region. Heat transfer in the liquid th in film is dominated by conducti on in the radial direction. There is no temperature gradient in the vapor phase in the radial direction. PAGE 63 63 3.2 Mathematical Modeling 3.2.1 Mass and Energy Balance The mass balance and energy balance between 0z and an upstream local point, zz in the twophase zone as shown in Figure 32 yield, .02,R l ll evap RmrurzdrMz (31) where ,0 lm is the liquid mass flow rate at 0z evap M z is the total mass evaporation rate between 0 z and zz at the liquidvapor interface defined by evap lvQz Mz h (32) and takes a positive quantity. lvh denotes the latent h eat of vaporization. Qz is the total heat transfer rate taking place at liquidvapo r interface due to phasechange between 0 z and zz and can be expressed by the following thermal energy balance ,, 0 0()2()z wl p l ll p l l zzQzRqzdzmcTmcT (33) where ()wqz is the imposed local heat flux on the so lidliquid interface at any location and a negative value. Note that wiTT and therefore ()Qzwill be a positive value. For constant heat flux, the total heat transfer rate can be defined by ,, 0 0()2z w lplllpll zzQzRqdzmcTmcT (34) Based on the assumption of dominant conduction in the radial direction, the local heat flux at the wall, for constant wall temperature, can be represented by the following equation, ,ln/wi wzlTT qk RRR (35) PAGE 64 64 where iT is the local temperature of the liquidvapor interface and is the liquid film thickness. The following differential equation for Qz can be obtained from Eq. (33) ,2()wl p l ldQ d R qzmcT dzdz (36) Substituting Eq. (35) into Eq. (36), Eq. (36) becomes the follows: ,2 ln/wi ll p l lTT dQ d km c T dz dz RR (37) For constant heat flux, the gradient of to tal heat transfer rate can be written ,2wl p l ldQ d qRmcT dz dz (38) 3.2.2 Liquid Pressure Gradient Neglecting the inertial term due to lower mass flow rate, the axial momentum equation for the liquid phase can be written as, 11ll ldpu r dzrrr (39) The boundary conditions of the solution to Eq. (39) at the solid wall and the liquidvapor interface are r=R=0 lu (310) ,1 li li rR liud T d rd T d z (311) The boundary condition of Eq. (311) indicates that the liqui d velocity gradient at the interface is affected by th e interfacial shear s tress and the su rface tension gradient due to surface tension effects. Th e solution of Eq. (39) with boundary conditions of Eqs. (310) and (311) gives th e velocities of the liquid phase and the liquidvapor interface, PAGE 65 65 2 121 () ln 4l l ldp urrCrC dz (312) 2 ,121 ln 4l li rR ldp uRC RC dz (313) with 2 1, 2 2111 2 1 ln 4li li lli l ldp dT d CR R dz dTdz dp CRCR dz (314) Substituting Eq. (312) into Eq. (31) and carrying out the integration, the liquid pressur e gradient with respect to th e axial direction can be determined, 7 ,0, 661()1li lli lvlidp C dT Qz d m dzChCdTdz (315) where 2 2 3 4 2 2 5 42 42 65 3 2 7541 ln 22 ln 11 ln ln 22 42 8 2l l l l lR CRR CRR CRRRR CRRRCRC CRCRC (316) 3.2.3 Vapor Pressure Gradient The vapor pressure gradient ca n be written from the axial momentum balance on the vapor phase as follows: 2 2 ,,1 2v vvvvvvli v vvlividp uARfuu dzAz R uuv (317) PAGE 66 66 We assumed that vapor follows the idealgas law to account for the compressibility of vapor flow, v v g v p R T (318) The vapor density gradient for the axial direction is then given as, 211vvvv gvvddppdT dzRdzTTdz (319) The bulk temperature of saturated vapor is re lated to the saturated vapor pressure. Thus, the bulk temperature gradient of saturated va por can be written by the ClausiusClapeyron equation, 2 g v vv vlv R T dTdp dzdzph (320) 3.2.4 Interfacial Shear Stress The liquid interfacial shear stress,,li is affected by the friction force due to the velocity difference between the liquid and vapor phases and the shear stress force caused by mass flux exchange of evaporation. It is determined base d on the relationship suggested by Blangetti et al.(1980) 2 ,,, ,1 2livvvlivvlivi f uuuuv (321) where v f is the Fanning friction factor for a laminar flow with a smooth surface and is given by, 16 Rev vf (322) and the Reynolds number, Rev, is defined as, ,2 Revvli v v R uu (323) PAGE 67 67 3.2.5 Liquid Temperature Variation The liquid bulk temperature, lT, is defined by, 2222()R l R lrTrdr T RR (324) ()lTr in the right term of Eq. (324) is the local liquid film te mperature profile given by the radial conductiondominant temperature dist ribution in a round tube shown below, () ln lnwi liTT r TrT R R R (325) Substituting Eq. (325) into Eq. (324) and integrating th e re sult, the liquid bulk temperature yields 2 2 2 22 ln 11 1 ln 22 4wi liTT TzT RR RR R RR R (326) The bulk temperature gradient of liquid with re spect to the axial dire ction can be obtained below, 2 2 2 21 ln/ 21 ln 224liwidTdTdTdT dzdzdzdz RR R RR R RR (327) Knowing the interfacial temperature and the liquid film thickness, the wall temperature for constant heat fl ux can be obtained by ln/wiw l R RRz TzTzq k (328) PAGE 68 68 3.2.6 Pressure Difference at the LiquidVapor Interface The pressure difference between the liquid and vapor phases is due to the capillary effect and the disjoining pressure, vlcdd p pppKp (329) where K is the local curvature of the liquidvapor interface defined by, 3/2 1/2 22 2 21 11 ddd K dzdzRdz (330) and d p is the disjoining pressure defined by the in termolecular force betw een the liquid and the solid. For nonpolar liquid, the disjoining pressu re can be written as (Wayner et al, 1976), 3dA p (331) where A is the dispersion coefficient and takes to be positive. However, pure water is strongly polar liquid and then the disj oining pressure can be define d by (Holm and Goplen, 1979) lnB dlgipRTA (332) with 1.49 A and 0.0243 B The disjoining pressure decrea ses rapidly as the liquid film thickness increases. The gradient of the disjoini ng pressure for the nonpolar and polar liquids are 43ddp Ad dzdz (333) lnlgi B di lgBRT dp dT d RA dz dzdz (334) 3.2.7 Liquid Film Thickness Differentiating Eq. (330) with Eq. (329), the thirdorder differential equation for the liquid film thickness PAGE 69 69 1 2 2 32 32 2 2 21 2 3/2 2 2 231 1 1 1 11 1vlddddd dzdzdzdz dddd RR dzdzdzdz dpdpdp d dzdzdzdz dd dzRdz 20i idT d dTdz (335) To obtain the liquid film thickness, new variables are defined by: and dd dzdz (336) Combining Eq. (336) and Eq. (335), the firstorder differential equation for Eq. (335) can be obtained as follows: 1 22 122 3/2 2 231 1 1 1 11 1vld i id dz RR dpdpdp dzdzdz dT d RdTdz (337) 3.2.8 Interfacial Temperature The evaporation of the liquidvapor interface is related to the temperature difference and pressure difference at the liquidvapor interface. Using kinetic theory, the evaporative mass flux through the interface is defined by the temperature and pressure jumps at the liquidvapor interface as follows (Wayner, 1991): evapiv lvmaTTbpp (338) with PAGE 70 70 1/2 1/22 22 2 22vlv uiuvi lv uiui p Mh M a R TRTT Vp M b RTRT (339) where is the accommodation coefficient, M is the liquid molecular weight, u R is the universal gas constant, iT is the interfacial temperature, v p is the saturated pressure at vT lvh is the latent of vaporization per unit mass at iT and lV is the liquid molar volume at iT The interfacial heat flux due to heat conduction through the liquid thin film is equal to the wall heat flux; 22we v a p l vqRmhR The interfacial temperature, iT for constant wall temperature, can be obtained by combining Eqs. (35), (329), (332), and (338) as follows: lnlwv i B llgkTGaTbK T kGabRA (340) with ()ln/lvGhRRR For constant heat flux, the interfacial heat flux approximately equals the wall heat flux. Th erefore, the interfac ial temperature is lnw v lv i B lgq R aTbK hR T abRA (341) 3.2.9 Mean Velocities of Liquid and Vapor Phases For steady state evaporation, total mass flow rate at any position in the evaporation region must be constant which requires, ,, tlzvzmmm (342) where ,lzm and ,vzm are the liquid and vapor mass flow rate s, respectively, in the evaporation region defined by, PAGE 71 71 ,, 0 0 0()lzllllll lvQz muAuA h (343) ,, 0 0 0()vzvvvvvv lvQz muAuA h (344) where 2 vAR is the crosssectional area of the vapor, 2 2 lARR is the crosssectional area of the liquid. The mean velocities of liquid of vapor phases are 2 2 lz lz lm u RR (345) 2 vz vz vm u R (346) The evaporative mass flux satisfies the mass continuity at the liquidvapor interface ,,evapllivvimvv (347) where ,liv and ,viv are the blowing velocities due to vapo rization at the liquidvapor interface and take a negative sign. The evapor ative mass flux can be defined by 1 2evap lvdQ m R hdz (348) The blowing velocity can be obtained from Eq. (347) and Eq. (348) as follows: ,1 2vi vlvdQ v R hdz (349) 3.3 Solution Procedures The fluid flow and heat transf er of microchannel flow with constant wall temperature and constant heat flux have been investigated by so lving the eight firstorde r differential governing equations, Eqs (37), (315), (317), (319), (320), (336), and (337), including eight variables: PAGE 72 72 Q, l p v p v vT To solve the governing equations the eight boundary conditions at 0z are needed as follows: 0 (350) 0 and 0 (351) 0Q (352) ,vvoTT (353) vsatvo p pT (354) 0 lvd p pp R (355) v v g v p R T (356) The gradient of addition heat Eq. (38) instead of Eq. (37) has been solved with Eqs. (315), (317), (319), (320), (336), and (337) for constant heat flux. The governing equations with the boundary conditions given by Eqs. (350) (356) have been simultaneously solved by using fourth order RungeKutta method. With additional conditions, lvdpdp dzdz, 0ddp dz, 0 dQ dz, 0idT dz, 0vdA dz, ,0viv at 0z the initial liquid film thickness0 has been obtained before starting the numerical proc edure as shown in Figure 34. For constant heat flux, the constant wall heat flux along the circular tube has been obtained from the energy balance as shown in Figure 33. HSHS w tqA q NDL (357) where N is the total number of the microchannel, H Sq is the heat flux imposed on the heat sink, H SA is the heated surface area, and tL is the total length of the microchannel. The wall PAGE 73 73 temperature wT at z=0 is an unknown variable and should be first obt ained before solving the governing equations. Therefore, th e initial wall temperature at z=0 has been obtained from the energy balance for the unit channel with the saturated temperature s atT ,02HS wsat cqW TT hR (358) where W is the width of unit channel, ch is the heat transfer coefficient between the channel walls and the fluid. In this st udy, the heat transfer coefficien t correlation proposed by Kandlikar and Balasubramanian (2004) was used to obtai n the initial wall temperature because the correlation showed good agreement with experi mental data of water flow boiling in microchannels. The heat sink temperature can be obtained from onedim ensional heat conduction between the channel bottom and the heat sink ,0 H Sw HSw sqH TT k (359) wherewH is the distance from the heat sink botto m wall to the microchannel bottom wall and s k is the thermal conductivity of the heat sink. The liquid mass flow rate at z=0 for constant wall temperature and heat flux can be obtained from the mass balance based on Eqs. (343) and (344) 2,evap lol zL lvQ mm h (360) where 2,lLm should be equal to total mass flow rate tm. The liquid velocity at the end of twophase region is 22t l zL lm u R (361) The mean velocity of vapor at z=0 can be obtained from the mass balance as follows: PAGE 74 74 ,0 ,0 2 0v v vm u R (362) The disjoining pressure,d p can be estimated from Eq. (332) and the blowing velocity, ,viv, has been given by Eq. (349). The liquid temperature, lT, given by Eq. (326) has been solved with the liquid f ilm thickness the wall temperature wT, and the interfacial temperature iT. The interfacial temperature of Eq. (340) or Eq. (341) for the constant wall temperature and the constant heat flux has been estim ated by us ing the root finding method at every point along the axial direction until th e tolerance approaches to 6110 After obtaining the liquid film thickness and the total heat transfer rate Qz, the mass flow rates and velocities of the liquid and vapor phases have been obtained from Eqs. (343) (344), and Eqs. (345) (346), respec tively. The previous iteration values of firstorder derivatives have been used to obtain the values of firstorder derivatives for every marching step. The evaporative mass flux evapm and the blowing velocities ,viv with Eq. (37) or Eq. (38) have been solved from Eq. (348) and Eq. (349), respectively. The solution procedure has b een itera ted until the film thickness approaches to the radius of the microchannel or the mass balance is satisfied as shown in Figure 34. 3.4 Results and Discussion A numerical code was developed to implemen t the physical model and was used to obtain the solutions. The mainly objectives of this st udy were to investigate the heat transfer and pressure drop in microchannel heat sinks with boiling. Water used as the working fluid. The presented results in this study we re compared with the experime ntal data conducted by Qu and PAGE 75 75 Mudawar (2004). Numerical solution s were obtained for annular film evaporation of water as the working fluid with constant wall te mperature and constant heat flux. 3.4.1 Constant Heat Flux The twophase flow and heat transfer char acteristics were numerically studies in a microchannel with a hydraulic diameter of 348 m. For different heat fluxes, the results of the liquid film thickness, the pressu re drop of the liquid and vapor phases, and the heat transfer coefficient at the liquidvapor interface were obta ined at different inle t temperatures with Tsat = 380.73 K: Tin = 30 C and 53.84010/tmk g s. The effects of heat flux on the liquid film th ickness are shown in Figure 35. The results show that the liquid film thickne ss initially increases along the axial directi on and then sharply increases. The main reason is the increased total heat transfer rate through the liquidvapor interface with the increased heat flux. It leads to increasing the liquid pressure gradient which results in increasing the liqui d film gradient given by Eq. (337). As the heat flux increase, the evaporative m ass flux at the liqui dvapor interface sharply increas es. It causes the liquid film thickness to drastically increase. The length of the liquid film thickness becomes smaller with increasing the heat flux. It is because of higher heat flux resulted in highe r evaporation rates. The effects of heat flux on the evaporative mass flux are presented in Figure 36. It is found that the evaporative ma ss flux increases as the imposed heat flux increases. The evaporative mass flux sharply increases at wq = 1.043MW/m2. It leads to decreasing the evaporating length of the liquid thin film and the liquid thin film thickness. As the total heat transfer rate through the liquidvapor interface increases with the imposed heat flux, the liquid pressure gradient increases to maintain the mass fl ow rate in the liquid film. It leads to sharply PAGE 76 76 increasing the liquid film th ickness and, given by Eq. (348), the evaporative mass flux affected by the liquid film thickness and total heat transfer gradient becomes larger. Figure 37 shows the gradient of liquid pressu re for different heat fluxes. The maximum liquid pressure gradient increase s as the heat flux increases. Th e reason is that more liquid mass flow rate from the bulk region to the evaporating region is required to maintain the mass balance due to the increased evaporative mass flux. The larger liquid pressure gradient provides more mass flow into the evaporating thin film region to be continuously va porized. Also, it is found that the liquid pressure gradient does not change as the capillary pressure increases. It means that the heat transfer through the liquidvapor interface is dominant in evaporating thin film region rather than in the extensive meniscus region gove rned by the capillary pressure only. The results are in good agreement with ones of Wee et al. (2005). The pressure difference of the liquid and vapor phases along the axial direction is shown in Figure 38. The magnitude of the liquid and vapor pressure increases along the axial directions. The liquid pressure is affected by the interfacial shear stress and the surfa ce tension effects as indicated in Eq. (315). The liquid pressure sharply increases in the evap orating region and remains unchanged in the extensive meniscus re gion as shown in Figure 38A and the vapor pressure sharply increases at the end of two phase region as illustrated in Figure 38B. The sharply increased vapor pressure results from th e interfacial shear stre ss as indicated in Eq. (317). Although the vapor pressure change is very sm aller than that of the liquid pressure, the vapor pressure may not be ignored due to the effects of the interfacial shear stress. The results point out that the effect of the liquid pressure plays a significant role in the mass flux exchange through the liquidvapor interface rath er than the vapor pressure. PAGE 77 77 The mean velocity variations of the liquid and vapor phases are shown in Figure 39 and Figure 310. The mean velocities of the liqui d and vapor phases are obtained from Eqs. (345) and (346). While the mean velocity of the liquid phase approach es the inlet velocity, 0.420 m /s, as illustrated in Figure 39, the vapor mean velocity abruptly drops as the liquid film thickness sharply increases along th e axial direction as shown in Figure 310. The magnitude of mean velocities of liquid and vapor phases incr eases slowly along the ax ial direction. They increase sharply toward the extensive meniscus region. The reason is that the liquid film thickness drastically increases in this region. It leads to the fact that the magnitude of the vapor velocity is higher than that of the liquid velocity. The interfacial shear stress variations at different heat fluxes are s hown in Figure 311. The results show that the interfacial shear stress increases as the h eat flux increases. The interfacial shear stress mainly affects the vapor pr essure gradient as indicated in Eq. (317). The abruptly increased in terfacial shear stress affects the liqui d film shape and results in the increasing vapor pressure. Figure 312 shows the interfacial temperature variations with diffe rent heat flux. The interfacial temperature is mainly affected by the capillary pressu re and disjoining pressure based on Eq. (341) and decreases along the axial direction due to the increas ed liquid film thickness. It is caused by the increasing capillary pressure an d the decreasing disjoining pressure which result in the sharply dropping interfacial temperature. Figure 313 shows the comparison of the wall temperature obtained from the numerical results to the experimental results (Qu and Mudawar, 2004). The calculated wall temperature was used to the initial wall temp erature at z = 0. The discrepancies between the predicted values and the experimental results of heat sink bottom temperature at 287.9/qWcm and PAGE 78 78 2240/Wcm are 0.56 % and 0.17 %, respectively. Th e numerical results using the predicted initial wall temperature may be valid. The effects of heat flux on the local heat tran sfer coefficients are shown in Figure 314. The local heat transfer coefficient can be obtained by ,ln/l czk h RRR (363) The result shows that with increasing heat flux, the maximum value of the heat transfer coefficient increases. The local heat transfer coe fficient is larger in the evaporating thin film region, while it becomes smaller as the liquid film thickness is larger. The reason is that the thermal resistance increases as the liquid film thickness increases. The liquid film thickness decreases in the evaporating thin film region with increased wall temperature resulted from the inlet temperature increased. It leads to decreasing the thermal resistance which causes higher heat transfer coefficient. The results point out that most eva poration through the liquidvapor interface occurs in the evaporat ing thin film region rather than in the meniscus region. Figure 315 presents the average heat transfer coefficient with differ ent heat flux. It is shown that as the heat flux increases, the averag e heat transfer coefficient increases. The reason is that the increased heat flux results in incr easing the evaporative mass flux at the liquidvapor interface. The increased evapora tive mass flux decreases the liquid film thickness and the length of the evaporating region. Thus, the change of the liquid film thickness caused by the imposed heat flux leads to the change of the heat transfer coefficient. Figure 316 shows the heat transfer coefficient with thermodynamic equilibrium quali ty. The result shows that the heat transfer coefficient increases with increasing quality. The current results are in good agreement with the heat transfer coefficient correlation proposed by Kandlikar and Balasubramanian (2004), while being in poor agreement with the experimental results by Qu and Mudawar (2004). Their results PAGE 79 79 showed that the heat transfer coefficient de creases with increasing thermodynamic equilibrium quality. They considered a larg e number of liquid droplets entrai ned into a vapor core, while a pure vapor phase without liquid droplets was consid ered in the current model. The entrainment of the liquid droplets into the vapor core may increase the liquid film thickness. The increased liquid film thickness due to the entrainment of the liquid droplets in the vapor phase causes thermal resistance to increase and then it leads to decreasing the h eat transfer coefficient with increasing thermodynamic equilibrium quality. 3.4.2 Constant Wall Temperature The transport phenomena through the liquidvapor interface in a microchannel having a hydraulic diameter of 348 m have been numerically investigat ed for different wall temperatures at Tsat = 380.73 K and 53.82310/tmk g s. For water as the working fluid, the effects of constant wall temperature on the liquid film thickness are shown in Figure 317. The liquid film thickness initially increases along the axial direction and then sharply increases. It is mainly due to the larger gradient of the liquid pressure based on Eq. (337). The larger gradient of the liquid pressu re causes the liquid film thickness to increases to maintain the mass balance from the bulk liquid into the liquid film. The length of the liquid film thickness decreases with increasing th e wall temperature. The results are primarily due to the higher wall temperature resulted in high er evaporative mass flux. It leads to decreasing the length of the evaporating thin film. The effects of constant wall temperature on the evaporative mass flux are presented in Figure 318. The evaporative mass flux can be obtained from Eq. (348). The maximum evaporative m ass flux exists within the evaporating thin film region. The maximum location approximately corresponds to the location of th e maximum liquid pressure gradients shown in PAGE 80 80 Figure 319. The evaporation through the liquidvapor interface is facilitated by the increased wall temperature. With increased evaporative ma ss flux, more liquid mass flow is required to maintain the mass balance in the liquid film. Figure 319 shows the gradient of liquid pressure based on Eq. (315) for different wall tem perature. It is shown that the maximum liquid pressure gradient increases with increasing the wall temperature. The reason is that as the evap orative mass flux increases by the increased wall temperature, the liquid pressure gradient become s larger to provide more mass flow into the liquid film. The pressure difference of the liquid and vapor phases along the axial direction is shown in Figure 320. The magnitude of the liquid and vapor pressure increases al ong the axial directions. The liquid pressure is affected by the interfacial shear stress and the surfa ce tension effects as indicated in Eq. (315). The liquid pressure sharply increases in the evap orating region and remains unchanged in the meniscus region as sh own in Figure 320A and the vapor pressure sharply increases at the end of the extensive meni scus region as illustrated in Figure 320B. The sharply increased vapor pressure results fr om the interfacial shear stress given by Eq. (317). Although the change of the vapor pressure is ve ry sm aller than of th e liquid phase, the vapor pressure cannot be ignored due to the effects of the interfacial shear stress. It is shown that the effect of the liquid pressure is dominant in the mass transfer through the liquidvapor interface rather than that of the vapor pressure. The mean velocity variations of the liquid and vapor phase are shown in Figure 321 and Figure 322. The mean velocities of liquid and vapor phases are obtained from Eqs. (345) and (346). While the mean velocity of the liquid phase approach es the inlet velocity, 0.420 m/s, as depicted in Figure 321, the vapor mean veloc ity sharply drops as the liquid film thickness PAGE 81 81 significantly increases along the axial direction as illustrated in Figure 322. The magnitude of mean velocities of the liquid and vapor phases slowly increases along th e axial direction and increases sharply toward the extensive meniscus region. The reason is that the liquid film thickness drastically increases in this region. It leads to the fact that the magnitude of the vapor velocity is much higher than that of the liquid velocity. The interfacial shear stress distributions at di fferent wall temperatures are shown in Figure 323. The results indicate that the interfacial shear stress increases as the wall temperature increases. The interfacial shear stress mainly aff ects the vapor pressure gr adient as indicated in Eq. (317). The sharp change of interfacial shear str ess affects the liquid film shape and results in the increasing vapor pressure. Figure 324 illustrates the interfacial temp erature variations with different wall temperatures. The interfacial temperature is mainly affected by the capillary pressure and disjoining pressure based on Eq. (340). It decreases along the axial direction as the liquid film thickness increases. It is caused by the increasing capillary pressure and the decreasing disjoining pressure. Figure 325 shows the local heat transfer coeffi cient with increasing wall temperature. The local heat transfer coeffi cient can be obtained by ,ln/l czk h RRR (364) As the wall temperature increases, the maximum value of heat transfer coefficient increases. The local heat transfer coefficient is very large value in thin film re gion, while that in the extensive meniscus region is very small. The reason is that the heat transfer coefficient decreases as the thermal resistance increases because the li quid film thickness increases based on Eq. (364). PAGE 82 82 The average heat transfer coefficient profiles with different wall temperatures are shown in Figure 326. It is shown that as the wall temp erature increases, the average heat transfer coefficient increases. The increased wall temperat ure results in the increased evaporative mass flux at the liquidvapor interface. It causes the liquid film thickness and the length of evaporating region to decrease. The heat tran sfer is facilitated by the decreased thermal resistance and liquid film thickness. 3.4.3 Pressure Drop Calculations based on the current model were performed to generate results for comparison with the experimental results of Qu and M udawar (2004) where diffe rent heat flux with 22400.1/ and 401.9/Gkgmskgmswere used in the experiment used by microchannels heat sinks as shown in Figure 327. They used re ctangular microchannels comprising twenty one channels covering a heated area of 4.48 cm2. Each channel is 4.48 cm long with a hydraulic diameter of 348 m. The deionized water enters the channels at 30 oC and 60 oC with the same mass velocity and the outlet pr essure maintains 117 kPa. The thermodynamic equilibrium quality can be obtained from ,wwb eq lchlvqPL x mh (365) where the wall heat flux wHSHSsqqANA bL is the length of twophase region, wP is heated perimeter per channel and ,lchm is the mass flow rate per channel. The length of twophase region can be obtained from ,,lchplsatin bt wwmcTT LL qP (366) Total pressure drop in microchannel can be obtained by PAGE 83 83 totalsptpo p ppp (367) where s p p is single phase pressure drop due to fr ictional force and accelerational force caused by specific volume change tp p is two phase pressure drop due to wall frictional force and accelerational force caused by phase change, and o p is the pressure drop of the outlet region. To compare the pressure drop of the current model and the experimental data of Qu and Mudawar, the single phase pre ssure drop before nucleate boili ng has been obtained by the general pressure drop relation (Carey, 1992). Figure 328 shows the comparison of pressure drops for the respective maximum pressure drop for different heat flux for the same ma ss velocity. Qu and Mudawar (2004) reported maximum pressure drops of 6.62 kPa at 155.5 W/cm2, 12.04 kPa at 202.97 W/cm2, and 17.35 kPa at 240 W/cm2 for 2400.1/Gkgms and Tin = 30 oC. The deviation of pressure drops between the current model and the experimental result at 2142.71/effqWcm and 2240.0/effqWcm is 21.26 % and 33.13 %, respectively, as shown in Figure 328A. The comparison of pressure drops for 2401.9/Gkgms and Tin = 60 oC is shown in Figure 328B. The reported maximum pressure drops for 299.54/effqWcm and 2204.39/effqWcm are 4.10 kPa and 17.38 kPa, respectively. The discre pancy between the cu rrent model and the experimental results for 299.54/effqWcm and 2204.39/effqWcm is 25.28 % and 39.16 %, respectively. This difference may mainly result in the microchannel geometry and the heated region. It may be caused by the fact that the cont raction pressure losses and expansion pressure recoveries in the plenums are ignored in the current model. Although the used geometry of a microchannel is different, the results show that the calculated pressure drop using the current model is approximately close to the experimental results of Qu and Mudawar (2004). The current PAGE 84 84 model considering the surface tens ion effects will provide the good information for the pressure drop in microchannel heat sinks rather than twophase separated model used in the microchannels in which the surface te nsion does not play critical role. PAGE 85 85 linm linT q Figure 31. Physical model for evaporation phe nomena through the liquidvapor interface in a microchannel. 2L 0 2L Figure 32. Geometry and coordinate system used in a numerical simulation. W Lt HwH e a t s i n k HSq zHch Figure 33. Microchannel heat sink unit cell used in a numerical simulation. PAGE 86 86 Set Boundary Conditions Calculate Recalculatethermophys ical properties End Yes No ,,,,,,,lvdvvpppQT Calculate thermophysical properties and initial conditions Calculate the initial film thickness Read ,,,,,,wcoolsatvintTTTTQDm Calculate the first order derivatives, the evaporative mass flux, the blowing velocity, and shear stress R Calculate ,,ilcTTh Calculate ,,,lvlvmmuu constwT Calculate ,,,wilcTTTh Yes No Figure 34. Solution procedure for obtaining unk nown variables. PAGE 87 87 02468101214 0 20 40 60 80 100 120 Liquid film thickness (m)z (m) qw = 0.652 MW/m2 qw = 0.782 MW/m2 qw = 0.869 MW/m2 qw = 1.043 MW/m2 Figure 35. Effects of heat flux on the liquid film thickness. 0.00.20.40.60.81.01.2 5 6 7 8 9 10 Evaporative mass flux x102 (kg/m2s)Heat flux (MW/m2) Figure 36. Evaporative mass flux vari ations at different heat fluxes. PAGE 88 88 02468101214 0 2 4 6 8 qw = 0.652 MW/m2 qw = 0.782 MW/m2 qw = 0.869 MW/m2 qw = 1.043 MW/m2 Liquid pressure gradient x109 (Pa/m)z (m) Figure 37. Liquid pressure gradient variations at different heat fluxes. 02468101214 0 1000 2000 3000 4000 5000 6000 qw = 0.652 MW/m2 qw = 0.782 MW/m2 qw = 0.869 MW/m2 qw = 1.043 MW/m2 pl (Pa)z (m)A Figure 38. Effects of heat fl ux on pressure difference along the ax ial direction. A) Liquid phase. B) Vapor phase. PAGE 89 89 02468101214 0 10 20 30 qw = 0.652 MW/m2 qw = 0.782 MW/m2 qw = 0.869 MW/m2 qw = 1.043 MW/m2 pv (Pa)z (m) B Figure 38. Continued 02468101214 300 250 200 150 100 50 0 qw = 0.652 MW/m2 qw = 0.782 MW/m2 qw = 0.869 MW/m2 qw = 1.043 MW/m2 Liquid mean velocity (m/s)z (m) Figure 39. Mean velocity variations of the liquid phase at different heat fluxes. PAGE 90 90 02468101214 600 500 400 300 200 100 0 qw = 0.652 MW/m2 qw = 0.782 MW/m2 qw = 0.869 MW/m2 qw = 1.043 MW/m2 Vapor mean velocity (m/s)z (m) Figure 310. Mean velocity variations of the vapor phase at different he at fluxes. 02468101214 0 400 800 1200 1600 2000 qw = 0.652 MW/m2 qw = 0.782 MW/m2 qw = 0.869 MW/m2 qw = 1.043 MW/m2 Interfacial shear stress (N/m2)z (m) Figure 311. Interfacial shear stress va riations at differ ent heat fluxes. PAGE 91 91 02468101214 385 390 395 400 405 410 qw = 0.652 MW/m2 qw = 0.782 MW/m2 qw = 0.869 MW/m2 qw = 1.043 MW/m2 Interfacial temperature (K)z (m) Figure 312. Effects of heat flux on interfacial temper ature variations. 80 120 160 200 240 380 390 400 410 420 430 G=400.1kg/m2s Temperature (K)Heat flux (W/cm2) Experimental results (Qu and Mudawar, 2004) THS,i THS,o Tw,i Tw,o Figure 313. Comparison of channel wall and he at sink bottom temperature to experimental results investigated by Qu and Mudawar (2004). PAGE 92 92 01234 0 2 4 6 8 10 Heat transfer coefficient x106 (W/m2K)Z qw = 0.652 MW/m2 qw = 0.782 MW/m2 qw = 0.869 MW/m2 qw = 1.043 MW/m2 Figure 314. Local heat tran sfer coefficient variations for different heat fluxes. 0.6 0.7 0.8 0.9 1.0 1.1 30 40 50 60 70 Average heat transfer coefficient (kW/m2K)Wall heat flux (MW/m2) mass flow rate = 3.823x105 kg/s Tin = 30 oC Figure 315. Average heat transfer coeffici ent variations with various heat fluxes. PAGE 93 93 0.000.040.080.120.160.20 20 40 60 80 100 Current model Kandlikar's correlation Channel exit Tin=30oC Heat transfer coefficient (kW/m2K)Thermodynamic equilibrium quality mass flow rate = 3.823x105 kg/s Figure 316. Heat transfer coefficient variat ions with different th ermodynamic equilibrium quality. 020406080100120 0 20 40 60 80 100 120 Tw = 391.24 K Tw = 391.91 K Tw = 392.32 K Tw = 393.06 K Liquid film thickness (m)z (m) mass flow rate = 3.823x105 kg/s Figure 317. Liquid film thickness varia tions for different wall temperatures. PAGE 94 94 020406080100120 0.00 0.05 0.10 0.15 0.20 mass flow rate = 3.823x105 kg/s Tw = 391.24 K Tw = 391.91 K Tw = 392.32 K Tw = 393.06 K Evaporative mass flux (kg/m2s)z (m) Figure 318. Effects of wall temp erature on evaporative mass flux. 020406080100120 0.0 4.0x1078.0x1071.2x1081.6x1082.0x108 mass flow rate = 3.823x105 kg/s Tw = 391.24 K Tw = 391.91 K Tw = 392.32 K Tw = 393.06 K Liquid pressure gradient (Pa/m)z (m) Figure 319. Liquid pressure gradient variations for different wall temperatures. PAGE 95 95 020406080100120 0 1000 2000 3000 4000 mass flow rate = 3.823x105 kg/s Tw = 391.24 K Tw = 391.91 K Tw = 392.32 K Tw = 393.06 K pl (Pa)z (m)A 020406080100120 0 200 400 600 800 1000 mass flow rate = 3.823x105 kg/s Tw = 391.24 K Tw = 391.91 K Tw = 392.32 K Tw = 393.06 K pv (Pa)z (m)B Figure 320. Effects of wall te mperature on pressure difference along the axial direction. A) Liquid phase. B) Vapor phase. PAGE 96 96 020406080100120 8 6 4 2 0 Tw = 391.24 K Tw = 391.91 K Tw = 392.32 K Tw = 393.06 K Liquid mean velocity (m/s)z (m) mass flow rate = 3.823x105 kg/s Figure 321. Mean velocity profiles of the liquid phase at different wall temperatures. 020406080100120 600 450 300 150 0 mass flow rate = 3.823x105 kg/s Tw = 391.24 K Tw = 391.91 K Tw = 392.32 K Tw = 393.06 K Vapor mean velocity (m/s)z (m) Figure 322. Mean velocity profiles of the vapor phase at different wall temperatures. PAGE 97 97 020406080100120 0 400 800 1200 1600 2000 mass flow rate = 3.823x105 kg/s Tw = 391.24 K Tw = 391.91 K Tw = 392.32 K Tw = 393.06 K Interfacial shear stress (N/m2)z (m) Figure 323. Influences of wall temp erature on interfacial shear stress. 020406080100120 380 382 384 386 388 390 mass flow rate = 3.823x105 kg/s Tw = 391.24 K Tw = 391.91 K Tw = 392.32 K Tw = 393.06 K Interfacial temperature (K)z (m) Figure 324. Interfacial te mperature profiles for diffe rent wall temperatures. PAGE 98 98 020406080100120 0.0 4.0x1048.0x1041.2x1051.6x105 mass flow rate = 3.823x105 kg/s Tw = 391.24 K Tw = 391.91 K Tw = 392.32 K Tw = 393.06 K Heat transfer coefficient (W/m2K)z (m) Figure 325. Effects of wall temperatur e on local heat transfer coefficient. 391.0391.5392.0392.5393.0393.5 0.0 2.0x1044.0x1046.0x1048.0x104 mass flow rate = 3.823x105 kg/s Average heat transfer coefficient (W/m2K)Wall temperature (K) Figure 326. Average heat transfer coefficien t variations for differe nt wall temperatures. PAGE 99 99 ,,,,linlinTpG wq 0eqx eqeqo x x Figure 327. Microchannel heat sinks operating at s ubcooled condition. PAGE 100 100 140160180200220240 0 10 20 30 40 G=400.1kg/m2s, Tin=30oC p (kPa)Heat flux (W/cm2) Exp. results (Qu and Mudawar, 2004) Current modelA 100120140160180200220 0 10 20 30 40 G=401.9kg/m2s, Tin=60oC p (kPa)Heat flux (W/cm2) Exp. results (Qu and Mudawar, 2004) Current modelB Figure 328. Comparison of pressure drops of th e current model with the experimental results observations by Qu and Mudawar (2004) fo r different inlet temperatures. A) Tin = 30 C. B) Tin = 60 C. PAGE 101 101 CHAPTER 4 NUMERICAL SIMULATION OF FORCED CONVECTIVE BOILING IN A MICROCHANNE L Heterogeneous boiling is a phasechange process which produces vapor bubbles on a heated surface in a superheated liquid and results in a highly efficient heat transfer behavior. It occurs in many applications in power generation industries, electronic cooling systems, and heat exchangers, etc. Recently, importance of boili ng heat transfer in a microchannel becomes increasingly stressed due to high heat dissipation rates with th e emergence of microand nanoelectromechanical systems. The recent studies ha ve investigated high heat transfer performance of flow boiling in microchannels (Zhang et al 2000; Jiang et al., 2001; Wu et al., 2006). However, a number of researches focused on the experiments with various geometries and different operational conditions to analyze bubble dynamics such as onset of nucleate boiling, bubble growth and departure. In addition, physical understanding of bubble dynamics under convective flow boiling in a microchannel is st ill limited because of the complexities of twophase flow boiling in a confined space. In this chapter, a numerical simulation me thod of twophase flow boiling in a microchannel has been developed to investigate the physic al phenomena of bubble dynamics and heat transfer. This simulation techni que uses a finitevolume numerical method in conjunction with the Volume of Fluid (VOF) algorithm for tracking the moving and changing liquidvapor interface associated with bubbles in the twophase flow. The obtained results have been compared to the experimental data in the open literature. 4.1 Background and Literature Review Bubble dynamics should be analyz ed to predict and understand the fundamental nature of flow boiling heat transfer in channel flows. Wh en nucleate boiling is first initiated on a heated surface, bubble growth is mainly governed by the momentum exchange between the bubble and PAGE 102 102 liquid. After the initial stage of the bubble growth, the bubble is growing basically due to considerable evaporation through the liquidvapor interface. The latent of vaporization is provided by the heat diffusion from the liquid near the liquidvapor interface. Thus, the majority of the ebullition process is primarily cont rolled by the heat and mass diffusion. Ebullition process of bubble nucleation, growth and departure from a heated surface is shown in Figure 41. Bubble growth on a heated surface is much more complicated due to a thermal boundary layer near the wall. Before bubble grows, the bubble may stay on a cavity mouth until the thermal boundary layer is esta blished. Hsu (1962) noted that the liquid temperature at the bubble tip is higher than the vapor temperat ure inside bubble in the thermal boundary layer. The bubble in the cavity will rapidly grow when a bubble growth criterion is satisfied. Evaporation rate at the liquidva por interface may be nonuniform during the bubble growth. It becomes larger near the heating wall and smaller near the bubble tip. A microlayer beneath a bubble is created by the initially rapid bubble growth in a conventional channel (Moore and Mesler, 1961). The ev aporation within the microlay er has a significant effect on bubble growth from a heated wall. A bubble departs from a heated surface after it gr ows to a certain size. For flow boiling in a conventional channel, the major forces act ing on a bubble include buoyancy, inertia, surface tension, and drag forces. Inerti al and surface tension forces tend to attach the bubble to the wall, while as the bubble becomes larger, the detachment of the bubble from the wall is governed by buoyancy, drag and lift forces. The bubble det aches from the wall when the lifting effects become large enough to overcome the retaining effect s of inertia and surface tension forces at the contact line. A number of researches have an alytically and experime ntally studied bubble dynamics of boiling in conventional channels in the past years, while fewer studies of flow PAGE 103 103 boiling in microchannels are considered in past years (Hapke et al., 2000; Qu and Mudawar, 2002; Wu and Cheng 2003c; Li and Pe terson, 2005; Kuo et al., 2006). 4.1.1 Onset of Nucleate Boiling (ONB) In forced flow boiling, the onset of nucleat e boiling represents the initiation of nucleate boiling. The first occurrence of bubbles distinguish es the transition from a singlephase to a twophase flow regime. The transition causes the signif icant change of heat transfer and pressure drop. Thus, prediction of the ONB is of important factors to understand flow boiling phenomena. A number of studies for the ONB model in conventional chan nels were conducted in the past decades. Hsu (1962) first proposed a mini mum superheat criterion for the ONB in pool boiling as follows: 212.8llvvwsat ONB satkhTT q T (41) He assumed that the bubble nuc leus would grow only if the minimum temperature surrounding the bubble is equal to the saturation temperatur e of the vapor inside the bubble. Sato and Matsumura (1963) proposed the heat flux for the ONB in terms of wall superheat: 28llvvwsat ONB satkhTT q T (42) Bergles and Rohsenow (1964) extended Hsus model and proposed a graphical solution to predict the incipient heat flux in flow boiling: 0.02342.16 1.15610821.8p ONBwsatqpTT (43) Hino and Ueda (1985) proposed the incipient he at flux of the ONB in terms of a cavity size: 2 ,max,max2ll s a t l v ONB satONB cl v ckk T v qT rh r (44) PAGE 104 104 The cavity size can be obtained from (Collier and Thome, 1994): 1/2 ,2satl ccrit lvONBTk r hq (45) Davis and Anderson (1966) provided an analytic al model of the approach of Bergles and Rohsenow (1964), and introduced the contact angle as a variab le in the prediction of ONB. 281cosllvvwsat ONB s atkhTT q T (46) Kandlikar et al. (1997) numeri cally computed the temperature at the location of the stagnation point around the bubble, which was used as the minimum temperature in the ONB criterion: 28.8llvvwsat ONB satkhTT q T (47) Celata et al. (1997) investig ated the onset of subcooled boi ling in the forced convective flow of water and recommended Thoms co rrelation (1965) for its good match with the experimental data: 20.00195exp0.023ONB wsatqTTp (48) Basu et al (2002) postulated the dependence of the available cavity size on the contact angle and proposed a correlation for the incipient heat flux: ONBspwsatspsatlqhTThTT (49) The onset of boiling in minichannels and microchannels has been studied by several investigators in recent years. Ghiaasiaan and Chedester (2002) proposed a semiempirical method to predict the ONB in turbulent flow in microtubes. Qu and Mudawar (2002) measured the incipient boiling heat flux in a microchannel heat sink and developed a mechanistic model to PAGE 105 105 incorporate both mechanical and thermal c onsiderations. Li and Cheng (2004) employed nucleation kinetics to derive the wall superhea t at the ONB, and included a consideration of the effects of contact angle, dissolved gas, and th e existence of microcavit ies and corners in the microchannels on ONB. Lee et al. (2004) investig ated bubble dynamics in a single microchannel. They used the ONB model proposed by Hino and Ueda (1985) to estimate the incipient heat flux. Li et al. (2004) experimentally studied the bubbly dynamics in two parallel trapezoidal microchannels. The ONB of nucleate boiling in a microchannel heat sink was experimentally studied by Liu et al. (2005). They also developed an analytical model to predict the incipient heat flux and the bubble size at the ONB. 4.1.2 Bubble Departure Diameter Bubble departure diameter on a nucleation site has been pred icted by several correlations. Fritz (1935) suggested the depart ure diameter correlated with the buoyancy and surface tension forces as follows: 0.0208d lvD g (410) where is the contact angle. Cole and Rohsenow (1969) proposed the bubble departure diameter with a modified Jakob number: d lvEo D g (411) with 2 *54EoCJa (412) *lplsat vlvcT Ja h (413) PAGE 106 106 where the coefficient C is 41.510 for water and 44.6510 for other liquids. Kocamustafaogullari (1983) pr oposed the bubble diameter at departure including highpressure data: 0.9 52.6410lv d lvvD g (414) Gorenflo et al. (1986) suggest ed the correlation of the depa rture diameter of bubble with thermal diffusivity of liquid and th e Jakob number at high heat fluxes as 43 42 12 11 3l dJa DC gJ a (415) where 1C has different values for differe nt liquids and the Jakob number Ja is defined by lplwsat vlvcTT Ja h (416) 4.1.3 Bubble Departure Frequency The complete process of liqui d heating, nucleatio n, bubble growth, and departure referred to as the bubble ebullition process is the main mechanism of heat transf er from a superheated wall during nucleate boiling. The bubbl e growth rates tend to be rela ted to the departure diameter and frequency of release. The inverse of the freque ncy of release equals to the sum of the waiting time, wt, and growth time, g t as shown Figure 42. 1wgtt f (417) The waiting time and growth time is difficult to determine quantitatively. The waiting time depends on the temperature fields on a heated su rface and in the liquid ne ar the nucleation site; the growth time, however, depends on the evap oration rate at the bubbl e base, around the bubble, and the bubble diameter at departure. PAGE 107 107 The frequency of bubble release in a conventio nal channel correlates with the departure diameter as follows: .d f Dconst (418) The predictive relations of the bubble frequency of release on a h eated wall have been developed by several researchers (Jakob and Fritz, 1931; Peebles and Barber, 1953; Zuber, 1963). Peebles and Barber (1953) suggested the relation as follows : 14 21.18g lv d wglt g fD tt (419) Zuber (1963) proposed th e following relation: 14 20.59lv d lg fD (420) Malenkov (1971) suggested that 1 1 1d d dvlvU fD Uhq (421) where 2 2lvd d lvdlvgD U D (422) Kennedy et al. (2000) investigated the onset of nucleation boiling and onset of flow instability in minicha nnels. Ghiaasiaan and Chedester (2002) examined the applicability of existing macroscale models on the onset of nucle ation boiling in minichannels and proposed a semiempirical method. Lee et al. (2004) studi ed experimentally bubble dynamics in a single trapezoid microchannel with a hydraulic diameter of 41.3 m. They pointed out that there are still not correlations of the prediction of bubble frequency in microchannels. Dhir (1998) PAGE 108 108 conducted a review on various aspects of boiling h eat transfer including bubble dynamics such as bubble growth, bubble departure, and bubble depart ure frequency. Kandlikar (2002) reviewed the literature on bubble grow th in minichannels. Bubble behaviors in a microchannel are quite different from those in a macrochannel during flow boiling heat transfer Bubble growth in a microchanne l is mainly confined by the channel wall in the streamwise direction and experi ences a very large pressure gradient in that direction. Thus, the bubble size is limited by the channel size and geometry. The buoyancy force may be negligible during the bubble departure process. In contrast, the drag force should be very significant due to a large pressure drop through the channel. The main objective of this study is to invest igate effects of channe l size and heat flux on nucleate boiling. Bubbly dynamics of bubble gr owth, bubble departure size, and departure frequency are studied using computational simula tion. Periodic flow boiling in a microchannel is discussed. The results of this study provide fundamental understanding of bubble dynamics and heat transfer during flow boiling in a microchannel. 4.2 Forced Convection Flow Boiling Numerical Simulation Method A finite volume method (FVM) has been adopted to discretize the gove rning equations for simulating bubble dynamics and heat transfer of forced convective flow boiling in a microchannel. Heat and mass transfer due to phase change between th e liquid and vapor phases have been identified and the volume of fluid (VOF) algorithm is used to capture the moving interface between the liquid and vapor phases. 4.2.1 Governing Equations The VOF method uses the NavierStokes equation s to simulate the twophase flow boiling. The gravitational effects are negligible because the surface effects are dominant compared to them in a microchannel. The contin uity and momentum equations are PAGE 109 109 0 V t (423) M V VVpVS t (424) The energy equation is pET cVTkTS t (425) The source terms SM and SE which denote interfacial mass, momentum, and energy interfacial exchange are included to simulate twophase flow boiling in a microchannel. The source terms will be discussed in the volume of fluid method. 4.2.2 Numerical Method 4.2.2.1 Finite volume method (FVM) The finite volume method is reliable and robust by using the conservation laws and treating the discontinuity at the interface. The finite volume method uses the volume integration for the conservation equations of mass and momentum, and energy (Eqs. (423) (425)) (Patankar, 1980): 0csVndS (426) M CV CS CS CS CVVdVVVndSpndSVndSSdV t (427) pp E CV CS CS CVcTdVcTVndSkTndSSdV t (428) where CV and CS represent the control volume and surface of the control volume, respectively, and n is the outward normal vector from the cont rol volume. Using the finite volume method, PAGE 110 110 the governing integral equations are discretized into the di scretized cell to solve the computational domain numerically as shown in Figure 43: 0faceVA (429) M facesfacesfacesfacesVVVVApAVASV t (430) pP E faces faces facescTVcTVAkTASV t (431) The discretized equations satisfy the conserva tion of appropriate properties for each cell. The scalar values such as pressure and temperature are stored at the cell center (C1, C2, and so on) and the velocity fields at the cell faces are us ed because the influence of pressure field is not properly represented in the moment um equation when the velocity fields are defined at the cell centered. The nodal values are evaluated by the weighted aver age values in th e surrounding cells. For example, the value at node 8 is determined fr om the weighted average of the values in the surrounding cells (C2, C3, C6, and C7). A firstorder implicit scheme is used for advancing the solutions of the integral gove rning equations in time. 1 1 n nn M faces faces facest VVVVApAVASV V (432) 1 1 n nn pp P E faces facest cTcT cTVAkTASV V (433) where n denotes the previous time step and n+1 current time step. The convection terms in Eq. (432) are discretized using QUICK (Quadratic Upstream Interpo lation for Convective Kinetics) scheme with thirdorder accuracy due to numerical diffusion errors of loworder upwind schemes. Th e diffusion terms are discretized with central PAGE 111 111 differencing formula having secondorder accurate. The fully implicit scheme with firstorder in time is applied to the time integration since small time steps are used to ensure numerical accuracy during the simulation. In twophase flow boiling, the ratio of liquid to vapor is very large. Thus, pressure interpolation scheme uses the PRESTO (Pressure Staggering Option) dealing with the discontinuity of pressure gradients for flows w ith rapidly changing densities at the liquidvapor interface. Pressu revelocity coupling calculation uses PISO (PressureImplicit with Splitting of Operators) based on SIMPLE (SemiImplicit Method for PressureLinked Equations) method (Patankar, 1980). PISO involves one predictor step and tw o corrector steps to improve the efficiency of the calculation b ecause the limitation of the SIMPLE and SIMPLEC (SIMPLEConsistent) method is that new veloc ities and corresponding fluxes do not satisfy the momentum balance after the pressure correction equations is solved. The predictor step of PISO solves the moment um equations with the intermediate pressure field p to give velocity components *uand *v using the SIMPLE algorithm. When the correction of the pressure field p is used, the velocity fields *uand *v satisfies continuity. At the first corrector step, new corrected values **u and **v are calculated using the pressure correction p as follows: *** p pp (434) ***uuu (435) ***vvv (436) The second corrector step is solved to obtain the second pressu re correction field p and a twicecorrected pressure field and ve locity fields are obtained from: ***** p ppppp (437) PAGE 112 112 *****uuu (438) *****vvv (439) The calculation will be repeated until the balance is satisfied. The energy equation will be solved using the corrected pressure and velocity fields. Iterative timeadvancement scheme is used for the timeadvancement scheme. All alge braic equations obtained from Eqs. (432) and (433) are solved by the linebyline m ethod using the TDMA (TriDiagonal Matrix Algorithm). For a given timestep, the discretized equa tions are iteratively solved as follows: Define intermediate iV and iT which are initialized with the current nV and nT. Iteratively solve the implicit equations un til the intermediate solutions satisfy the equations: 1,1,2,3, (i = iteration number)in iit VVFVVi V (440) 1 1,1,2,3, (i = iteration number)nn ii ppt cTcTFTTi V (441) The intermediate iV and iT come to be the solution for new time step n+1 as the intermediate solutions converge. The segregated solution method linearizes the discrete, nonlinear governing equations to make a system of equations for the dependent va riables in every computational cell and it has a robustness to solve twophase flow boiling. The resu ltant linear system is solved sequentially as follows: Based on the current solution and/or the initial ized solution, fluid properties are updated. The momentum equations are solved using current values for pressure and face mass fluxed before updating the velocity field. A pressure correction equation is then solved to obtain the corrected pressure and velocity fields and the face mass, which satisfy the continuity equation; since the velocities obtained from the previous step may not sati sfy the continuity equation, the pressure correction equation is derived fr om the continuity and the linearized momentum equations. PAGE 113 113 The energy equation is solved using the previ ously updated values of the other variables. These steps are continued until the convergence criteria are me t at a current time step. The momentum, energy equations at a new tim e step are calculate d if the convergence criteria are satisfied. 4.2.2.2 Volume of Fluid (VOF) A volume of fluid model was chosen to analyze twophase flow boiling with the motion of large bubbles in a liquid phase. The VOF technique solves a single set of momentum equations and determines the interface shap e of the two immiscible fluids by tracking the volume fraction of each phase in each computational cell. Th e sum of the volume fraction of both phases is unity in each control volume. In this model, the vapor phase will be tracked and the volume fraction of the vapor phase is v Each computational cell is within three conditions as follows: 0v : liquid phase only 1v : vapor phase only 01v : twophase mixtures in the cell The interface shape between both phases is track ed by solving the continuity equation for the volume fraction of phases. The continui ty equation can be written as follows: vv vvVS t (442) where v is the volume fraction for the vapor phase and S is the source term which represents the mass transfer accounting for phase change between the liquid and vapor phases. The volume fraction equation cannot be solved for the primary phase. In the current model, the liquid phase is set to be the primary one; thus, the volume frac tion of the primary phase can be obtained from 1lv (443) PAGE 114 114 The discretized volume fraction equation with th e firstorder implicit scheme is obtained by using the finite volume approach: 1 1 n nn vvvv vv facest VASV V (444) Based on the void fraction, a single momentum equa tion is solved for the entire computational domain and the resulting velocity field is shared among the phases as follows: s vVVVpVgF t (445) where p is the pressure, 0, gg is the gravitational acceleration, s vF is the surface tension force per unit volume and and are the density and viscosity. Equation (445) depends on the volum e fractions of all phases through the properties and These properties are obtained by the presence of the computational phases in each control volume. The density and viscosity at the liquidvapor interface vary with the volume fraction tracked as follows: 1vlvv (446) 1vlvv (447) The source term in Eq. (445) is the surface tension force per unit volum e and the gravitational forces are negligible in current study. The c ontinuum surface force model of Brackbill et al. (1992) is used in Eq. (445). Figure 44 shows a typical twophase ce ll with an immersed phase interface. The surface curvature is defined in terms of the divergence of the unit normal vector Kn (448) where K is the surface curvature and n is the unit normal vector of the interface given by n n n (449) PAGE 115 115 where n is the surface normal defi ned as the gradient of v as vn (450) Combining Eqs. (449) and (450) with Eq. (448), the curvature of th e interface can be rewritten as v vK (451) The source term in Eq. (445) can be expressed as a volume force as 12v sv lvK F (452) where is obtained from Eq. (446). The contact angle w is used to adjust the surface normal in cells near the wall. It results in the adjustment of the curvature of the su rface near the wall. Thus, the surface normal at the cell near the wall is defined as cossinwwwwnnt (453) where wn and wtare the unit normal and tangential vectors to the wall, respectively, and w defines the angle between the wall and the tangent to the interface at the wall. The contact angle w measures inside the first phase as shown in Figure 45. Eq. (453) with the normally calcu lated surface normal one cell away from the wall determines the local curvature of the surface. To capture the liquidvapor interface during the simulation, convection and diffusion fluxes through the control volume faces are calcula ted and balanced with source terms within the control volume in the FVM approach. Thus, in th e VOF method, the interface is reconstructed at the beginning of a time step for more accurate simu lation of its evolution as shown in Figure 46. When the cell is near the inte rface between two phases, the geom etric reconstruction scheme is PAGE 116 116 used. This scheme represents the interface between two phases using a piecewiselinear approach and uses this linear shape for calculation of th e advection of fluid th rough the cell faces. The procedures of the scheme are The orientation of the curve within each twophase cell is determined by Eq. (449) and (450). Given the orientation of the planar that represents the interface in a cell, determine the position nd of the linear interface relative to the center of each twophase cell based on the volume fraction of the vapor as shown in Figur e 44. These steps are often referred to as the interface reconstruction st eps (Rider and Kothe, 1998). The mass fluxes through each face are computed using the linear interface representation and information about the normal and tangentia l velocity distribution on the face as shown in Figure 47. Once the mass fluxes are calcul ated in one direction, the interface is reconstructed before the fluxes are de termined in the other direction. The volume fraction in each cell is determin ed by using the balance of fluxes calculated during the previous step. The energy equation for the mixture is given by E E VEpkTS t (454) where ES is the energy source term, and describes th e phasechange due to the latent heat of vaporization through the li quidvapor interface and k is the thermal conductivity defined by 1vlvvkkk (455) The VOF model treats energy and temper ature as massaveraged variables: 1 1 n qqq q n qq q E E (456) where q E for each phase is based on the specific heat of that phase and the shared temperature. The liquid and vapor phase energies are PAGE 117 117 and llvvEhEh (457) where h is the sensible enthalpy of theqth phase defined by ,refT qpq ThcdT (458) The reference temperature used to calculate the enthalpy of each fluid is 298.15 K in Eq. (458). The discretized form s of momentum and energy equation are given by Eqs. (432) and (433). The boundary conditions are 0 and specified at 0 and ppuv xLyL (459) 0 and 0 at and tptpuv LLxLyL xy (460) 0 at and g age tp t p p LLxLyL (461) at 0 and in p pTTxLyL (462) at and y=specifiedeffp tpT kqLxLL y (463) where tL is the total length of a microchannel and pL is the plenum length of inlet and outlet. The boundary conditions at the wall as 0 at the wall uv (464) 0 and =0 expect the bottom wall TT xy (465) In current study, FLUENT software having a ro bustness and stability is used to solve the continuity, volume fraction equation, and momentum and energy equations and to analyze the forced convective boiling in a microchannel. PAGE 118 118 4.3 Problem Formation and Description Since extensive computational resources are needed to resolve the forced convection boiling systems, both twodimensional and threedimensional simulations are proposed. As the twodimensional simulation is much less demanding in computing times than the threedimensional case, so strategically, the former is used for general physical trends and parametric evaluations while the latter is adopted for more specific cases such as th e flow boiling instability in a microchannel in Chap. 5. 4.3.1 TwoDimensional Model The schematic of the twodimensional copper microchannel to be i nvestigated in this chapter is shown in Figure 48. The microchannel has a length of 60 mm. The channel dimensions are listed in Table 41. Fi ve different channel heights are used: Hch = 1 mm, 700 m, 500 m, 300 m, and 30 m. The length to height aspect ratio ranges from 160 to 1,111. Therefore, the channel length is long enough to provide both entr y and fully developed phenomena. Three constant heat fluxes of 200 kW/m2, 387 kW/m2, and 500 kW/m2 are used and the constant heat flux condition is placed only at the outside surface of the bottom channel wall which has a thickness of 3 mm thickness. Initially subcooled water as th e working fluid with a given inlet temperature of Tin = 30 C and various velocities in the range of 0.015 m/s to 0.5 m/s enter from the inlet plenum to the channel inlet. A mixture of liquidvapor flow leaves the channel as heat is being added to the worki ng fluid. The thermophysical properties of the working fluid used in the simula tion are listed in Table 42. 4.3.2 Numerical Infrastructure The finite volume algorithm and the VOF method presented and discussed above for obtaining the simulation solutions are implemented using the Fluent software package. The computational model was created and meshed using the GAMBIT software available in the PAGE 119 119 Fluent package. The grid system and the boundary condition used to simula te the twophase flow boiling are illustrated in Figure 49A. The grid density of the cha nnel using the fine mesh for Hch = 30 m is shown in Figure 49B. The model has the similar dimensions as some of the actual channels in the literature excep t that the plenum and channel length was modeled. The grid spacing for each case is in the range of 3 m to 50 m to be sufficiently smaller than the bubble departure diameter. A grid density for each case is given by Table 43 and larger meshes are only used in the portion of in let and outlet plenums. The boundary conditions used in the nu merical simulation are as follows: At the liquid inlet, the fluid velocities a nd subcooled temperature are specified and the volume fraction of the vapor phase v is set to be zero. At the outlet, the pressure as gauge pressure is set to zero: FLUENT uses internally the gauge pressure. Absolute pressure is the sum of the operating pressure and the gauge pressure in FLUENT. In the model, the defa ult value of 1 atm (101.325kPa) is assumed as the operating pressure and the absolute pressure at the outlet is assumed as 1 atm. No slip condition is assumed on the wall. Surface tension effects are considered and the contact angle w is specified. Constant heat flux is imposed on the heat sink bottom wall only. The other channel walls are set to an adia batic condition except the bottom wall in a microchannel. Mass transfer and phasechange effects th rough the liquidvapor interface are included during the simulation. The source terms including mass transfer be tween the liquid and vapor phases and vaporization through the liqui dvapor interface in Eqs. (442) and (454) are set up by User Defined Functions (UDFs). The m ass transfer terms for each phase are equal and have opposite signs: vlSS (466) PAGE 120 120 vEl vSSh (467) where lS andvS denote the mass transfer for each phas e through the interf ace, respectively. The mass transfer at the onset of nucleate boili ng was obtained from the incipient heat flux proposed by Sato and Matsumura (1963) 2 ,1 8llvsatONB v satkT S T (468) The wall superheat s atONBT for the onset of boiling and the superheated liquids, s atT for evaporation through the interface (Ka ndlikar, 2004) was calculated by ,4 11 2satlvsp llvsub satONB llv satlvspTvh khT T kh Tvh (469) and s ubsatinTTT is the subcooling. Eq. (469) can be written in terms of the hydraulic diam eter hD ,4 11 2satlv lvhsub satONB lvh satlvTvChDT T hD TvC (470) with sph lhD NuC k (471) where s ph is the single phase heat transfer coefficien ts. In this study, the fully developed Nusselt number for a threeside heated rectangular channel wa s used to obtain the wall superheat at the onset of boiling (Kandlikar, 2006). Phase change through the interface in superheat liquids was considered by the empirical model proposed Lee (1979) 1l s at vl s atT SC T (472) PAGE 121 121 where C is the relaxation parameter which controls the convergent speed a nd solution stability. The solution using the higher relaxation parame ter becomes unstable, while using the lower value may slowly be convergent. The relaxation pa rameter of 0.6 was used by trial and error to ensure the stable and fa st convergent solution. The local curvature is determined from the fi rst phase as shown in Figure 45. The first phase is set to be the vapor phase The contact angle approaching 175o was specified for highlywetting liquids; if the defined contact angle is in correct, the bubble will not detach from the wall and instead spread along it. Th e bubble with the small contact angle slowly moves toward the downstream under the inertial forces. The computational procedure including the implementation of the UDFs is shown in Figure 410. The UDFs accounting for vaporization are applied to the onset of boiling and twophase regions. These UDFs were written in C programming language and implemented in the FLUENT software. A segregated unsteady solver to describe temporal vari ations of bubbles was used in the simulation. The PRESTO (pressure staggering option) scheme was used for the pressure interpolation, the QUICK (Quadratic In terpolation for Convective Kinetics) scheme was adopted for the momentum and energy equations. The PISO (pressureimplicit with splitting of operators) scheme was chosen for the pressurevelocity coupling, th e Courant number which defines the ratio of the time step to the cell resistance time wa s set to be 0.25 for the volume fraction calculation. The QUICK formulation is th irdorder accurate which helps to mitigate the unfavorable effect of artificial diffusion that can occur when using loworder upwind schemes. The simultaneous computations of the energy an d momentum equations require a large number of iterations and for some cases one and half month was requir ed to obtain the computational results. PAGE 122 122 4.4 Results and Discussion 4.4.1 Onset of Nucleate Boiling Figure 411 shows the wall temper ature profiles when nucleat e boiling first appears near the upstream region for the cha nnel heights of 1 mm and 700 m with q = 500 kW/m2 and m = 3.4486 kg/s. Nucleation, for channel heights of 1 mm and 700 m, occurs at 67 ms and 58.5 ms, respectively. Nucleate boiling first initiate s when the wall temperature exceeds the onset of nucleate boiling superheat value. The wall temperature drops sharply near the nucleation position and the wall temperature distribution becomes approximately cons tant in the downstream region of nucleation site as shown in Figure 411A a nd B. The reason is that when the bubble is generated on the heated wall, the wall temperatur e drops due to the large amount of latent heat consumption for facilitating vaporization. However, in experiments, the position of nucleation is not easily determined because bubble embryos occu rred on the heat surface are extremely small (Hino and Ueda, 1985). For this reason, the onset of nucleate boiling position has been assumed to correspond to the position where the wall te mperature abruptly decreases by experimental researchers. However for numerical simulations, the boiling nucleation is readily detected as shown in Figure 415 for the case of Figure 411A. The results as illustrated in Figures 411A and B indicate that the wall temperature for Hch = 1 mm is higher than that for Hch = 700 m at the given mass flow rate, m = 3.4486 kg/s. The velocities at the microchannel inlet for Hch = 1 mm and 700 m are 0.015/ms and 0.0214/ms, respectively. Thus, for the case of Hch = 1 mm, the Reynolds number is 20 as compared to 26 for the case of Hch = 700 m. The lower Reynolds number results in a higher superheat requirement or higher wall temperature for nucleation for the case of Hch = 1 mm as evidenced in Figure 411. Figure 411A shows two bubbles initiated on the heated wall for Hch = 1 mm, while for Hch = 700 m, only one bubble PAGE 123 123 generated on the heated wall as indicated in Figure 411B. The reason is that the higher convection rate for the case of Hch = 700 m results in a larger heat removal rate that does not leave enough heat for supporting a second bubble nucleation. Additionally, since the second bubble blocks some liquid flow, the wall temperat ure at the upstream of the bubble immediately increases as shown in Fi gure 411A for all two b ubble formation sites. Figure 412 illustrates the wall temperature pr ofiles when the firs t bubble occurs in a microchannel of Hch = 500 m at q = 500 kW/m2 for two different mass flow rates. The nucleate boiling initiations for m = 3.4486 kg/s and m = 1.7245 kg/s take place at times of 53.8 ms and 52 ms, respectivel y in the downstream region due to higher Reynolds numbers of 33 and 162 at the channel inlet as compared to both of Hch = 1 mm and 700 m cases. As shown in Figure 412, similar to previous cases, the wall temperature drastica lly decreases when the bubble generates on the heated surf ace due to the vaporization of superheated liquids near the wall and the transient conduction from the wall to the bulk liquid. The onset of nucleate boiling location slightly moves toward the channel exit with increasing mass flow rates. The results show that the wall temperature during nucleation for Hch = 500 m is lower than those of Hch = 1 mm and 700 m during the bubble nucleation on the heated wall. It is caused by the higher convective effects to remove heat more efficiently from the surface than Hch = 1 mm and 700 m cases. The higher singlephase heat transfer coefficients for cases of Hch = 500 m with Re33hD and Re162hD significantly affect the wall temper ature levels as compared to the results of Hch = 1 mm and 700 m with Re20hD and Re26hD respectively. The wall temperature variations during inci pient boiling for differe nt heat fluxes of q = 387 kW/m2 and 500 kW/m2 for Hch = 300 m with the same mass flow rate of m = 3.4486 PAGE 124 124 kg/s are shown in Figure 413. The first bubble is observed further downstream region similar to the case for Hch = 500 m due to relatively high er Reynolds number, Re45hD resulted from the small channel height. For different heat fl uxes with the same mass flow rate and same channel height, the nucleation time for q = 500 kW/m2 is much shorter than that for q = 387 kW/m2 because more heat conducts from the mi crochannel bottom wall to the channel wall. However, the nucleation position rema ins the same for both cases. For the smallest channel of Hch = 30 m at q = 500 kW/m2 and m = 3.4486 kg/s, the nucleation phenomenon is more wide spread spontaneously as many sites become active simultaneously as illustrated in Figure 414. Seve ral smallest bubbles are nucleated almost at the same time on the entire microchannel heated su rface after a very short simulation time t = 16ms as oppose to one or two nucleation sites in the larger channel cas es (see Figures 411 413). The heater surface temperature distribution that shows a lot of small fluctuations (.2 oC) due to numerous bubble formations is plotted in Figure 414. The bubble distributi on pattern is given in Figure 421A. The reason is that the bulk liquids are rapidly s uperheated relatively uniformly due to very high heat transfer ra tes in a very small channel. In addition, the wall temperature at nucleation for Hch = 30 m is lower than for Hch = 1 mm, 700 m, 500 m, and 300 m. It is mainly due to the extremely small channel hei ght that results in lower degree of superheat requirement for onset of nucleate boiling as given in Eq. (470). 4.4.2 Effects of Channel Height The range of channel height in the current study goes from 30 m to 1 mm and it covers both microscale and miniscale boiling phenomena. The effects of channel height are summarized in Table 4, where we note that the heat flux and mass flow rate for all the five cases are kept the same. As a result, the Reynolds num ber increases with decreasing channel height. PAGE 125 125 The ONB time decreases steadily with decreasing ch annel height while the first nucleation axial location increases with decreasing channel height with the only exception that the nucleation covers almost the entire channel length as explained before. 4.4.3 Bubble Coalescence and Departure Information for bubble growth, coalescence and departure on the superheated wall is important for the understanding of microchanne l boiling. With the same conditions of q = 500 kW/m2 and m = 1.7245 kg/s, three channel sizes of Hch = 700 m, 500 m, and 300 m were investigated. The results are give n in Figures 416, 417 and 418 for Hch = 700 m, 500 m, and 300 m, respectively. Even though for three different channel sizes, the bubble growth coalescence and departure are relatively similar. For Hch = 700 m case, there are two bubbles generated at to + 0 ms (see Figure 416A) and then they c ontinue to grow respectively until to +31.1 ms when the coalescence takes place (as shown in Figure 416 B, C and D). The coalescenced bubble remains on the heater surface and continue to grow (as ill ustrated in Figure 416 E, F and H). Finally it reaches to a size large enough for the buoyancy and flow inertia force to overcome the surface tension force for the bubble to liftoff and depart from the heater surface (as shown in Figure 416I and J). As can be seen that a new bubble is nucleated at the same spot from where the old bubble departed (as shown in Figure 416K and L) As mentioned earlier, all three size cases involve similar mechanisms of bubble growth, merging, departure a nd nucleation from the departure spot. The bubble departure size versus the channel height at q = 500 kW/m2 and m = 1.7245 kg/s is provided in Figure 419. For Hch = 700 m, the departure size is 105 m. It is shown that the bubble departure si ze increases with the channel he ight. The main reason is that PAGE 126 126 since the Reynolds number decreases with the channel height and the depa rture size decreases with increasing Reynolds number as the flow inertia force that works for the departure is directly proportional to the Reynolds number. As a result, bubbles depart earlier with a smaller size when the Reynolds number is larger. 4.4.4 TwoPhase Flow Patterns Development of the twophase flow pattern w ith time is given. Figure 420 illustrates the temporal flow pattern evolution for Hch = 1 mm at q = 500 kW/m2 and m = 3.4486 kg/s. At t = 500 ms, the flow is basically filled with is olated small bubbles as shown in Figure 420A. As time increases, more bubbles are generated that results in bubblebubb le interactions and coalescence (as illustrate d in Figure 420B and C). In gene ral, the results show bubbly flow, slug/plug flow, and the elongated slug/semiannular flow at the simulation time t = 670 ms (see Figure 420D). Stratified flow is not observed during this simulation due to the dominant surface tension forces. The obtained flow patterns show some similarity to those of Kew and Cornwell (1997). Their results present thr ee flow patterns; isolated bubbles, confined bubbles and an annular slug flow regime. Their confined bubble is equivalent to the elongated bubble or a slug flow. Figure 421 shows the twophase fl ow structure for the case of Hch = 30 m at q = 500 kW/m2 and m = 3.4486 kg/s. At t = 16 ms, an array of bubbles are nucleated almost simultaneously due to the extremely small channel height that results in very high heat source density for the coolant as shown in Figure 421A At t = 21 ms, the bubbles have grown to the size of the channel height that establishes a pe riodic boiling condition in the microchannel (see Figure 421B). In summary, the bubbly flow was observed from t = 21 to 22.4 ms (Figure 421B and C). The elongated slug flow occurred at t = 22.9 ms as depicted in Figure 421D. The bubbly PAGE 127 127 flow was detected again at t = 24 ms as illu strated in Figure 421E. Churn flow due to the entrapped liquids between the wall and bubble base was observed at t = 26 ms in Figure 421F. The elongated slug flow occurred again at t = 31 ms as illustrated in Figure 421G. The numerical results in Figure 21 show qualitatively the same tr ends with those found in the experimental data provided by Wu and Cheng (2003c ). They observed the periodic flow boiling: bubbly flow, churn flow, singlephase liquid flow, elongated slug flow, and churn flow. 4.4.5 Heat Transfer and Flow Associated with Elongated Bubbles When the bubbles grow to the size of the cha nnel height, they would continue to elongate in the flow direction. These elongated bubbles may reach local dryout condition. For the two elongated bubbles seen in Figure 421D, the upstream bubble as shown in Figure 422A has reached a local dryout condition where a spike in the wall temperature is associated with the dryout point located at axial position 74.71 mm. The corresponding velocity of the liquid phase as shown in Figure 422B displays a dip at the dryout point. On the other hand, the downstream bubble shown in Figure 423A has not reached any local dryout point as the wall temperature profile does not have any spikes. The relatively hi gher temperature near th e tail of the bubble is caused by the slowdown of the liquid velocity in the region as plotted in Figure 423B. The velocity vector contour plot is shown in Figure 424B for the channel portion that contains the two elongated bubbles from Figure 424A. Basically, the vapor velocity profiles inside the two bubbles and the liquid velocity pr ofile between the two bub bles are relatively parabolic in nature. The velocity magnitudes ar e very small between th e bubble and the channel wall. 4.4.6 Heat Transfer The average heat transfer coefficient as a f unction of the twophase quality for different channel heights at q = 500 kW/m2 and m = 3.4486 kg/s is provided in Figure 425. It is PAGE 128 128 shown that the heat transfer coefficient in the twophase region is much greater than in the singlephase region. It is simila r to results reported in the expe rimental data of Qu and Mudawar (2004). The heat transfer coefficien t sharply increases at lower qualities (~0.11). It is due to the subcooling effects leading to e nhancing the heat transfer. In the saturated boiling region with increasing quality, the heat tran sfer coefficient decreases. Th e transition with nucleate boiling beneath elongated slug flow, local dryout and rewetting process may affect the decrease of the heat transfer coefficients. It is consistent with the experimental observations by several researchers (Ravigururajan, 1998; Kandlikar an d Balasubramanian, 2004; Lee and Mudawar, 2005b; Diaz and Schmidt, 2007). The h eat transfer coefficient is st rongly affected by the channel heights and vapor quality with gi ven mass flow rate and heat flux around the quality of 0.4. It is caused by the constraints of the channel heights. The diameter of vapor bubbles approaches the channel height due to the cons traints and it leads to nucleate boiling, bubbly flow, elongated slug flow as shown in Figure 421. The heat transfer mechanism of nucleate boiling at lower quality is mainly dominant in a microchannel. On the other hand, the effects of vapor quality on the heat transfer are nearly small at highe r qualities. It is caused by the fact that with increasing quality, elongated slug flow becomes dominant pattern. The heat tran sfer is mainly affected by evaporation through the liquid film entrapped by the elongated slug bubble, as already reported by Jacobi and Thome (2002) and Thome et al. (2004). The effects of heat flux on the he at transfer coefficient at Hch = 300 m and m = 3.4486 kg/s are shown in Figure 426. It is shown that the heat transfer coefficients significantly depend on the heat flux and decrease with incr easing quality around 0.11. This behavior is consistent with the experimental investigations by Lee and Mudawar (2005b) and Diaz and Schmidt (2007). As the above mentioned, the sharp decrease of the heat transfer PAGE 129 129 coefficient may be caused by the nucleate boi ling beneath elongated slug bubble, local dryout and rewetting process. It is noted that the heat transfer coeffici ents are a strong function of heat flux; however, the effects of th e quality are very small. PAGE 130 130 s atT Figure 41. Ebullition mechanis m of nucleation and bubble growth and depa rture at an active cavity site. Departure diametertimewt g t g twtt waiting time growth time Figure 42. Bubble release frequency w ith the waiting time and the growth time. Figure 43. Control volume used to discretize the governing equations. PAGE 131 131 n nd Figure 44. Twophase cell with an immersed phase interface. w Figure 45. Measurement of the contact angle near a heated wall (FLUENT, 2005). A 1.0 1.01.0 1.0 1.0 1.0 0.87 0.310.91 0.0 0.32 1.0 0.0 0.0 0.31 0.87 B Figure 46. Liquidvapor interface calculation using a volume fraction in each cell. A) Actual interface shape. B) Interf ace shape represented by the geometric reconstruction scheme using piecewiselinear interpolation. PAGE 132 132 n ut Figure 47. Calculation of volume flux thr ough the computational face (Welch and Wilson, 2000). ,,lininuT wq Figure 48. Physical model of a microcha nnel used in a computational simulation. PAGE 133 133 A B Figure 49. Grid system and boundary conditions used to simula te twodimensional twophase flow boiling model. A) Grid meshes with boundary conditions. B) Fine grid meshes in a microchannel for Hch = 30 m. Figure 410. Computational proc edure including mass transfer a nd phase change calculations. PAGE 134 134 15.215.616.016.416.8 391.6 392.0 392.4 392.8 393.2 Wall temperature (K)Axial direction (mm) mass flow rate = 3.448x106kg/s, ReDh=20 Hch=1mm, q=500kW/m2t=67ms Boiling positionA 15.215.616.016.416.8 384.0 384.5 385.0 385.5 386.0 386.5 387.0 Wall temperature (K)Axial direction (mm) mass flow rate = 3.448x106kg/s, ReDh=26 Hch=700m, q=500kW/m2Boiling position t=58.5msB Figure 411. Wall temperature prof iles at incipient boiling for di fferent channel heights. A) Hch = 1 mm. B) Hch = 700 m. PAGE 135 135 73.073.574.074.575.0 380.0 380.5 381.0 381.5 382.0 382.5 383.0 t=53.8ms mass flow rate = 3.448x106kg/s, ReDh=33 mass flow rate = 1.724x105kg/s, ReDh=162Wall temperature (K)Axial direction (mm) t=52ms Boiling position Hch=500m, q=500kW/m2 Figure 412. Effects of mass flow rate on wall temperature profile during incipient boiling at Hch = 500 m. 73.073.574.074.575.0 378.0 378.5 379.0 379.5 380.0 380.5 381.0 381.5 382.0 t=382.9ms q=387kW/m2 q=500kW/m2Wall temperature (K)Axial direction (mm) mass flow rate = 3.448x106kg/s, ReDh=45 Boiling position t=51.1ms Hch=300m Figure 413. Effects of heat flux on wall temp erature profile during incipient boiling for Hch = 300 m. PAGE 136 136 15 30 45 60 75 376.0 376.2 376.4 376.6 376.8 377.0 Wall temperature (K)Axial direction (mm) Hch=30m, q=500kW/m2mass flow rate = 3.448x106kg/s, ReDh=92 t=16ms Figure 414. Wall temperature fluctuations while bubbles generate everywhere in a microchannel for Hch = 30 m. A B Figure 415. Onset of nucleate boiling on a heat ed surface at different mass flow rates for Hch = 1mm. A) Upstream region at t = 67 ms for m = 3.4486 kg/s. B) Downstream region at t = 68 ms form = 1.7245 kg/s. PAGE 137 137 A B C D E F G H I J K L Figure 416. Bubble departure pro cedure in a microchannel with Hch = 700 m at q = 500 kW/m2, m = 1.7245 kg/s and Re = 126 for different simulation times. A) t0 + 0 ms, B) t0 + 30 ms, C) t0 + 30.5 ms, D) t0 + 31.1 ms, E) t0 + 40 ms, F) t0 + 90 ms, G) t0 + 119 ms, H) t0 + 121 ms, I) t0 + 121.2ms, J) t0 + 121.3 ms, K) t0 + 122 ms and L) t0 + 126 ms. PAGE 138 138 A B C D E F G H I J Figure 417. Bubble departure pro cedure in a microchannel with Hch = 500 m at q = 500 kW/m2, m = 1.7245 kg/s and Re = 162 for different simulation times. A) t0 + 0 ms, B) t0 + 6 ms, C) t0 + 46 ms, D) t0 + 56 ms, E) t0 + 65 ms, F) t0 + 66.5 ms, G) t0 + 67 ms, H) t0 + 67.1 ms, I) t0 + 69 ms and J) t0 + 71 ms. PAGE 139 139 A B t t+51ms C D t+51.5ms t+51.6ms E F t+51.7ms t+51.8ms G H t+51.9ms t+56ms Figure 418. Bubble departure pro cedure in a microchannel with Hch = 300 m at q = 500 kW/m2, m= 1.7245 kg/s and Re = 223 for different simulation times. A) t0 + 0 ms, B) t0 + 51 ms, C) t0 + 51.5 ms, D) t0 + 51.6 ms, E) t0 + 51.7 ms, F) t0 + 51.8 ms, G) t0 + 51.9 ms and H) t0 + 56 ms. 0.20.40.60.81.01.2 1 10 100 ReDh=460 ReDh=223 ReDh=162 ReDh=127 Departure diameter (m)Channel height (mm) mass flow rate = 1.724x105kg/s q=500kW/m2ReDh=96 Figure 419. Bubble departure diameter variations for different cha nnel heights at same heat flux and mass flux. PAGE 140 140 A B C D Figure 420. Contour plots of vapor volume fraction at va rious simulation times for Hch = 1 mm at q = 500 kW/m2 and m = 3.4486 kg/s. A) Isolated small bubbly flow at t = 500 ms, B) Bubblebubble interac tion at t = 580 ms, C) Bubbly flow with coalescence at t = 600 ms and D) Bubbly and elongated sl ug/semiannular flow at t = 670 ms. A B C D E F G Figure 421. Contour plots of vapor volume fraction at va rious simulation times for Hch = 30 m at q = 500 kW/m2 and m = 3.4486 kg/s. A) Nucleate boiling at t = 16 ms, B) Bubbly flow at t = 21 ms, C) Bubbly flow at t = 22.4 ms, D) Elongated slug flow at t = 22.9 ms, E) Bubbly flow at t = 24 ms, F) Churn flow at t = 26 ms and G) Elongated slug flow at t = 31 ms. PAGE 141 141 74.6974.7074.7174.7274.7374.7474.7574.76 400 402 404 406 408 410 412 t=22.9ms Axial direction (mm) mass flow rate = 3.448x106kg/s, ReDh=92 Hch = 30m, q = 500kW/m2 Wall temperature (K)A 74.6974.7074.7174.7274.7374.7474.7574.76 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 mass flow rate = 3.448x106kg/s, ReDh=92 Hch = 30m, q = 500kW/m2 Liquid velocity (m/s)Axial distance (mm)B Figure 422. Wall temperature a nd velocity profile underneath a bubble with local dryout near the downstream region for Hch = 30 m at q = 500 kW/m2 and m= 3.4486 kg/s. A) Wall temperature variations. B) Liquid velocity variations underneath a bubble. PAGE 142 142 74.8074.8274.8474.8674.8874.90 402 403 404 405 406 407 mass flow rate = 3.448x106kg/s, ReDh=92 Hch = 30m, q = 500kW/m2t=22.9ms Axial direction (mm) Wall temperature (K)A 74.8074.8274.8474.8674.8874.90 0.5 1.0 1.5 2.0 2.5 3.0 3.5 mass flow rate = 3.448x106kg/s, ReDh=92 Hch = 30m, q = 500kW/m2Liquid velocity (m/s)Axial direction (mm) B Figure 423. Wall temperature and velocity profile underneath a slug bubble near the downstream region for Hch = 30 m at q = 500 kW/m2 and m= 3.4486 kg/s. A) Wall temperature variations. B) Liquid velo city variations un derneath a bubble. PAGE 143 143 A B Figure 424. Contour and vector plots of vapor bubbles in a microchannel with Hch = 300 m for q = 500 kW/m2 and m = 1.7245 kg/s at t = 375 ms. A) Vapor volume fraction. B) Vector plot around bubbles. 0.0 0.2 0.4 0.6 0.8 20 40 60 80 100 120 q = 500 kW/m2 Hch = 700 m Hch = 300 m Hch = 30 m Heat transfer coefficient (kW/m2K)Quality mass flow rate = 3.448x106 kg/s Figure 425. Heat transfer coefficient variations with quality for differe nt channel heights at q = 500 kW/m2 and m= 3.4486 kg/s. PAGE 144 144 0 00 20 40 60 8 0 20 40 60 80 100 q = 500 kW/m2 q = 200 kW/m2Heat transfer coefficient (kW/m2K)Quality Hch = 300 m mass flor wate = 3.448x106kg/s Figure 426. Heat transfer coefficient variati ons with quality for diffe rent heat fluxes at Hch = 300 m and m= 3.4486 kg/s. PAGE 145 145 Table 41. Dimensions of the microcha nnel used in computational simulation Wch, mm Hch, mm Hs, mm H p mm Lch, mm L p mm 231.0 1.0, 0.7 0.5, 0.3 0.03 3.0 15.0 60.0 15.0 Table 42. Thermophysical prope rties of working fluid and heat sinks used in computational simulation Water, Tin = 303.15 K, Tsat = 373.15 K l, v, kg/m3 958.35, 0.5982 hlv, kJ/kg 2256.41 cp,l, cp,v, J/kgK 4125.67, 3768.16 kl, kv, W/mK 0.6791, 0.0251 ,lv Pas 42.817910 51.22710 N/m 0.05891, 0.05696 Copper s, kg/m3 8978 cp,s, J/kgK 381 ks, W/mK 387.6 Table 43. Grid density for each cas e used in computational simulation Channel height (mm) Inlet plenum Channel Outlet plenum 1.0 157,1520,157 10020,520,49920 3818,3820,3818 0.7 157,1520,157 10020,520,49920 3818,3820,3818 0.5 157,1520,157 10020,520,49920 3818,3820,3818 0.3 157,1515,157 10015,515,49915 3818,3815,3818 0.03 157,1510,157 2000010 4020,4010,4020 PAGE 146 146 Table 44. Effects of cha nnel height on nucleate boiling at q = 500 kW/m2 and m = 3.4486 kg/s Channel height (mm) Simulation time (ms) Boiling position (mm) Re 1.0 67.0 15.9, 16.5 20 0.7 58.5 16.0 26 0.5 53.8 74.8 33 0.3 51.1 74.9 45 0.03 16.0 15 75 92 PAGE 147 147 CHAPTER 5 TWOPHASE FLOW INSTAB ILITIES IN A TRAPEZ OI DAL MICROCHANNEL Fluid flow and heat transfer in microscal e channels are used in many applications, including microchannel heat sinks micro heat exchanger, and microfluidic devices. Especially, twophase microchannel heat sinks are widely studied due to their very high heat removal performance with a large surface area to volume ratio, compact dimensions, and low flow rate requirements. For these reasons, numerous studies of twophase flow boiling in microchannels have been extensively investigat ed in the last decades; the studies, however, were mainly limited to experimental works (Bowers and Mudawar, 1994a; Coleman and Garimella, 1999; Zhao and Bi, 2001; Qu and Mudawar, 2003a and 2003b; Kuo et al., 2006). Hence, there is a need to numerically analyze twophase flow boiling through a microchanne l as considered in Chap. 3 and Chap. 4, which would help us better analyze the available experimental data and also understand the heat transfer mechanism of twophase flow boiling such as vaporization through the liquidvapor interface in liq uid thin films and bubble dynamics dealing with nucleate boiling and bubblebubble interaction. The flow pattern transition from bubbly flow to slug flow causes undesired flow instabilities. Twophase flow instabilities lead to dryout and mechanical failure due to wall temperature fluctuations relate d with flow pattern transition. It can be also observed in microchannels as in macrochannels. The flow inst abilities are very complicated because several effects such as density wave oscillation, pr essure drop and temperature fluctuation occur simultaneously during flow boiling (Boure et al ., 1973; Tadrist, 2007; Kakac and Bon, 2008). Recently, these flow instab ilities in microchannels have been studied by several authors (Brutin et al., 2003; Brutin and Tadrist, 2004; Qu a nd Mudawar, 2003a and 2004; Wu and Cheng, 2004; Hetsroni et al., 2006; Muwanga et al., 2007; Wang et al 2007; Huh et al., 2007). Their PAGE 148 148 investigations are quite constrai ned in experimental works and mo re comprehensive study is still needed to investigate the fundame ntal issues of flow instabil ities to provide safe operating conditions of twophase flow boiling in microchannels. In this chapter, threedimensional numerical modeling of flow boiling in a trapezoidal microchannel has been discussed to analyze flow instabilities. Through the current study, an alternative predictive method to avoid an experi mental error which can occur by experimental works will be provided to determine stable and unstable flow boiling regimes. 5.1 TwoPhase Flow Instabilities 5.1.1 Static and Dynamic Flow Instabilities Twophase flow instabilities leading to mechanical damage and system control problems are observed in many thermal management systems using forced or natural convection boiling in a heated channel. Their mechanis tic behavior is very complicate due to simultaneous interactions between channel geometries and flow conditions such as subcooled liquid, mass flux, heat flux, and inlet pressure (Qu and Muda war, 2003a; Wu and Cheng, 2003c and 2004). Flow instabilities are of important factors to provi de thermal safety limits for su ch systems under certain operating conditions. Flow instabilities are categoriz ed as static and dynamic inst ability (Bergles, 1981). Static instabilities are related to the steadystate characteristics of boiling system. It may change the boiling system to a different steadystate operati ng condition or lead to a cyclic oscillating behavior relative to flow excursion, boiling cr isis, and flow pattern transition. The periodic oscillation is caused by a small perturbation. Thus the steadystate system tends to be unstable under specific flow conditions. Dynamic instabi lities are caused by simultaneous interaction between density change, mass flux, and pressure drop. Rapid vapor bubbl e generation in boiling system causes oscillations in relation with dynamic instabili ties: density wave, pressure drop, PAGE 149 149 thermal, and acoustic. Twophase flow instabiliti es are briefly reviewed in this section. More detailed information for flow in stabilities can be found in review papers of Boure et al. (1973), Ishii (1976), Bergles (1981) and Kakac and Bon (2008) for macrochannel s, Tadrist (2007) for narrow channels, and Cheng and Wu (2006) for minichannels and microchannels. 5.1.1.1 Static instabilities Ledinegg (1938) first observed flow instabil ity in twophase boiling systems. Ledinegg instability, flow excursio n, is characterized by a sudden change in low mass flow rates. It may observe when the slope of pressure drop (channel characteristic curve) is negative and steeper than external characteristic curve with increasi ng mass flow rate as shown in Figure 51. The portions having negative gradient of the pressure drop can be unstable. More pressure drop to sustain the flow results from a sm all negative disturbance to a lower flow. It leads to a shift to point p when the flow rate decreases. The part of minimum pressure drop on the channel characteristic curve is referred to as the onset of flow instability (OFI) (Ghiaasiaan, 2008). Further decrease in the mass flow rate may lead to an increas ing pressure drop beyond the OFI point. Steadystate operation require s that the external and channel characteristic values be the same: p pp Points p and p are stable because a perturbation in the flow rate at these points cause an imbalance between th e channel and external characteri stic curve. It tends to bring the system back to the original steady state. Point p is unstable. When the system operates at point p a small positive perturbation leads the sy stem to the stable and steady condition p ; a small negative perturbation in the fl ow rate, however, moves to the point p the stable and steady condition. Kakac and Bon (2008) and Gh iaasiaan (2008) proposed the slope of the external characteristic curve be stee per than that of the internal a nd/or the external characteristic curve is nearly a vertical line to avoid Ledinegg instability. PAGE 150 150 Flow pattern transition instabilit y results from cyclic flow patt ern transition and flow rate variations due to mainly different pressure dr op characteristics caused by various flow patterns (Boure et al.1973; Kakac and B on, 2008). Flow pattern may change from bubbly flow to slug/annular flow when the void fraction increases if the flow rate reduces slightly from the operating point due to a result of a small perturbati on. The pressure drop in the slug/annular flow is less than that in bubbly flow. The flow ma y be caused to increase and the void fraction may reduce. As a result, the flow pattern changes ag ain to bubbly flow. It is similar to the flow excursion referred to Ledinegg inst ability. Nayak et al. (2003) studi ed the flow pattern transition instability in a natural circulation heavy water re actor. They developed a mathematical model to analyze the flow pattern transiti on instability and used the drift flux model to estimate the void fraction in the twophase region. Their results sh owed that the instabil ity occurs at a higher power than the Ledinegg instability and increases with a decrease in pressure or an increase in inlet subcooling. 5.1.1.2 Dynamic instabilities Dynamic instabilities occur if there is a sufficient interaction and delayed feedback between thermal or hydrodynamic inertia and compressibility of the twophase mixture or other feedback effects between flow rate, pressure drop and density change as a result of the vapor generation in a heated channel. Static instabilitie s may be on the threshold of these instabilities. The dynamic instabilities can be classified as : density wave oscillations, pressure drop oscillations, thermal oscillations, acoustic oscillations. Density wave oscillations are among the most common instabilities in twophase flow boiling channels. These oscillations occur due to a result of phase lag and feedback among the flow rate, the pressure drop, and the vapor gene ration rate. When boiling takes place in heated channels, a decrease in mass flow rate leads to an increase in en thalpy of the fluid and in turn a PAGE 151 151 decrease in the fluid density due to the increasing residence time in the heated channels. The propagation time of the pressure dr op disturbance is much faster than that of the fluid enthalpy and density. As a result, the phase lag will be induced between the pressure drop and the inlet flow. It can lead to the oscillations of the inlet mass flow rate. The characteristics of the oscillations show low frequencies ( ~ 1 Hz); it is related to the residence time of the fluid in the channel (Tardist, 2007). The interaction between the heated channel and the compressible volume of the heated wall can lead to pressure dr op oscillations (Stenning, 1964). Boure et al. (1973) proposed pressure wave oscillation were observed in s ubcooled boiling, bulk boiling, and film boiling. They proposed that induced flow oscillations may lead to mechanical damages or system control problems. These oscillations occur when the pres sure drop decreases with increasing flow rate and flow excursion is on the thre shold of them. The frequencies of these oscillations are very low ( < 0.1 Hz). Thermal oscillations can occur when the wall temperature highly fluctuates due to flow pattern transition (Stenning, 1964; Kakac et al., 1 990). These oscillations lead to largeamplitude wall temperature oscillations and in turn dryout a nd boiling. In addition, th ey can deteriorate the local heat transfer performance causing wall temperature fluctuati on or critical heat flux. Kakac et al. (2003) studied the effect of inlet subcooling on twopha se flow dynamic instabilities in a vertical and horizontal channel ha ving 7.5 mm I.D. They observed that the thermal oscillations accompany the pressure drop oscillations and both are in phase; the magnitude of the pressure drop oscillations is always larger than that of the wall temperature oscillations. 5.1.2 Microchannel Forced Convection Flow Boiling Instabilities Twophase flow instabilities in macrochannels were extensively st udied by a number of researchers in the past decades. Yadigaroglu (1981) pointed out that twopha se flow instabilities PAGE 152 152 induced by unstable flow patterns, particularly slug flow, may amplify fl ow disturbances and lead to a boiling crisis and mechanical damage. These instabilities may also occur in minichannels and microchannels as the channel decreases and the imposed heat flux increases under flow boiling. Table 5.1 summarizes experi mental investigations for twophase flow instabilities of boiling in mi nichannels and microchannels. Dynamic flow boiling instability in microchannels has been experimentally studied by several researchers in recent years. Kandlikar (2 001) observed the flow patterns of boiling in six rectangular parallel cha nnels with 60 mm long and 1 mm hy draulic diameter. Severe pressure oscillations with 1 Hz frequencies were invest igated when boiling in the channels takes place. Reverse flows were also observed due to growth of slug bubbles in the microchannels. Hetsroni et al. (2002) performed expe riments on flow boiling of wa ter in 21 silicon triangular microchannels with a base of 250 m. They investigated peri odic pressure drop and outlet temperature oscillations as a re sult of vapor generation in the pa rallel microchannels. They noted that the fluctuations of the wall temperature corres pond to that of the pressure drop. Their results, also, showed that the amplitude and frequenc y of the pressure drop and wall temperature oscillations increases while th e heat flux and the vapor quality increase. Qu and Mudawar (2003a and 2004) observed experimentally pressure drop oscillations with largeamplitude and parallel channel instabilities in 21 parallel rectangular microchannels with 231 m wide and 713 m high. They identified severe pressure drop oscill ation and mild parallel channel instability. The results showed the severe pressure drop caused by interaction between vapor generation and compressible volume in the flow loop upstream is large periodic and large amplitude oscillations. The mild parallel channel instability produced mild flow fluctuations. The instability is the result of density wave oscillation within each channel and feedback inte raction between channels. They PAGE 153 153 noted that the severe pressure dr op oscillations may trigger prematu re critical heat flux. Brutin et al. (2003) investigated twophase flow instabilities in a rectangular microchannel using npentane with a hydraulic diameter of 889 m. They observed two types of behavior: a steady state characterized by pressure drop fluctuations with low amplitudes (from 0.5 to 5 kPa/m) and with no characteristic frequency; an unsteady st ate with high amplitudes (from 20 to 100 kPa/m) and frequencies ranging from 3.6 to 6.6 Hz. Th e steady and unsteady behaviors are governed by the two control parameters: heat flux and mass ve locity. The influence of the inlet condition on twophase flow instability in a rectangular minichannel of 0.54 mm2 was observed by Brutin and Tadrist (2004). They found steady and unste ady behaviors depend on the upstream boundary conditions: a constant velocity at the minichannel entrance and at the syringe outlet. The pressure loss fluctuations with relatively small amplit udes and higher frequencies between 6 and 18 Hz are observed, when a damper is not connected. Wu and Cheng (2003c and 2004) carried out expe rimental investigation on flow boiling of water in parallel trapezoidal silicon microc hannels with a hydraulic diameter of 82.8 m, 158.8 m, and 185.6 m, respectively. They first observed th ree kinds of oscillatory flow boiling modes in microchannels: liquid/tw ophase alternating flow (LTA F) at low heat flux and high mass flux, continuous twophase flow (CTF) at medium heat flux and low mass flux, and liquid/twophase/vapor alternating flow (LTVAF) at high he at flux and low mass flux. For the LTAF case, a bubbly flow pattern was dominant and a period of oscillations was 15.6s. The amplitudes of wall temperature and inlet pressure fluctuations were in the range of 20 44 C and 5 10 kPa, respectively. The oscillations of wall temperature and pressure were in phase, while those of pressure and mass flux were out of phase. For the CTF case, the oscillation period was 32s and the oscillation amplitudes of wall temp erature and pressure were to 723 C and 8 PAGE 154 154 kPa. The oscillations of inlet pressure a nd mass flux were in phase. The LTVAF period was approximately 53 s and the amplitudes of wall temp erature fluctuations changed from 66 C to 230 C. The oscillation amplitude of pressure wa s 44 kPa. The oscillations of pressure and mass flux were out of phase, while the oscillations of wall temperatur e and pressure were in phase. The oscillation behavior of wall temperature a nd inlet pressure showed a large amplitude and long period. They pointed out that the oscillation periods depend significantly on heat flux and mass flux. Hetsroni et al. ( 2005) investigated experimenta lly water boiling phenomena with periodic wetting and dryout in 21 triangul ar silicon microcha nnels a base of 250 m. They pointed out that occurrence of singlephase liquid and twophase liquidvapor flow leads to pressure drop fluctuations. Their results show ed dryout occurred afte r venting of elongated bubble due to very rapid expansion. Xu et al. (2005) studied the onset of flow instability (OFI) and dynamic unsteady flow in 26 re ctangular microchannels with 300 m 800 m. The working fluids were deionized water and methanol The onset of flow instability was observed at the outlet temperature of 93 96 C and two ty pes of oscillations were identified: large amplitude/long period oscillation (LALPO) superimposed with small amplitude/short period oscillation at the lower inlet liquid temperatur e and small amplitude/short period oscillation (SASPO) at the higher inlet temperature. They pointed out that thermal oscillations lead to two types of oscillations. More recently, Huh et al. (2007) carried out ex perimental investigations on water flow boiling in a horizontal recta ngular microchannel having a hy draulic diameter of 103.5 m, a width of 100 m, a height of 107 m, and a length of 40 mm. They measured three types of oscillations: wall temperature, pressure drop, a nd mass flux with a very long period (100 200 s) and large amplitude. These periodic fluctuations corresponded closely to the two flow pattern PAGE 155 155 transitions: a bubbly/slug flow and an elongated slug/semian nular flow. They noted that the instability of the flow pattern transition cause a cyclic behavior of wall temperature, pressure drop, and mass flux. Wang et al. (2007) obser ved flow boiling instab ilities of water in microchannels to identify stable flow boiling regime. The experiments were carried out with six parallel silicon microchannels and a single microchannel with a hydraulic diameter of 186 m and a length of 30 mm. They found stable and unstable flow boiling regimes depending on the mass flux. A flow boiling map in terms of heat fl ux with mass flux was proposed to show stable flow boiling and unstable flow boiling regime. Th eir results showed long period oscillations and short period oscillation in temperature a nd pressure. Chang and Pan (2007) studied experimentally the twophase flow instability for water boiling in 15 parallel rectangular microchannels with a hydraulic diameter of 86.3 m. Reversed slug or annular flows were alternatively observed in every channel for unstable flow boiling. They proposed a stability map of inlet subcooling number versus phase change number and considered the magnitude of pressure drop oscillations as a trigger mechanism of the revers ed flow. Muwanga et al (2007) observed flow boiling oscillations in two different silicon mi crochannel heat sinks with 269 m 283 m using water as working fluid; one is a straight standard microchannel with 45 straight parallel channels and the second contains crosslinks between the channels. Their results showed flow oscillations were similar in characteristic trends between the two configurations. The amplitude of inlet and outlet temperature fluctuat ions were up to 15 C (p eak to peak) and 21 C (peak to peak), respectively. The amplitude of in let pressure fluctuations were up to 12.5 kPa (peak to peak). The frequency of inlet temperatur e oscillations were dependent on heat flux, flow rate, and inlet subcooling. The amplitude of the inlet pressure fluctuations increased with increased heat flux and increased flow rate. PAGE 156 156 Configurations of inlet and outlet to reduce fl ow instabilities in mi crochannels have been studied by Kosar et al. (2005) a nd Wang et al. (2008). Kosar et al. (2005) observed geometrical effects of inlet orifices in parall el microchannels with a width of 200 m and a depth of 264 m to suppress flow instabilities. Wang et al. (2008) investigated influences of inlet and outlet configurations on flow boiling inst abilities in parallel microchannels with a hydraulic diameter of 186 m. They used three type configurations: TypeA with restriction for flow inlet and outlet, TypeB without restriction for fl ow inlet and outlet, and TypeC with restriction for flow inlet and without restriction for flow outlet. They point ed out that steady flow boiling prevails in the parallel microchannels with TypeC connection. The studies of flow boiling inst abilities in microchannels have been extensively relied on experimental works and only few studies were performed in a single microchannel (Huh et al., 2007; Wang et al., 2007) in recent years. There ha ve been no predictive to ols to investigate the flow boiling instabilities in th e open literatures as well. Although some literatures have shown experimentally the flow instability phenomena in parallel microchannels and in a single microchannel, the understanding of flow instabili ties and related heat transfer mechanisms in a microchannel is still very limited due to complex phenomena under flow conditions. 5.1.3 Research Objectives The purpose of the current study is to perform a direct numerical si mulation using a threedimensional model which expands the twodime nsional model of Chap 4 to observe flow instabilities occurred during twophase flow boiling. Twophase flow instabilities are complicated phenomena with several coupled eff ects under flow conditions. These instabilities affect harmfully the thermal systems. The stable and unstable ope rating conditions are PAGE 157 157 investigated to provide a better understanding of flow instabilities when boiling in a microchannel takes place. Thus, the objectiv es of current work are as follows: Develop a threedimensional numerical simu lation model to investigate flow boiling instabilities occurred in a single trapezoidal microchannel. Study the relationship between instability of fl ow boiling with system parameters such as pressure and system temperature with a fi xed inlet temperature and hydraulic diameter. Investigate pressure and wall temperature osci llations with various heat flux and mass flux conditions. The obtained results through numerical simulati ons will be qualitatively and quantitatively compared with existing experimental data available in the literature. 5.2 Problem Descriptions With the large ratio of length to hydraulic di ameter, microchannel heat sinks are used to remove effectively higher heat fluxes from hi gh power density electronic systems. Figure 52 illustrates the geometry of a trapezoidal microchanne l used to simulate flow instabilities of twophase flow boiling. The geometry and dimensions of the microchannel co rrespond exactly to the experimental apparatus conducted by Wang et al. (2007). The width of top and bottom, and depth of the trapezoidal microchannel are 427 m, 208 m, and 146 m, respectively. The angle at the base is 53 degree, the hydraulic di ameter and length of the channel are 186 m and 38 mm with inlet and outlet unheated region, respectively. The dimensions of the microchannel used in the computational simulation are given in Ta ble 52. The thermal boundary condition at the bottom is the constant heat flux boundary condition and assumes the adiabatic boundary condition at the top and two side walls as shown in Figure 52A. Referring to Figure 52B, wall heat flux, wq is imposed on the microchannel heat ed inside area and is obtained as: 2effcell w botcqW q NWH (51) PAGE 158 158 where effq is the imposed heat flux on th e bottom of the microchannel, cellW and botW are the width of the unit cell and mi crochannel, respectively, cH is the height of the trapezoidal microchannel heated, N is the number of the microchannel. Figure 53 shows the computa tional domain used in the threedimensional numerical simulation for flow boiling. Taking advantage of symmetry, half a channel is chosen for the current numerical simulation as illustrated in Figu re 53A. Water is used as the working fluid and inlet temperature, Tin of 99.85 C (373.0K), flows through th e inlet plenum with an adiabatic condition. The subcooled liquid approaches up to saturation temperat ure and vapor bubbles generates on the heated wall. Twophase mixt ure flows out of the outlet plenum with Tout and pout as shown in Figure 53B. The heat fl uxes in the rage of 297.8 495.4 kW/m2 are imposed on the heated wall and the mass fluxes of 679.5 3523.2 kg/m2s flow at the inlet region. The pressure outside of the outlet pl enum is assumed to be the atmospheric pressure. Wall temperatures, Tw1, Tw2, Tw3, Tw4, and Tw5, are obtained on the heated wall along the flow direction as depicted in Figure 53B. 5.3 Computational Model The 3D model geometry was created and meshed using the GAMBIT software. The grid of half a flow channel has been generated us ing hexahedral mesh el ements for a 3D flow channel geometry as shown in Figure 54 and the grid density used in the current simulation is listed in Table 53. The boundary conditions us ed to obtain computational results are listed follows: Symmetric condition is applied because only half a cha nnel is modeled to reduce computational time for this intensive simulation studies. The velocity and temperature of subcooled liquid are specified at the liquid inlet. PAGE 159 159 The primary and secondary phase is set to liquid and vapor phase, respectively. The volume fraction of vapor phase v is equal to zero because the pure liquid flows through the inlet. At the outlet plenum the pressure as gauge pres sure is set to zero: the gauge pressure is internally used in the FLUENT software. Ab solute pressure is the sum of the operating pressure and the gauge pressure. In this mode l, the default value of 1 atm (101.325 kPa) is assumed as the operating pressure and the absolu te pressure at the out let is assumed as 1 atm. At all channel walls, no slip boundary condition is imposed. Surface tension effects are considered and the contact angle is specified to adjust the curvature of the surface near the wall. Constant wall heat flux is imposed at the bot tom and side wall of the heated channel. The top of the channel and in let and outlet plenum are assumed to be an adiabatic condition. Mass transfer through the liquidvapor inte rface accounting for phasechange effects are included during the solution procedure. A volume of fluid model (VOF) was chose to an alyze the flow instabilities of twophase flow boiling with the motion of large bubbles in liquid phase. The VOF technique computes a single set of momentum equations and determines the interface shape of the two immiscible fluids by tracking the volume fraction of each phase in each comput ational cell. In this study, vapor phase was tracked as the secondary phase and the volume fraction of the vapor phase was obtained from the continuity e quation. UDFs (User Defined Functi ons) are used to take account for phasechange between the liquid and vapor phas e. A segregated unsteady solver was used in the simulation due to having a robustness to so lve twophase flow boiling with a large bubble motion. The density ratio of liquid to vapor is ve ry large in twophase flow boiling. It causes the discontinuity of pressure gradient with rapidly changing densities at the liquidvapor interface. For the pressure interpolation, the PRESTO (p ressure staggering option) scheme was used to deal with large pressure gradie nts at the cell face. For the mo mentum and energy equations, the PAGE 160 160 QUICK (Quadratic Interpolati on for Convective Kinetics) tec hnique was adopted. With thirdorder accuracy, the QUICK scheme can minimize th e resultant artificial diffusion errors that may occur when using loworder upwind scheme. The PISO (pressureimplicit with splitting of operators) scheme based on the SIMPLE algorithm was used for the pressurevelocity coupling with an additional corrector level that satisfi es both continuity and momentum equations in contrast with the SIMPLE scheme. It has faster convergence of the solution for every timestep and robustness. The solution procedure is similar to that given in Chap. 4. The simultaneously computational procedure between the momentum and energy equations for threedimensional case requires much more computational times than for the twodimensional case in Chap. 4. The simultaneous computations of the energy and mo mentum equations for threedimensional case require a large number of iterations and three an d half months were required to obtain the computational results. In addi tion, due to insufficient compu ting power, high inlet temperature was used in this simulation. It leads to less simulation times than lower inlet temperature. 5.4 Results and Discussion Flow instabilities in a single trapezoidal mi crochannel have been i nvestigated for three different cases as listed in Ta ble 54 for various wall heat fluxes and working fluid mass fluxes. Case I (effq = 297.8 kW/m2 and G = 3523.2 kg/m2s) represents low heat flux and high mass flux; Case II (effq = 307.6 kW/m2 and G = 711.8 kg/m2s) is for medium heat flux and medium mass flux; and Case III (effq = 495.4 kW/m2 and G = 679.5 kg/m2s) stands for high heat flux and low mass flux. PAGE 161 161 5.4.1 Flow Boiling Patterns Flow boiling patterns predicted in the middle of the trapezoidal microchannel for Case I are shown in Figure 55. When a liquid flow with a very small degree of subcooling enters a heated channel, nucleate boiling in itiates very quickly near the downstream of the inlet as the heated channel wall temperature exceeds the amount of superheat required for boiling incipience. As shown in Chap. 4, the threshold of nucleat e boiling strongly depends on the imposed heat flux and mass flux at a given inlet temperature and mass flux. Thus, vapor bubbles grow as nucleate boiling is initiated along the cha nnel and then depart from the heated surface when the bubble size is large enough. As more and more isolated single bubbles are entrai ned in the flow, the flow pattern becomes the bubbly flow at t = 12 ms as depicted in Figure 55A. At t= 15 ms, the bubbles receive more heat from the surrounding superheated liquid, then grow larger, and coalesce with neighboring bubbles as shown in Figur e 55B. As a result, the flow pattern evolves from bubbly flow to slug flow. At t = 33.86 ms, the bubble slugs receive more heat and grow more along the channel direction such that the flow pattern becomes an elongated slug flow as shown in Figure 55C. The flow patterns in the microchannel display a se quential and periodic behavior: nucleate boi ling, bubbly flow, bubbly/slug flow, and elongated slug flow as shown in Figure 55D. Figure 56 shows flow boiling pa tterns predicted in the middl e section of the microchannel for Case II. Generally, the flow patterns of Case II are similar to those of Case I for some portions of the flow. After boiling incipience, the flow pattern is also th e bubbly flow (at t = 5 ms) as shown in Figure 56A. Bubbl y flow transition to slug flow transition at t = 7 ms and then to elongated slug flow at t = 14.92 ms are provid ed in Figure 56C and E. Periodic flow boiling pattern, also, is observed for Case II (Figure 56G ). The flow patterns of Case II obtained from the numerical simulation are consistent with the experimental observations by Wang et al. PAGE 162 162 (2008). They observed isolated bubbles, bubble co alescence, and elongated bubbles as illustrated in Figure 56B, D and F. The transition from bubbly flow to slug bubble flow resulted from expansion and coalescence trigge rs flow instabilities in the microchannel. The elongated slug bubble prevents liquid mass from flow ing through the flow channel. It causes the inlet pressure and wall temperature to increase and then oscill ate. The increased inlet pressure forces the elongated slug bubble to move toward the channel exit. Small bubbles occur in the channel and the wall temperature decreases due to a rewetting process. It shows period ic behavior as shown in Figures 55D and 56G. The flow instabilities cause mechani cal damage and operating control problem due to local dryout and periodic wall te mperature fluctuations. The flow patterns for CASE III are very similar to those of CASE II, thus they will not be shown here. 5.4.2 Wall Temperature Fluctuations With shortperiod oscillations, temporal variatio ns of wall temperatures for Case I at lower heat flux and higher mass flux are depicted in Figure 57. The temperature fluctuations are mainly initiated by the transition from bubbly to sl ug flow patterns. The temporal variations of wall temperature for Case I start at t = 15 ms wh en the transition begins (as shown in Figure 57A). The results show that the amplitudes of wall temperature at Tw4 is much larger than those at Tw3. Tw3 displays a very small scale oscillati ons during the simulation period, while Tw4 shows higher oscillatory behavior. The reason is that the flow structure around Tw3 contains basically isolated bubbles (Fig. 56A) with a lower quality than that near Tw4 with slug flows as the large amplitude oscillation is basically caused by wet (liquid) and dry (v apor) alternati ng conditions on a constant heat flux heater surface. Another im portant finding is that the temporal mean temperatures for both locations are relative ly constant. Although the imposed heat flux, subcooling at the inlet, and channel geometry are moderately diffe rent, the obtained wall temperature fluctuations show qualitatively and quantitatively good agreement with the PAGE 163 163 experimental results by Qu a nd Mudawar (2004) as shown in Fi gure 57B. Their results also show that the wall temperature oscillations at Ttc4 have higher amplitudes than those at Ttc3. The quantitative discrepancy between th e current results and the experime ntal data may be due to the large difference in inlet subcooling levels They used highly subcooled liquid ( Tsub = 70 C), while lower subcooled liquid ( Tsub = 1.5 C) was used in this study. Higher subcooled liquid condenses some of the generated bubbles in the channel and before departing and the bubble resides in the channel for a longer period, while lower subc ooling leads to fast bubble departure and coalescence. The trend of the obtained results also agrees with the experimental data by Xu et al. (2005). They pointed out that periodic oscillati ons become shorter with increasing inlet temperature. Figure 57C shows bubbl y flow at t = 12 ms before st arting oscillation (Region I in Figure 57A). At t = 15 ms, bubbly to slug flow transition causes flow oscillation to start as illustrated in Figure 57D (region II in Figure 57A). Periodic boiling consisted of bubbly and elongated slug flows as depicted in Figure 57D is dominant in region III in Figure 57A. It leads to wall temperature fluctuations with a cyclic rewetting process. Temporal wall temperature fluctuations at Tw3 and Tw4 for Case II are shown in Figure 58. Comparing to Case I, Case II has a similar wall heat flux level but the liquid inlet mass flux is only one fifth of that for Case I. As a result, the flow development and transition time are much quicker that Case I as shown in Figures. 55 and 56. In contrast with the results of Case I (as shown in Figure 57A), the temporal wall temperature fluctuations at Tw3 are similar to those at Tw4. For both cases, the mean temperatures are ri sing with time. The reason is that nucleate boiling is initiated near the inlet for Case II due to much smaller liquid inlet mass flux that also causes the bubbles to expand and coalesce near Tw3 rather than to be condensed. Wang et al. (2007) reported the same results to Case II. For the flow patterns, the bubb ly flow is dominant in PAGE 164 164 Region I as shown in Figure 58B At Region II, transition from bubbly flow to slug flow leads to increasing wall temperature fluctuations. Pe riodic flow boiling with elongated slug flow causes the higher wall temperature fluctuations in Region III as depicted in Figure 58D. The peaktopeak amplitude and oscillation period of wall temperature fluctuations are 15.44 C and 1.6ms, respectively, as depicted in Figure 58A. For Case III, the wall heat flux is 61% higher than that in Case II while the liquid mass flux at the inlet is about the same. Similar to Case II (as shown in Figure 58A), the temporal wall temperature fluctuations at Tw3 are at the same levels to those at Tw4 as shown in Fig. 59. But the two are offphase to each other. For bot h cases, again the mean temperatures are rising with time. Similarly with Case II, the wall temperature fluctuations are mainly affected by the periodic flow boiling pattern. The peaktopeak amplitude an d oscillation period of wall temperature fluctuations are 17. 79 C and 1.0 ms, respectively, as shown in Figure 59. The amplitude and oscillation period of both Cases II and III show large amplitude and very shortperiod as a result of increasing heat flux and decreasing mass flux comparing to Case I. The current results are qualitatively consistent with the experimental investigation by Wang et al. (2007). In the current results, th e amplitude of wall temperature fluctuations increases and the period of oscillation becomes shorter with in creasing heat flux and d ecreasing mass flux. Wang et al. (2007) also showed that the wall temperatures at Tw3 and Tw4 increase as heat flux increases and mass flux decreases. 5.4.3 Pressure Fluctuations Figure 510 shows temporal variations of pressure fluctuations in the inlet and outlet regions for Case I. Pressure fluctuations are induced by the alternating flow pattern between bubbly and elongated slug flow as shown in Fi gure 55B and C. The following provides the physics for the fluctuating nature of the inlet and outlet pressures. The transition between the PAGE 165 165 bubbly and slug flows at t=15 ms (Fig. 55) leads to increasing inlet pressure as shows in Figure 510A. The reason is that liquid mass flow is pa rtially blocked by the elongated slug bubble in the flow channel after the bubbly/sl ug flow transition takes place and it leads to an instantaneous increase in the inlet pressure. In turn, the increased inlet pressure then pushes the elongated slug bubbles through the channel that causes the inlet pressure to decrease. The ou tlet pressure fluctuations are lower in magnit ude but still show a periodic behavior (Figure 510B). The peaktopeak amplitude and oscillation period of inlet pressure fluctuations is 6.5 kPa and 1.9 ms, respectively. It is noted that for Case I, the period of in let pressure fluctuations at 1.9 ms is very close to that of wall temperatures at 1.8 ms as shown in Figure 57A. The inlet and outlet pressure oscillations are also in phase as de picted in Figure 510C. As compared with the experimental data of Qu and Mudawar (2003a) as shown in Figure 510D, it is seen that qualitatively the simulated results have the same trends as those of the experimental data, while different subcooling levels were used. The reason for the similar trends may be that the flow patterns resemble each other in th at they all start with a bubbly fl ow and then made a transition to elongated slug flow. The temporal variations of the inlet pressure for Case II are illustra ted in Figure 511. The peaktopeak amplitude and oscillation period of inlet pressure ar e 37.32 kPa and 1.1 ms, respectively (Figure 511A). The re sults show that much higher pre ssure fluctuations than those of Case I are predicted as shown in Figure 510A. The reason is again that the generation rate of vapor bubbles for Case II is five times greater than that for Case I due to the drastic reduction in the inlet liquid mass flux as listed in Table 54. The oscillation period investigated by Wang et al. (2007) is 11.0 ms as shown in Figure 511B. The discrepancy between the current result and the experimental data by Wang et al. (2007) may be due to the di fferent inlet temperatures that PAGE 166 166 provide different inlet subcooling levels. As above mentioned, th e oscillation period is shorter with increased inlet temperatur e. Thus, the oscillation period of the current study with higher inlet temperature is much shorter than that of Wa ng et al. (2007) with a lower inlet temperature. In addition, the difference between the current results and the experimental results may also be caused by the compressible volume near upstream and the annular/mist flow observed in their experiment. Figure 512 shows inlet pressure fl uctuations for Case III. It is shown that the peaktopeak amplitude and oscillation period of inlet pressure history are 76.70 kPa and 0.73 ms, respectively, as shown in Figure 512A. The oscillation period reported by Wang et al. (2007) is 6.4 ms as illustrated in Figure 512B. Although th e discrepancy between the numerical result and the experimental data exists, both cases show that the amplitude and oscillation period of the inlet pressure is strongly affected by the increas ed heat flux. The amplitude of inlet pressure fluctuations is proportional to im posed heat flux and the oscillati on period of those is inversely proportional to mass flux as listed in Table 54. 5.4.4 Stability Map Figure 513 shows the twophase forced c onvective boiling stability map reported in Chang and Pan (2007) in terms of the subcooling number, Nsub, and the phase change number, Npch, as defined by i sub lvvh N h (52) pch lvvQ N mh (53) where ih is the inlet subcooling, lvh the latent heat of vaporization, the difference between liquid and vapor density, Q the total h eat input rate to the channels and m the total mass flow PAGE 167 167 rate to the channels. Chang and Pan (2007) constructed the stability map for a single microchannel based on the twophase thermal instab ility theory developed by Saha et al. (1976) with experimental data from Qu and Mudawar (2004) in addition to their own results. The different sections on the map are divided by a solid line and a dotted line. The solid line represents that the exit quality of the flow is equal to unity which means all the flows that fall to the left of this line are singlephase liquid flow s. The dotted line denot es the stability boundary, so the triangular area bounded by the two lines is for stable twophase flow s while the area to the right of the dotted line repres ents the unstable flows. The three cases simulated in th e current research are also placed in the map of Chang and pan (2007). Case I is in the stable regime while Cases II and III are in the unstable regime. Since the subcooling numbers of the current cases are very small due to the low inlet subcooling, however the bulk of the experimental data from Chang and Pan (2007) and Qu and Mudawar (2004) are larger due to the higher inlet subcoo ling. The phase change numbers of the current numerical results are higher because the phase change number is proportional to the imposed heat flux and inversely propor tional to the mass flux. It is important to point out that the author concludes the stability map developed by Chang and Pan (2007) is valid for the current simula ted results based on the following analysis. As mentioned in Section 5.4.2 on wall temperature fl uctuations, the mean temperatures for both Tw3 and Tw4 locations are relatively constant (Fig. 57A) for Case I while they are on a rising trend for both Cases II and III (Fig. 58A and Fig. 59). Based on th e wall temperature fluctuation patterns, Case I is considered stab le and Cases II and III are unstable. PAGE 168 168 p Flow rate p pp External characteristic Channel characteristic External characteristic OFI ONB Figure 51. Ledinegg instabilit y based on pressure difference and mass flow rate (Kakac and Bon, 2008). Lch Hcelleffq effq A Figure 52. Physical diagram of a trapezoidal microchannel. A) Dimensions of a trapezoidal microchannel. B) Dimensions of a twodimensional microchannel. PAGE 169 169 HchWtopWwWbot HwWcell effq HcellB Figure 52. Continued A Symmetry plane qw qw y z x Wtop/2 Wbot/2 B qw LchLinLout linT, linm outTout p z y 1 wT2 wT3 wT4 wT5 wT xL2L3L4L5 L1L6 Figure 53. Symmetric fl ow channel geometry used in a computational simulation. A). Half a microchannel. B) Flow region in a microchannel for a computational simulation. PAGE 170 170 Figure 54. Grid system of a threedimensional CFD model. A Flow direction B C D Flow direction Figure 55. Flow patterns in the middle of a trapezoidal micr ochannel for Case I. A) Bubbly flow at t = 12 ms, B) Transition from bubbly flow to slug flow at t = 15 ms, C) Elongated Slug flow at t = 33.86 ms, and D) Periodic flow boiling. PAGE 171 171 A B C D E F G Figure 56. Flow patterns in the middle of a trapezoidal micr ochannel for Case II. A) Bubbly flow at t = 5 ms, B) Isolated bubble (Wang et al., 2008), C) Bubblyelongated transition at t = 7 ms, D) Bubble coalescen ce (Wang et al., 2008), E) Elongated slug flow at t = 14.92 ms, F) Elongated bubble (Wang et al., 2008) and G) Periodic flow boiling. PAGE 172 172 101520253035404550 100 110 120 130 140 150 160 III IIII I Wall temperature (oC)Simulation time (ms) Tw3 Tw4 q=297.8kW/m2G=3523.2kg/m2s Start oscillation A B Figure 57. Wall temperature fluctuations with small amplitude and shortperiod oscillation and flow patterns for Case I. A) Current result, B) Experimental investigation by Qu and Mudawar (2004), C) Bubbly flow at t = 12 ms (Region I), D) Bubblyelongated transition at t = 15 ms (Region II), E) Elongated slug flow at t = 33.86 ms (Region III). PAGE 173 173 C Flow direction D E Figure 57. Continued PAGE 174 174 6 8 10 12 14 16 100 110 120 130 140 150 160 III II Wall temperature (oC)Simulation time (ms) Tw3 Tw4 q=307.6kW/m2G=711.8kg/m2s Start oscillation I A B C D Figure 58. Wall temperature fluctuations with large amplitude and shortperiod oscillation and flow patterns for Case II. A) Wall temperature fluctuati ons, B) Isolated bubbly flow at t = 5 ms, C) Transition from bubbly flow to slug flow at t = 7 ms, and D) Elongated slug flow at t = 14.92 ms. PAGE 175 175 67891 01 1 120 130 140 150 160 170 Wall temperature (oC)Simulation time (ms) Tw3 Tw4 q=495.4kW/m2G=697.5kg/m2s Start oscillation Figure 59. Wall temperature fluctuations with large amplitude and shortperiod oscillation for Case III. PAGE 176 176 A 1520253035404550 1.15 1.20 1.25 1.30 1.35 G = 3523.2 kg/m2s q = 297.8 kW/m2 Inlet pressure (bar)Simulation time (ms) 1520253035404550 0.90 0.95 1.00 1.05 1.10 1.15 1.20 Outlet pressure (bar)Simulation time (ms) G = 3523.2 kg/m2s q = 297.8 kW/m2 B C 46 47 48 49 50 0.9 1.0 1.1 1.2 1.3 1.4 pin Pressure (bar)Simulation time (ms) poutq = 297.8 kW/m2G = 3523.2 kg/m2s D Figure 510. Inlet and outlet pressure fluctu ations with small amplitude and shortperiod oscillation for Case I. A) Inlet pressure oscillations, B) Outlet pressure oscillation, C) Inlet and outlet pressure in phase and D) Parallel channel instability (Qu and Mudawar, 2003a). PAGE 177 177 6 8 10 12 14 16 1.00 1.25 1.50 1.75 2.00 Inlet pressure (bar)Simulation time (ms) q=307.6kW/m2G=711.8kg/m2sA Time (s)Inlet pressure (bar) B Figure 511. Inlet pressure fluctuations with large amplitude and shortperiod oscillation for Case II. A) Current result. B) Experime ntal data obtained by Wang et al. (2007). PAGE 178 178 567891011 1.00 1.25 1.50 1.75 2.00 2.25 2.50 Inlet pressure (bar)Simulation time (ms) q=495.4kW/m2G=697.5kg/m2sA B Figure 512. Inlet pressure fluctuations with large amplitude and shortperiod oscillation for Case III. A) Current result. B) Experime ntal data obtained by Wang et al. (2007). PAGE 179 179 0 400 800 1200 1600 0 100 200 300 400 Case III Case II Case I Current result, stable (Case I) Current result, unstable (Case II and Case III) Chang and Pan (2007), stable (singlephase) Chang and Pan (2007), stable (twophase) Chang and Pan (2007), unstable (twophase) Qu and Mudawar (2004), Tin = 30 oC (unstable) Qu and Mudawar (2004), Tin = 60 oC (unstable)NsubNpchxe = 0Stability boundary Figure 513. Stability map in subcooling numbe r versus phase change number based on Chang and Pan (2007). PAGE 180 180 Table 51. Summary of investigation on flow instabilities in minichannels and microchannels Author Geometry and fl ow conditions Remarks Kennedy et al. (2000) 23cm long tubular test section; 1.17 mm and 1.45 mm in diameter with a 16cm long heated length Working fluid: Subcooled water Tin = 50 oC; G = 800 4500 kg/m2s q = 0 40 MW/m2; pou t = 3.44 10.34 bar Investigated the onset of flow instability in uniformly heated horizontal microchannels. Generated pressure drop versus mass flow rate curves for fixed wall heat flux and channel exit pressure. Kandlikar (2002) Six parallel mi nichannels Observed flow instability w ith reversal flow. Hetsroni et al. (2002) 21 parallel triangular microchannels with a base of 250 m G = 148 290 kg/m2s; q = 32 and 36 kW/m2 pout = 1.0 bar Observed experimentally the pressure and temperature fluctuations for fixed wall heat flux and mass flux. Found the amplitude and frequency of the pressure and temperature fluctuations increases as the increasing heat flux and vapor quality. The temporal behavior of temperature fluctuations corresponds to that of pressure fluctuations. Wu and Cheng (2003) 8 and 15 parallel silicon micr ochannels of trapezoidal crosssection Working fluid: Subcooled water Dh = 158.8 and 82.8 m G = 137 251 kg/m2s; q = 5.78 13 W/cm2 Observed first largeamplitude/long period boiling fluctuations can be sustained when the fluctuations of pressure drop and mass flux have phase differences. Observed smallamplitude/shortp eriod fluctuations in inlet temperature and in instantaneous mass flux due to the bubble dynamic instabilities in the twophase flow. The fluctuation periods are found to be dependent on channel size, heat flux, and mass flux. Qu and Mudawar (2003) 21 parallel rectangular microchannels with 231713 m Working fluid: Subcooled water Tin = 30 and 60 C G = 134.9 400.1 kg/m2s q = 100.4 and 203.0 W/cm2; pout = 1.17 bar Investigated two types of twophase flow instability, namely severe pressure drop oscillations and mild parallel channel instability. Pressure drop oscillation produces fairly periodic, largeamplitude for oscillations, which are the resu lt of interaction between vapor generation in channels and compre ssible volume in the flow loop upstream of the heat sink. Parallel channel instability produces only mild flow fluctuations, which are the result of density wave oscillation within each channel and feedback interaction between channels. The severe pressu re drop oscillation may tr igger premature CHF. PAGE 181 181 Table 51. Continued Author Geometry and fl ow conditions Remarks Brutin et al. (2003) Vertical rectangular microchannels with 0.540 and 0.550 mm3 Working fluid: npentane G = 240, 480, and 960 kg/m2s q = 200 500 kW/m2 Observed a steady state characterized by pressure drop fluctuations with low amplitudes (from to 0.5 to 5 kPa/m) and no characteristic frequency; a nonstationary state of twophase flow. For unsteadyflow boiling, inlet pressure signal presents a regular pattern of high amplitude (8 kPa) and low frequency (~4 Hz and high peak amplitude 3000 rms). Derived a stability map drawn in nondimensional parameters: outlet vapor quality and Reynolds number. Wu and Cheng (2003) 8 and 15 parallel silicon microchannels Dh = 158.8 and 82.8 m Working fluid: Subcooled water G = 137 251 kg/m2s; q = 9.6 155 kW/m2 Observed four flow patterns: bubbly flow, slug flow, churn flow, and other peculiar flow patterns. Investigated first alternative slingphase liquid flow and twophase flow causing largeamplitude/long period fluctuations in fluid temperature and wall temperatures. Found the fluctuations periods de pend on channel size, heat flux, and mass flux. Brutin and Tadrist (2004) A rectangular minichanne l with of 0.5 4 mm2 Dh = 889 m Working fluid: npentane Tin = 25 C; Q = 15.7 125.6W Re = 250 9000 Investigated a twophase pressure dr op characteristic for several heat fluxes. Found a stability criterion depends on the two controlled parameters (heat flux and massflow rate). The amplitude fluctuations become higher as the inlet Reynolds number becomes smaller. Balasubramanian and Kandlikar (2005) 6 rectangular minichannels Dh = 333 m Working fluid: Subcooled water G = 120 kg/m2s; q = 208 316 kW/m2 Obtained the pressure drop at a given mass flux and for increasing heat fluxes. Showed the magnitude of the pressure drop incr eases with an increase in surface temperature for a given mass flux. Observed the peaktopeak variation for pressure drop remains almost constant for small change s in the surface temperature. Slug flow was observed to be dominant for higher surface temperatures, and in some cases the resulting back flow extended into the inlet manifo ld, causing severe fl ow maldistribution. PAGE 182 182 Table 51. Continued Author Geometry and fl ow conditions Remarks Xu et al. (2005) 26 rectangular microchannels with 300 m width 800 m depth Dh = 436 m Working fluid: Subcooled water and methanol Tin = 30, 50, and 70 oC; G = 20 1200 kg/m2s Q = 100 450 W Observed fine bubbles are found close to the exit of the microchannels when the OFI appears. Investigated the OFI moved in the di rection of larger mass flux with increasing inlet temperature and power input. Found the onset of flow instability occurs at the outlet temperature of 93 96 Identified two types of oscillati ons: large amplitude/long period oscillation (LALPO) s uperimposed with small amplitude/short period oscillation, small ampl itude/short peri od oscillation (SASPO), and thermal oscillation. Observed temperature and pressure oscillations were caused by the liquidvapor interface moving upst ream and downstream in the microchannels. Hetsroni et al (2006) 5 parallel triangular microchannels Dh = 100220 m Working fluid: Water and ethanol G = 32 200 kg/m2s; Q = 120 270 kW/m2 Observed oscillation freq uency is the same for the pressure drop, for the fluid temperature at the outle t manifold, and fo r the mean and maximum heater temperature fluctua tions, and all these fluctuations are in phase. At a constant value of mass flux, the oscillation amplitude was found to increase with increase in heat flux. Wu et al., (2006) 8 parallel silic on trapezoidal microchannels Dh = 72.7 m Working fluid: Subcooled water Tin = 30 oC; G = 65.1 2620.7 kg/m2s q = 0 257.2 kW/m2; pout = 1 bar Found the onset of nucleate boiling depends on both the amount of heat flux and mass flux. Observed boiling inst ability flow modes: the liquid/twophase alternating flow and the liquid/tw ophase/vapor alternating flow. Proposed the occurrence of LTAF and LTVAF, as well as their induced largeamplitude/longperi od oscillations of various measurements, are owing to the reversed flow of vapor core because of the limited expansion space in the microchannels. Chang and Pan (2007) 15 parallel silicon rect angular microchannels Dh = 86.3 m Working fluid: Deionized water G = 22 110 kg/m2s; q = 7.68 87.7 kW/m2 Demonstrated significantly differe nt twophase flow patterns under stable or unstable conditions. For the stable cases bubble nucleation, slug flow or annular flows appear sequentially in the flow direction. For the unstable cases, forward or reversed slug/annular flows appear alternatively in every channel. Showed the magnitude of pressure drop oscillations may be used as an index for the appear ance of reversed flow. PAGE 183 183 Table 51. Continued Author Geometry and fl ow conditions Remarks Huh et al. (2007) A rectangular microchannel Dh = 103. 5 m Working fluid: Water G = 170 360 kg/m2s; q = 200 530 kW/m2 Observed the heated wall temperatur e, pressure drop, and mass flux all fluctuated with a long period and large amplitude at very small mass and heat flow rate conditions. These periodic fluctuations exac tly matched the transition of two alternating flow patterns inside th e microchannel: a bubbly/slug flow and elongated slug/semiannular flow. The flow pattern transition inst ability in the single microchannel caused a cyclic behavior of the wall temperature, pressure, and mass flux, and this behavior had a very long period (100200 s) and large amplitude. Muwanga et al. (2007) Two silicon microcha nnel heat configurations: 45 straight parallel channel and crosslinked paths at three locations with 269 m and 283 m Working fluid: Water Tin = 70 and 80 oC; G = 91.0 228.0 kg/m2s q = 4.1 7.9 W/cm2 Flow oscillations are found to be similar in characteristic trends between the two configurations, showing a decreasing frequency with increasing heat flux. The oscillation amplitude are relatively large and identical in frequency for the inlet temperatur e, outlet temperature, inlet pressure, and pressure drop. Wang et al. (2007) 8 parallel silicon trapezoidal microchannels and a single microchannel with a hydraulic diameter of 186 m and a length of 30 mm Working fluid: Subcooled water Tin = 35 oC q = 226.9 497.8 kW/m2 for parallel channels q = 84.5 495.4 kW/m2 for a single channel Found stable and unstable flow boiling regimes existed, depending on the mass flux. Identified two unstable flow boilin g regimes, with longperiod oscillation (more than 1 s) and shor tperiod oscillation (less than 0.1 s) in temperature and pressure. Wang et al. (2008) 8 parallel s ilicon trapezoidal microchannels Dh = 186 m Working fluid: Subcooled water Tin = 35 oC Investigated effects of inlet/outle t configurations on flow boiling instabilities in parallel microchannels. Found amplitudes of temperature an d pressure oscillations in the TypeB connection are much smaller than those in the TypeA connection under the same heat fl ux and mass flux oscillations. Observed nearly steady flow boiling exists in the parallel microchannels with the TypeC connection under the experimental conditions. PAGE 184 184 Table 52. Dimensions of the microcha nnel used in computational simulation (unit: mm) Wcell Wtop Wbot Ww Hcell Hch Hw Lin Lout Lch L1,6 L2,,5 0.708 0.427 0.208 0.141 0.525 0.146 0.379 4.0 4.0 30.0 2.0 7.5 Table 53. Grid density used in computational simulation Inlet plenum Flow channel Outlet plenum 202040 2020500 202080 Table 54. Summary of the amplitude and oscilla tion of inlet pressure under different heat fluxes and mass fluxes Case q, kW/m2 G, kg/m2s Amplitude (PeaktoPeak), kPa Period, ms I 297.8 3523.2 6.5 1.9 II 307.6 711.8 37.32 1.1 III 495.4 679.5 76.70 0.73 Table 55. Summary of stable and unstable flow boiling regime under different heat fluxes and mass fluxes effq kW/m2 297.8 307.6 495.4 G, kg/m2s 3523.2 711.8 679.5 effqG kJ/kg 0.08 0.43 0.73 Remark Stable with small amplitude and short period Unstable with large amplitude and short period Unstable with large amplitude and short period PAGE 185 185 CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS Forced convection flow boiling in a m icrochannel is the subject of this study. Heat transfer phenomena through the liquidvapor interface, nu cleation and bubble dynamics, twophase flow patterns and flow instabilities in the microcha nnel have been theoretic ally and numerically studied. The solutions for the heat transfer ra tes through the liquidvapor interface and pressure drop, twophase flow structures and flow instabilities in the microchannel were obtained by a semianalytical model and using the FLUENT so ftware with user supplied models. Water was used as the working fluid and a circular and a trapezoidal microchannel were used in the simulation. Bubbly, slug, elongated slug, and annular flow pattern we re mainly predicted in most microchannels (Kandlikar 2002a; Xu et al., 200 5; Revellin et al., 2006). The elongated slug bubble flow is the dominant flow pattern in the mi crochannel. In the semianalytical model, fluid flow and heat transfer for an elongated sl ug bubble were studied. Major conclusions and recommendations for future research are provided in this chapter. 6.1 Summary and Conclusions The results of this research can l ead to the following main conclusions: 1. From the developed semianalytical model incorporated with surface tension effects, heat transfer is mainly governed by evaporation at the vaporliquid thin film interface. The length of the evaporation region decreases with increasing heat flux and wall temperature. The results show a qualitatively good agreemen t with those of Qu and Ma (2002) and Wee et al. (2005). 2. The pressure gradient of liquid phase incr eases with increasing heat flux and wall temperature for the elongated bubble model. The numerical results show that the pressure drop of vapor phase is ten times less than th at of liquid phase affected by surface tension and interfacial stress, however, the pressure drop of the vapor phase may not be ignored due to the interfacial stress. The predicted pr essure drops show a reasonable accuracy with differences of 21.26 % to 39.16%, as compared with the experimental results by Qu and Mudawar (2004). The discrepancy may be mainly due to the channel geometry, and inlet and outlet plenum. The proposed model used a circular microchannel rather than a rectangular channel (Qu and Muda war, 2004) and did not consider the pressure loss at the plenums at inlet and outlet. The trend of the pr edicted pressure drop is consistent with the PAGE 186 186 experimental data reported by Qu and Mudawar (2004). The results of this study provide a basic predictive tool to estimate pressure dr op and information of heat transfer phenomena in microchannel heat sinks. 3. In forced convection flow boiling in a micr ochannel, periodic flow boiling with bubbly, churn, and elongated slug flow are investig ated from the numerical simulations. These flow patterns were also observed in the expe rimental results by Wu and Cheng (2003). The transition from bubbly to elonga ted slug flow was found to be the source of flow instabilities. 4. For flow boiling in a microchannel, the numeri cal results show that the wall temperature drops sharply at the nucleation site and the wall temperature becomes mostly constant in the downstream region when nucleate boiling oc curs. The reason is that heat from the heated surface is removed by heat diffusion n ear a bubble due to the large amount of latent heat consumption required for facilitating vapo rization. The wall temperature at nucleation site for the very small channel with 30 m in height is lower than that for the large channels with 1 mm, 700, 500, and 300 m in heights. It is because that the extremely small channel leads to lowest degree of s uperheat requirement for the onset of nucleate boiling (ONB). 5. The heat transfer coefficients at lower qual ities are much greater than those at higher qualities in a microchannel. The results show that the mechanism of nucleate boiling is dominant in a microchannel. A microlayer beneath a bubble may play important role in the heat transfer as well. In a microchannel, the effect of channel height s on the heat transfer coefficient is significant. The h eat transfer coefficient is a strong function of heat flux and quality in a nucleate boiling region, while the ef fect of the quality is nearly as small in saturated boiling region. 6. Wall temperature and inlet pressure oscillations caused by flow instabilities are investigated in a microchannel w ith a hydraulic diameter of 186 m. The numerical results show large amplitude and short period oscill ation for wall temperature and inlet pressure fluctuations, while Wang et al. (2007) showed large amplitude and long period oscillations in their experiments. The discrepancy between the current results and the experimental results by Wang et al. (2007) may be due to co mpressible volume at an upstream region, higher subcooled liquids, and annular/mist flow observed in the experiment. The current study observed only bubbly and el ongated slug flow patterns and used very low subcooled liquid. Those periodic flow pa tterns and lower subcooling lead to short period oscillations. The amplitudes of wall temperature and inlet pre ssure fluctuations are in the range of 12 17.79 C and 6.5 76.70 kPa, respectively. The obtained results point out that the amplitude and oscillation period are strongly affected by th e imposed heat flux and mass flux. 7. Twophase flow boiling instabilities in a mi crochannel were investigated. Stable and unstable cases were identified and they all fi t correctly in the flow boiling stability map developed by Chang and Pan (2007). PAGE 187 187 6.2 Recommendations for Future Research Numerical simulations have been conducted to provide the fundament al understanding of forced convection flow boiling in a microchannel and significan t progress has been achieved. However, due to the complex phenomena of twopha se flow boiling, there are several key issues which remain unsolved and the following are recommended for future investigations: 1. Surface tension effects and nonisothermal interface for the elongated slug bubble have been considered in the semianalytical model. A bulk liquid region has not been investigated in the modeling. It is suggeste d to couple the semianalytical model with the bulk liquid region to consider convecti on effect over the entire channel. 2. The effects of the channel size, heat flux, and mass flux on nucleation and bubbly dynamics have been studied in this work. Various operating pressure effects are not included and needs further study. 3. The current research only considered the single channel case for stability study. Parallel channels have been used in many experimental works, therefore it is useful to conduct numerical simulations for those cases. 4. It took three to four months to finish threedimensional modeling of flow instabilities with only limited results in this research. The numerical simulation of flow instabilities with large subcooling needs more computing power; however, it may provide more results to compare with experimental data. PAGE 188 188 APPENDIX A DERIVATION OF FILM THICKNESS AND E QUILIBRIUM NONEVAPORATING FILM THICKNESS A.1 Liquid Film Thickness The difference of liquid and vapor pressure can be obtained using the modified YoungLaplace equation vlcd p ppp (A1) where v p is the vapor pressure, l p is the liquid pressure, c p is the capillary pressure, and d p is the disjoining pressure. The capil lary and disjoining pressures aff ect the liquid film shape and the heat transfer phenomena in the liquid thin film region. The disjoining pressure, d p for nonpolar liquids is 3dA p (A2) where the Hamaker constant A is determined by the solidliqui d properties and is taken to be positive. For polar liquids such as water with the strong van der Waals force, the disjoining pressure is defined by lnB dlgipRTA (A3) with 1.49A and 0.0243B The gradients of the disjoining pressure for nonpolar and polar liquids are 43ddp Ad dzdz (A4) lnlgi B di lgBRT dpdT d RA dz dzdz (A5) The capillary pressure related with the local meniscus curvature is defined by c p K (A6) PAGE 189 189 where the local curvature,K, is defined by 3/2 1/2 22 2 21 11dd d K dzdzRdz (A7) The local curvature defined in Eq. (A7) in the extended meniscus region can be written 2 21 d K dzR (A8) because the liquid film slope ma y be small. Combining Eqs. (A1), (A6), and (A7), Eq. (A1) can be written 3/2 1/2 22 2 21 11vl ddd d p pp dzdzRdz (A9) Differentiating Eq. (A1) with Eq. (A6) with respect to z, the local curvature grad ient is obtained as 1vl ddpdpdp dK d K dzdzdzdzdz (A10) Differentiating Eq. (A7), the results yields 3/2 5/2 2 22 32 32 1/2 3/2 22 2 21 2131 11 dKdd ddd dzdzdz dzdzdz dd ddd RR dzdz dzdzdz (A11) In the extended meniscus region, differentiating Eq. (A8) gives a relationship as shown: 3 2 31 dKdd dzdz dz R (A12) Combining Eq. (A10) and Eq. (A11), we obtain as follows: PAGE 190 190 3/2 5/2 2 22 32 32 3/2 2 2 1 2 2131 1 ddddd dzdzdzdzdz ddd R dzdzdz d R d 1/2 21 1 0vl dd zdz dpdpdp d K dzdzdzdz (A13) Equation (A13) is multiplied by 3/2 21d dz and arranging th ese results, 1 2 2 32 32 2 2 21 2 3/2 2 231 1 1 1 1vlddddd dzdzdzdz dddd RR dzdz dzdz dpdpdp d dzdzdzdz Kdd dzdz 3/20 (A14) and substituting Eq. (A7) into Eq. (A14), the thirdorder nonlinear differe ntia l equation is obtained as follows: PAGE 191 191 1 2 2 32 32 2 2 21 2 3/2 2 2 231 1 1 1 11 1vlddddd dzdzdzdz dddd RR dzdz dzdz dpdpdp d dzdzdzdz dd dzRdz 20idT d dTdz (A15) In the extended meniscus, the thirdorder nonlinear differential equation is 3 2 2 3 21111vlddpdpdp dddd dz dzdzdzdzdzRdz R (A16) We let d dz (A17) d dz (A18) Combining Eq. (A15) Eq. (A18), the firstorder differen tial equation can be obtained as f ollows: 1 22 122 3/2 2 231 1 1 1 11 1vld id dz RR dpdpdp dzdzdz dT d Rd T d z (A19) 211 11vlddpdpdp dd dz dzdzdzRdz R (A20) PAGE 192 192 A.2 Equilibrium NonEvaporating Liquid Film With use of a kinetic theory to relate the net mass flux crossing a liquidvapor interface to interfacial temperature differential and an extende d Clapeyron equation to e quate the variation in the equilibrium vapor pressure with temperatur e and disjoining pressure the evaporation mass flux escaping from the liquidvapor interface is gi ven as a function of temperature and pressure jump at the interface evapiv lvmaTTbpp (A21) with 1/2 1/22 22 2 22vlv uiuvi lv uiui p Mh M a R TRTT Vp M b RTRT (A22) where is the accommodation coefficient, M is the liquid molecular weight, u R is the universal gas constant, iT is the interfacial temperature, v p is the saturated pressure at vT lvh is the latent of vaporization per unit mass at iT and lV is the liquid molar volume at iT In the equilibrium nonevaporating liquid film the liquidvapor interfacial temperature approaches and becomes equal to the wall temp erature. The equilibrium nonevaporating liquid film region thickness can be obtained by using assumptions, such as 0evapm iwTT 0 and 0 For polar liquids, Equation (A1) can be written as lnB lv lgwppKRTA (A23) where the local curvature K with 0 and 0 is 01 K R PAGE 193 193 With the above mentioned assumptions and Eq. (A23), the equilibrium nonevaporating liquid film thickness 0l nB wvlgwaTTbKbRTA can be obtained by 1/ 1/ 0expB wv B lgwbKaTT A bRT (A24) For nonpolar liquids, substituting Eq. (A2) and Eq. (A6) into Eq. (A1), the resulting equation yields 3lvA ppK (A25) Combining Eq. (A21) and (A25) with 0evapm wTT 0 and 0 the equilibrium nonevaporating liquid film thickness for nonpolar liquids 30wvaTTbKbA is 1/3 0wvbA aTTbK (A26) PAGE 194 194 APPENDIX B DERIVATION OF LIQUID AND INTERFACIAL TEMPERATURE B.1 Liquid Temperature Profile Liquid temperature profile can be deri ved by using the following assumptions: Steady state The energy transfer caused by the convection term is smaller than that resulted from the conduction in the liquid film because the thickness of the liquid film is very small. Thus, the heat transfer in the liquid f ilm results from the conduction only. Thermal diffusion in the axial direction is negligible No heat generation Incompressible flow Viscous dissipation is negligible For constant wall temperature, we get the energy equation as follows: 0lk T r rrr (B1) with boundary conditions are () ()w iTrRT TrRT (B2) Integrating Eq. (B1) twice gives 12()ln TrCrC (B3) Using the boundary conditions, Eq. (B2), we obtain 1 2ln ln lnwi wi iTT C RR TT CT R RR (B4) Substituting Eq. (B4) into Eq. (B3), the temperature profil e in the liquid f ilm is PAGE 195 195 ln lnwi iTT r TrT R RR (B5) The liquid bulk temperature, lT is given by 2 22()R l R lrTrdr Tz RR (B6) Substituting Eq. (B5) in Eq. (B6), the resulting equations are 2 22 ln lnR wi li RTTz r Tz rTz dr Rz RRz RR The integrand in the above equation yields 2 2 2 21 2 11 1 ln 22 4 lni wiTzRR TTz R RR R RR Combining the results and Eq. (B6), the bulk liquid temperature is 2 2 2 22 ln 11 1 ln 22 4wi liTTz TzTz RR RR R RR R (B7) The derivative of lTz with respect to axial direction is 2 2 2 221 ln 11 1 ln 22 4liwidTdTdTdT dzdzdzdz RR RR R RR R (B8) Note that the liquid film thickness change in the axial direction is assumed to be negligibly small. For constant heat flux, boundary conditions are as follows PAGE 196 196 ()lw rR iT kq r TrRT (B9) where wq takes a negative value. Combining Eq. (B9) and Eq. (B3), we get 1 2ln w l w i lq CR k q CTRR k (B10) The liquid temperature profile is ln w i lq r TrTR kR (B11) Combining Eq. (B11) and Eq. (B6), the results are 2 2 2 2 2 2 2 2 222 2222l n 2 lnln 2 2 ln ln ln ln 242422R w i R l l RR w i RR i w lq r rTRdr kR T RR q rTdrRrrrRdr k RR T RR RR RRR q RR R RR R RR k Arranging the results, the liquid bulk temperature is 2 2 2 22 1 ()() ln 224w li lR q RRR TzTz kR RR (B12) The film thickness change in the axial directi on is assumed to be negligibly small, the gradient of the liquid bulk temper ature along the axia l direction is PAGE 197 197 lidTdT dzdz (B13) The wall temperature can be obtained from Eq. (B11) ln w wi lq R TzTzR kR (B14) B.2 Interfacial Temperature The interfacial heat flux q is defined by ilvevapqhm (B15) The interfacial heat flux is approx imately equal to the wall heat flux wq Thus, the interfacial heat flux, 22iwqRqR yields ln/wi ilTT qk RRR (B16) Combing Eq. (B15) with Eq. (B16), the evaporative mass flux is ln/lwi evap lvkTT m h RRR (B17) The interfacial temperature variatio n along the liquidvapor interface is ()ln/lv iwevap lh TTmRRR k (B18) Differentiating Eq. (B18), the gradient of the interfacial temperature is 1lniwlv evap ldTdTh R d m dzdzkRdz (B19) For polar liquid, using the definition of evaporative mass flux, Equation (B17) can be written PAGE 198 198 ln () ln/B lwi iv lgi lvkTT aTTbKRTA hR RR and arranging the results 11 () () ln/ ln/ lnlw li lv lv B vl g ikT kT hR hR RR RR aTaTbKbRTA 1 () ln/ 1 ln () ln/lw v lv B l lg i lvkT aTbK hR RR k ab R A T hR RR We obtain the following equation as: 1 () ln/ 1 ln () ln/lw v lv i B l lg lvkT aTbK hR RR T k ab R A hR RR (B20) The interfacial temperature can be expressed by lnlwv i B llgkTGaTbK T akGbRA (B21) with()ln/fgGhRRR For nonpolar liquids, the interfacial temperature 311 ()() ln/ln/ lw li lv lv ivkT kT hR hR RR RR A aTaTbKb 31 () ln/ 1 () ln/lw v lv l i lvkT A aTbKb hR RR k aT hR RR PAGE 199 199 is obtained by 31 () ln/ 1 () ln/lw v lv i l lvkT A aTbKb hR RR T k a hR RR (B22) Equation (B22) can be written as 3 lw v i lkT A aTbKb G T k a G (B23) For constant heat flux, the interfacial heat flux is iwR qq R (B24) For polar liquids, the interfacial temperature lnB w ivlgi lvqR aTTbKRTA hR yields lnwv lv i B lgR qa TbK Rh T abRA (B25) For nonpolar liquids, the interfacial temperature 3 w iv lvqR A aTTbK hR yields 3 w v lv iqR A bKbaT hR T a (B26) PAGE 200 200 APPENDIX C NUSSELT NUMBER AND HEAT TRANSFE R COE FFICIENT CORRELATION C.1 Nusselt Number Fully developed Nusselt number fo r a threeside heated recta ngular channel (K andlikar, 2005) 2 024 23 1341cc cccAAA Nu AAA (C1) where cab is the channel aspect ratio and a is the unheated side in a threeside heated rectangular channel and 08.2321A 12.0263A 21.2771A 30.29805A 42.2389A and 3 66.532210A. C.2 Heat Transfer Coefficient Correlation For flow boiling in microchannels with Re100lo, modified heat transfer coefficients based on Kandlikar correlations (1990) is defi ned by (Kandlikar and Balasubramanian, 2004) 0.8 0.8 0.2 0.7 ,,,0.0668311058.01cNBD clo FlclohC oxhB oxF h (C2) 0.8 0.8 0.9 0.7 ,,,1.1361667.21cCBD clo FlclohCoxhBoxFh (C3) where the fluid surface parameter FlF is 1.0 for water, the convection number Co and boiling number B o are defined by 0.5 0.81v lx Co x (C4) w lvq Bo Gh (C5) The heat transfer coefficient of single phase flow is given by clo hk hNu D (C6) PAGE 201 201 LIST OF REFERENCES Agostini, B., Fabbri, M., Park, J. E., Wojtan, L., Thom e, J. R., and Michel, B., 2007, State of the Art of High Heat Flux Cooling Technologies, Heat Transfer Eng., 28, pp. 258281. Bahrami, M., Yovanovich, M. M., and Culham, J. R., 2006, Pressure Drop of Fully Developed, Laminar Flow in Rough Microtube s, J. Fluid Engineering, 128, pp. 632637. Baker, O., 1954, Simultaneous Flow of Oil and Gas, Oil Gas J., 53, pp. 183195. Barnea, D., Shoham, O., Taitel, Y., and Dukler, A. E., 1980, Flow Pattern Transition for Horizontal and Inclined Pipe s: Experimental and Comparison with Theory, Int. J. Multiphase Flow, 6, pp. 217225. Basu, N., Warrier, G. R., and Dhir, V. K., 2002, Onset of Nucleate Boiling and Active Nucleation Site Density During Subcooled Flow Boiling, J. Heat Transfer, 124, pp. 717728. Begg, E., Khrustalev, D., and Faghri, A., 1999, Complete Condensation of Forced Convection TwoPhase Flow in a Miniature Tube, J. Heat Transfer, 121, pp. 904915. Bergles, A. E., 1981, Instabilities in TwoPhase Systems, TwoPhase Flow and Heat Transfer in the Power and Process Industries, Hemisphere, Washingt on, DC, pp. 383411. Bergles, A. E., and Dormer, Jr., T., 1969, Subcooled Boiling Pressure Drop with Water at Low Pressure, Int. J. Heat Mass Transfer, 12, pp. 459470. Bergles, A. E., and Rohsenow, W. M., 1964, The Determination of ForcedConvection SurfaceBoiling Heat Transfer, J. Heat Transfer, 86, pp. 365372. Blangetti, F., and Naushahi, M. K., 1980, I nfluence of Mass Transfer on the Momentum Transfer in Condensation and Evaporation Phenomena, Int. J. Heat Mass Transfer, 23, pp. 16941695. Boure, J. A., Bergles, A. E., and Tong, L. S ., 1973, Review of TwoPhase Flow Instability, Nucl. Eng. Des., 25, pp. 165192. Bowers, M. B., and Mudawar, I., 1994a, High Flux Boiling in Low Flow Rate, Low Pressure Drop MiniChannel and MicroChannel Heat Sinks, Int. J. Heat Mass Transfer, 37, pp. 321332. Bowers, M. B., and Mudawar, I., 1994b, TwoP hase Electronic Cooling Using MiniChannel and MicroChannel Heat Sinks: Part 1Design Cr iteria and Heat Diffusion Constraints, J. Electron. Packag., 116, pp. 290297. Bowers, M. B., and Mudawar, I., 1994c, TwoPhase Electronic Coo ling Using MiniChannel and MicroChannel Heat Sinks: Part 2Flow Rate and Pressure Drop Constraints, J. Electron. Packag., 116, pp. 298305. PAGE 202 202 Brackbill, J. U., Kothe, D. B., and Zemach C., 1992, A Continuum Method for Modeling Surface Tension, J. Comput. Phys., 100, pp. 335354. Brauner, N., and MoalemMaron, D., 1992, Iden tification of the Range of Small Diameter Conduits, Regarding TwoPhase Flow Pattern Transition, Int. Comm. Heat Mass Transfer, 19, pp. 2939. Brutin, D., and Tadrist, L., 2004, Pressure Drop a nd Heat Transfer Analysis of Flow boiling in a Minichannel: Influence of the Inlet Condition on TwoPhase Flow Stability, Int. J. Heat Mass Transfer, 47, pp. 23652377. Brutin, D., Topin, F., and Tadrist, L., 2003, Experimental Study of Unsteady Convective Boiling in Heated Minichannels, Int. J. Heat Mass Transfer, 46, pp. 29572965. Carey, V. P., 1992, LiquidVapor PhaseChange Phenom ena: An Introduction to the Thermophysics of Vaporization and Condensation Processes in Heat Transfer Equipment, Hemisphere Publishing Corp., Washington, DC. Celata, G. P., Cumo, M., and Mariani, A., 1997, Experimental Evaluation of the Onset of Subcooled Flow Boiling at High Liquid Velo city and Subcooling, Int. J. Heat Mass Transfer, 40, pp. 28792885. Chang, K. H., and Pan, C., 2007, TwoPhase Flow Instability for Boiling in a Microchannel Heat Sink, Int. J. Heat Mass Transfer, 50, pp. 20782088. Cheng, P., and Wu, H. Y., 2006, Mesoscale and Mi croscale PhaseChange Heat Transfer, Adv. in Heat Transfer, 39, pp. 461563. Chisholm, D., 1983, TwoPhase Flow in Pipelines and Heat Exchangers, Longman Inc., New York. Chung, P. M. Y., Kawaji, M., Kawahara, A., and Shibata, Y., 2004, TwoPhase Flow Through Square and Circular Microcha nnelsEffects of Channel Geometry, J. Heat Transfer, 126, pp. 547552. Cole, R., and Rohsenow, W. M., 1969, Correla tion of Bubble Departure Diameters for Boiling of Saturated Liquids, AIChE Chemical engineering Progress Symposium Series, 65, pp. 211213. Coleman, J. W., and Garimella, S., 1999, Charact erization of TwoPhase Flow Patterns in Small Diameter Round and Rectangular Tube s, Int. J. Heat Mass Transfer, 42, pp. 28692881. Collier, J. G., and Thome, J. R., 1994, Convective Boiling and Condensation, Third ed., Oxford University Press, Oxford. Cornwell, K., and Kew, P. A., 1992, B oiling in Small Parallel Channels, Proceedings of CEC Conference on Energy Efficiency in Process Technology, Athens, October 1992, Paper 22, Elsevier Applied Sciences, pp. 624638. PAGE 203 203 DasGupta, S., Kim, I. Y., and Wayner, P. C., Jr., 1994, Use of the KelvinClapeyron Equation to Model an Evaporating Curved Microfilm, J. Heat Transfer, 116, pp. 10071015. DasGupta, S., Schonberg, J. A., Kim, I. Y., and Wyaner, P. C., Jr., 1993, Use of the Augmented YoungLaplace Equation to Model Equilibrium and Evaporating Extended Menisci, J. Colloid Interface Sci., 157, pp. 332342. Davis, E. J., and Anderson, G. H., 1966, The Incipience of Nucleate Boiling in Forced Convection Flow, AIChE J., 12, pp. 774780. Derjaguin, B. V., Nerpin, S. V., and Churay ev, N. V., 1965, Effect of Film Transfer upon Evaporation of Liquids from Capillaries, Bulletin Rilem, 29, pp. 9398. Derjaguin, B. V., and Zorin, Z. M., 1957, Optical Study of the Absorption and Surface Condensation of Vapors in the Vicinity of Saturation on a Smooth Surface, Proceedings of 2nd International Congr ess on Surface Activity, 2, London, pp. 145152. Dhir, V. K., 1998, Boiling Heat Transfer, Annu. Rev. Fluid Mech., 30, pp. 365401. Diaz, M. C., and Schmidt, J., 2007, Experimental Investigation of Tr ansient Boiling Heat Transfer in Microchannels, Int. J. Heat Fluid Flow, 28, pp. 95102. Dupont, V., and Thome, J. R., 2005, Evaporation in Microchannels: Influence of the Channel diameter on Heat Transfer, Microfluid Nanofluid, 1, pp. 119127. Fedorov, A. G., and Viskanta, R., 2000, ThreeDimensional Conjugate Heat Transfer in the Microchannel Heat Sink for Electronic Packaging, Int. J. Heat Mass Transfer, 43, pp. 399415. Feng, Z. P., and Serizawa, A., 1999, Visualization of TwoPhase Flow Patterns in an UltraSmall Tube, Proceedings of the 18th Multiphase Flow Symposium of Japan, Suita, Osaka, Japan, July 1516, pp. 3336. Feng, Z. P., and Serizawa, A., 2000, Measuremen t of SteamWater Bubbly Flow in UntraSmall Tube, Proceedings of the 37th National Heat Transfer Symposium of Japan, I, pp. 351352. FLUENT 6.2 Users Guide, 2005, Fluent Inc., Lebanon, New Hampshire. Friedel, L., 1979, Improved Friction Pressure Dr op Correlation for Horizont al and Vertical Two Phase Pipe Flow, Paper E2, European Two Phase Flow Group Meeting, Ispra, Italy, pp. 485492. Fukano, T., Kariyasaki, A., and Kagawa, M., 1989, Flow Patterns and pressure drop in isothermal gasliquid flow in a horizontal capillary tube, ANS Proceedings, 1989 National Heat Transfer conference, ISBN 0894481495, ANS, 4, pp. 153161. PAGE 204 204 Garimella, S. V., and Sobhan, C. B., 2003, Tra nsport in Microchannels A Critical Review, Annu. Rev. Heat Transfer, 13, pp. 150. Geisler, K. R. L., 2007, BuoyancyDriven Tw o Phase Flow and Boiling Heat Transfer in Narrow Vertical Channels, Ph.D. Disse rtation, University of Minnesota. Gernet, N. J., Wert, K. L., Baldassarre, G. J., Wilson, M. J., and Wattelet, J. P., Pumped Liquid Cooling System for Desktop Computers, Proceedings IMAPS 38th International Symposium on Microelectronics, Philadelphia, USA, 2005. Ghiaasiaan, S. M., 2008, TwoPhase Flow, Boiling, and C ondensation in Conventional and Miniature Systems, Cambridge University Press, New York, NY. Ghiaasiaan, S. M., and Chedester, R. C., 2002, B oiling Incipience in Microchannels, Int. J. Heat Mass Transfer, 45, pp. 45994606. Gorenflo, D., Knabe, V., and Bieling, V., 1986, Bubble Density on Surfaces With Nucleate Boiling Its Influence on Heat Transfer and Burnout Heat Fluxes at Elevated Saturation Pressures, Proceedings of 8th Internati onal Heat Transfer Conference, Hemisphere, Washington, 4, pp. 19952000. Hapke, I., Boye, H., and Schmidt, J., 2000, Onset of Nucleate Boiling in Minichannels, Int. J. Therm. Sci., 39, pp. 505513. Hetsroni, G., Mosyak, A., Pogrebnyak, E., and Segal, Z., 2005, Explosive Boiling of Water in Parallel MicroChannels, Int. J. Multiphase Flow, 31, pp. 371392. Hetsroni, G., Mosyak, A., Segal, Z., and Ziski nd, G., 2002, A Uniform Temperature Heat Sink for Cooling of Electronic Devices, Int. J. Heat Mass Transfer, 45, pp. 32753286. Hewitt, G. F., 2000, Fluid Mechanics Aspects of TwoPhase Flow, Hand Book of Phase Change Boiling and Condensation, Kandlikar, S. G, Shoji, M., and Dhir, V., eds., Taylor & Francis, Philadelphia. Hino, R., and Ueda, T., 1985, Studies on Heat Tr ansfer and Flow Charact eristics in Subcooled Flow Boiling Part 1. Boiling Characte ristics, Int. J. Multiphase Flow, 11, pp. 269281. Holm, F. W., and Goplen, S. P., 1979, Heat Tr ansfer in the Meniscus ThinFilm Transition Region, J. Heat Transfer, 101, pp. 543547. Hsu, Y. Y., 1962, On the Size Range of Active Nucleation Cavities On a Heating Surface, J. Heat Transfer, 84, pp. 207216. Huh, C., Kim, J., and Kim, M. H., 2007, Flo w Pattern Transition Instability During Flow boiling in a Single Microchannel, Int. J. Heat Mass Transfer, 50, pp. 10491060. Hwang, Y. W., and Kim, M. S., 2006, The Pr essure Drop in Microtubes and the Correlation Development, Int. J. Heat Mass Transfer, 49, pp. 18041812. PAGE 205 205 Jacobi, A. M. and Thome, J. R., 2002, Heat Transfer Model for Evaporation of Elongated Bubble Flows in Microchannels, J. Heat Transfer, 124, pp. 11311136. Jiang, P.X., Fan, M. H., Si, G. S., and Ren, Z. P., 2001, ThermalHydraulic Performance of Small Scale MicroChannel and PorousMedia Heat Exchangers, Int. J. Heat Mass Transfer, 44, pp. 10391051. Jiang, L., Wong, M., and Zohar, Y., 2001, Force d Convection Boiling in a Microchannel Heat Sinks, J. Microelectromech. Syst., 10, pp. 8087. Judy, J., Maynes, D., and Webb, B. W., 2002, Cha racterization of Frictional Pressure Drop for Liquid flows Through Microchannels, Int. J. Heat Mass Transfer, 45, pp. 34773489. Kakac, S., and Bon, B., 2008, A Review of Tw oPhase Flow Dynamic Instabilities in Tube Boiling Systems, Int. J. Heat Mass Transfer, 51, pp. 399433. Kakac, S., Cao, L., and Avelino, M. R., 2003, The Effect of Inlet Subcooling on TwoPhase Flow Dynamic Instabilities in Tube Boiling System, Low Temperature and Cryogenic Refrigeration, Kakac, S., Avelino, M. R., and Sm irnov, H. F., eds., Kluwer Academic Publishers, Dordrecht, Netherlands, pp. 131144. Kakac, S., Veziroglu, T. N., Ozboya, N., and Lee, S. S., 1990, Investigation of Thermal Instabilities in a Forced C onvective Upward Boiling System, Exp. Therm. Fluid Sci., 3, pp. 191201. Kandlikar, S. G., 1990, A General Correlation for TwoPhase Flow Boiling Heat Transfer Inside Horizontal and Vertical Tubes, J. Heat Transfer, 112, pp. 219228. Kandlikar, S. G., 2002a, TwoPhase Flow Pattern s, Pressure Drop, and Heat Transfer During Boiling in Minichannel Flow Passages of Co mpact Evaporators, Heat Transfer Eng., 23, pp. 523. Kandlikar, S. G., 2002b, Fundamental Issues Re lated to Flow Boiling in Minichannels and Microchannels, Exp. Therm. Fluid Sci., 26, pp. 389407. Kandlikar, S. G., 2004, Heat Transfer Mechanis ms During Flow Boiling in Microchannels, J. Heat Transfer, 126, pp. 816. Kandlikar, S. G., 2006, SinglePhase Liquid Flow in Minichannels and Mi crochannels, in Heat Transfer and Fluid Flow in Minichannels and Microchannels, Elsevi er Ltd., Kindlington, Oxford, UK., pp. 87136. Kandlikar, S. G., and Balasubramanian, P., 2004, An Extension of the Flow Boiling Correlation to Transition, Laminar and Deep Laminar Flows in Minichannels and Microchannels, Heat Transfer Eng., 25, pp. 8693. PAGE 206 206 Kandlikar, S. G., and Grande, W. J., 2003, Evolution of Microchannel Flow PassagesThermohydraulic Performance and Fabrica tion Technology, Heat Transfer Eng., 24, pp. 317. Kandlikar, S., G., Mizo, V., Cartwright, M., and Ikenze, E., 1997, Bubble Nucleation and Growth Characteristics in Subc ooled Flow Boiling of Water, National Heat Transfer Conference, HTD342, ASME, pp. 1118. Kandlikar, S. G., Steinke, M. E., Tian, S., Campbell, L. A., 2001, HighSpeed Photographic Observation of Flow Boiling of Wa ter in Parallel Minichannels, Proceedings of National Heat Transfer Conference, ASME, Anaheim, CA, pp. 765684. Kawahara, A., Chung, P. M. Y., and Kawaji, M., 2002, Investigation of TwoPhase Flow Pattern, Void Fraction and Pressure Drop in a Microchannel, Int. J. Multiphase flow, 28, pp. 14111435. Kawaji, M., and Chung, P. M. Y., 2004, Adiabat ic GasLiquid Flow in Microchannels, Microscale Thermophys. Eng., 8, pp. 239257. Kennedy, J. E., Roach, G. M., Jr., Dowling, M. F., AbdelKhalik, S. I., Ghiaasiaan, S. M., Jester, S. M., and Quershi, Z. H., 2000, The Onset of Flow Instability in Uniformly Heated Horizontal Microchannel s, J. Heat Transfer, 122, pp. 118125. Kenning, D. B. R., Wen, D. S., Das, K. S., and Wilson, S. K., 2006, Confined Growth of a Vapour Bubble in a Capillary Tube at Init ially Uniform Superheat: Experiments and Modelling, Int. J. Heat Mass Transfer, 49, pp. 46534671. Kew, P. A., and Cornwell, K., 1997, Correlations for the Prediction of Boiling Heat Transfer in SmallDiameter Channels, Appl. Therm. Eng., 17, pp. 705715. Kocamustafaogullari, G., and Ishi i, M., 1983, Interfacial Area a nd Nucleation Site Density in Boiling Systems, Int. J. Heat Mass Transfer, 26, pp. 13771387. Kosar, A., Kuo, C. J., and Peles, Y., 2005, Suppr ession of Boiling Flow Oscillations in Parallel Microchannels by Inlet Restri ctors, J. Heat Transfer, 128, pp. 251260. Kuo, C. J., Kosar, A., Peles, Y., Virost, S., Mishra, C., and Jensen, M. K., 2006, Bubble Dynamics During Boiling in Enhanced Surf ace Microchannels, J. Microelectromech. Syst., 15, pp. 15141527. Kuznetsov, V. V., and Shamirzaev, A. S., 1999, TwoPhase Flow Pattern and Flow Boiling Heat Transfer in Noncircular Channel with a Small Gap, in TwoPhase Flow Modeling and Experimentation, G.P. Celata, P. Dj Marco and R.K. Shah, eds., Edizioni ETS, pp. 249253. Lee, W. H., 1979, A Pressure Iteration Scheme for TwoPh ase Flow Modeling, Technical Paper No. LAUR79975, Los Alamos, Ne w Mexico, USA: Los Alamos National Laboratory. PAGE 207 207 Lee, P. S., and Garimella, S. V., 2006, Therma lly Developing Flow and Heat Transfer in Rectangular Microchannels of Different Aspect Ratios, Int. J. Heat Mass Transfer, 49, pp. 30603067. Lee, H. J., and Lee. S. Y., 2001a, Pressure Drop Correlations for Tw oPhase Flow within Horizontal Rectangular Cha nnels with Small Heights, Int. J. Multiphase Flow, 27, pp. 783796. Lee, H. J., and Lee, S. Y., 2001b, Heat Tr ansfer Correlation for Boiling Flows in Small Rectangular Horizontal Channe ls With Low Aspect Ratios, Int. J. Multiphase Flow, 27, pp. 20432062. Lee, J., and Mudawar, I., 2005a, TwoPhase Flow in HighHeatFlux Micr oChannel Heat Sink for Refrigeration Cooling Applicat ions: Part I Pressure Drop Ch aracteristics, Int. J. Heat Mass Transfer, 48, pp. 928940. Lee, J., and Mudawar, I., 2005b, TwoPhase Flow in HighHeatFlux Micr oChannel Heat Sink for Refrigeration Cooling Appli cations: Part II Heat Transfer Characteristics, Int. J. Heat Mass Transfer, 48, pp. 941955. Lee, R. C., and Nydahl, J. E., 1989, Numerical Calculation of Bubble Growth in Nucleate Boiling From Inception Through Depa rture, J. Heat Transfer, 111, pp. 474479. Lee, P. C., Tseng, F. G., and Pan, C., 2004, B ubble Dynamics in Microchannels. Part I: Single Microchannel, Int. J. Heat Mass Transfer, 47, pp. 55755589. Li, J., and Cheng, P., 2004, Bubble Cavitation in a Microchannel, Int. J. Heat Mass Transfer, 45, pp. 26892698. Li, H. Y., Tseng, F. G., and Pan, C., 2004, B ubble Dynamics in Microchannels. Part II: Two Parallel Microchannels, In t. J. Heat Mass Transfer, 47, pp. 55915601. Li, J. M., and Wang, B. X., 2003, Seize Effect on TwoPhase Regime for Condensation in Micro/Mini Tubes, Heat TransferAsian Res., 32, pp. 6571. Li, J., and Peterson, G. P., 2005, Microscal e Heterogeneous Boiling on Smooth SurfacesFrom Bubble Nucleation to Bubble Dynamics Int. J. Heat Mass Transfer, 48, pp. 43164332. Li, J., Peterson, G. P., and Cheng, P., 2004, ThreeDimensional Analysis of Heat Transfer in a MicroHeat sink with Single Phase Fl ow, Int. J. Heat Mass Transfer, 47, pp. 42154231. Liu, D., Lee. P. S., and Garimella, S. V., 2005, Prediction of the Onset of Nucleate Boiling in Microchannel Flow, Int. J. Heat Mass Transfer, 48, pp. 51345149. Lockhart, R. W., and Martinelli, R. C., 1949, P roposed Correlation of Data for Isothermal TwoPhase, TwoComponent Flow in Pipes, Chem. Eng. Progress, 45, pp. 3948. PAGE 208 208 Longtin, J. P., Badran, B., and Gerner, F. M., 1994, A OneDimensional Model of a Micro Heat Pipe During SteadyState Oper ation, J. Heat Transfer, 116, pp. 709715. Mala, G. M., and Li, D., 1999, Flow Characteristics of Water in Microtubes, Int. J. Heat Fluid Flow, 20, pp. 142148. Malenkov, I. G., 1971, Detachment Frequency as a Function of Size of Vapor Bubbles, Inzh.Fiz. Zh., 20, pp. 704708. Mandhane, J. M., Gregory, G. A., and Aziz, K ., 1974, Flow Pattern Map for GasLiquid Flow in Horizontal Pipes, In t. J. Multiphase Flow, 1, pp. 537553. Matinelli, R. C., and Nelson, D., 1948, Predictio n of Pressure Drop Duri ng Forced Circulation Boiling of Water, Trans. ASME, 70, pp. 695702. Mei, R., Chen, W. C., and Klausner, J. F., 1995, Vapor Bubble Growth in Heterogeneous BoilingPart I. Formulation, Int. J. Heat Mass Transfer, 38, pp. 909919. Mertz, R., Wein, A., and Groll, M, 1996, Experi mental Investigation of Flow Boiling Heat Transfer in Narrow Channels, Int. J. Heat Technol., 14, pp. 4754. Miller, C. A., 1973, Stability of Moving Surfaces in Fluid Systems with Heat and Mass Transport, AIChE J., 19, pp. 909915. Mirzamoghadam, A., and Catton, I., 1988, A Physic al Model of the Evaporating Meniscus, J. Heat Transfer, 110, pp. 201207. Mishima, K., and Hibiki, T., 1996, Some Char acteristics of AirWater TwoPhase Flow in Small Diameter Vertical Tubes, Int. J. Multiphase Flow, 22, pp. 703712. Moore, F. D., and Mesler, R. B., 1961, The Measurement of Rapid Surface Temperature Fluctuations During Nucleate Boiling of Water, AIChE J., 7, pp. 620624. Moosman, S., and Homsy, G. M., 1980, Evapora ting Menisci of Wetting Fluids, J. Colloid Interface Sci., 73, pp. 212223. Mudawar, I., 2001, Assessment of HighHeatF lux Thermal Management Schemes, IEEE. Trans. Compon. Pack. Technol., 24, pp. 122141. Mukherjee, H., and Brill, J. P., 1985, Empirical Equations to Predict Flow Patterns in TwoPhase Inclined Flow, Int. J. Multiphase Flow, 11, pp. 299315. Mukherjee, A., and Dhir, V. K., 2004, Study of Lateral Merger of Vapor Bubbles During Nucleate Pool Boiling, J. Heat Transfer, 126, pp. 10231039. Mukherjee, A., and Kandlikar, S. G., 2005, Numerical Simulation of Growth of a Vapor Bubble During Flow Boiling of Water in a Mi crochannel, Microfluid. Nanofluid., 1, pp. 137145. PAGE 209 209 Mukherjee, A., and Kandlikar, S. G., 2007, Numer ical Study of Single Bubbles with Dynamic Contact Angle During Nucleate Pool Boiling, Int. J. Heat Mass Transfer, 50, pp. 127138. Muwanga, R., Hassan, I., and MacDonald, R., 2007, Characteristics of Flow Boiling Oscillations in Silicon Microchanne l Heat Sinks, J. Heat Transfer., 129, pp. 13411351. Nayak, A. K., Vijayan, P. K., Jain, V., Saha D., and Sinha, R. K., 2003, Study on the FlowPatternTransition Instability in a Natural Circulation Heavy Water Moderated Boiling Light Water Cooled Reactor, Nucl. Eng. Des., 225, pp. 159172. Park, K., Noh, K. J., and Lee, K. S., 2003, T ransport Phenomena in the ThinFilm Region of a MicroChannel, Int. J. Heat Mass Transfer, 46, pp. 23812388. Patankar, S. V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Company, Washington, DC. Peebles, F, N., and Garber, H. J., 1953, Studi es on Motion of Gas Bubbles in Liquids, Chem. Eng. Prog., 49, pp. 8897. Peles, Y. P., and Haber, S., 2000, A Steady State, One Dimensional, Model for Boiling Two Phase Flow in Triangular MicroChannel, Int. J. Multiphase Flow, 26, pp. 10951115. Peng, X. F., Peterson, G. P., and Wang, B. X., 1994, Heat Transfer Ch aracteristics of Water Flowing Through Microchannels Exp. Heat Transfer, 7, pp. 265283. Philips, R. J., 1990, MicroChannel Heat Sinks, Advances in Thermal Modeling of Electronic Components, BarCohen, A., Kraus, A. D., eds., ASME Press, New York, 2, pp. 109184. Potash, M., Jr., and Wayner, P. C., Jr., 1972, Evaporation from a TwoDimensional Extended Meniscus, Int. J. Heat Mass Transfer, 15, pp. 18511863. Qu, W., and Ma, T., 2002, Effects of The Polarity of Working fluids on VaporLiquid Flow and Heat Transfer Characteristics in a Ca pillary, Microscale Thermophys. Eng., 6, pp. 175190. Qu, W., Mala, G. M., and Dongqing, L, 2000, P ressureDriven Water Flows in Trapezoidal Silicon Microchannels, Int. J. Heat Mass Transfer, 43, pp. 353364. Qu, W., and Mudawar, I., 2002, Prediction and Measurement of Incipient Boiling Heat Flux in MicroChannel Heat Sinks, Int. J. Heat Mass Transfer, 45, pp.39333945. Qu, W., and Mudawar, I., 2003a, Measurement a nd Prediction of Pressure Drop in TwoPhase Microchannel Heat sinks, Int. J. Heat Mass Transfer, 46, pp. 27372753. Qu, W., and Mudawar, I., 2003b, Flow Boiling Heat Transfer in TwoPhase MicroChannel Heat SinksI. Experimental I nvestigation and Assessment of Correlation Methods, Int. J. Heat Mass Transfer, 46, pp. 27552771. PAGE 210 210 Qu, W., and Mudawar, I., 2003c, Flow Boiling H eat Transfer in TwoPhase MicroChannel Heat SinksII. Annular TwoPhase Flow Model, Int. J. Heat Mass Transfer, 46, pp. 27732784. Qu, W., and Mudawar, I., 2004, Transport Phenomena in TwoPhase MicroChannel Heat Sinks, J. Electron. Packag., 126, pp. 213224. Qu, W., Yoon, SM., and Mudawar, I., 2004, T woPhase Flow and Heat Transfer in Rectangular MicroChannels J. Electron. Packag., 126, pp. 288300. Ravigururajan, T. S., 1998, Impact of Channel Geometry on TwoPhase Flow Heat Transfer Characteristics of Refrigerants in Microchannel Heat Exchangers, J. Heat Transfer, 120, pp. 485491. Renk, F. J., and Wayner, P. C., Jr., 1979, An Evaporating Ethanol Menisc us: Part II Analytical Studies, J. Heat Transfer, 101, pp. 5962. Revellin, R., Dupont, V., Ursenbacher, T., Thome, J. R., and Zun, I., 2006, Characterization of Diabatic TwoPhase Flows in Microchannels: Flow Parameter Results for R134a in a 0.5 mm Channel, Int. J. Multiphase Flow, 32, pp. 755744. Revellin, R., and Thome, J. R., 2007, A New Ty pe of Diabatic Flow Pattern Map for Boiling Heat Transfer in Microchannels, J. Micromech. Microeng., 17, pp. 788796. Rider, W. J., and Kothe, D. B., 1998, Reconstructing Volume Tracking, J. Comput. Phys., 141, pp. 112152. Saha, P., and Zuber, N., 1978, An Analytical Study of the Thermally Induced TwoPhase Flow Instabilities Including the Effect of Thermal NonEquilibrium, Int. J. Heat Mass Transfer, 21, pp. 415426. Sartre, V., Zaghdoudi, M. C., and Lallemand, M., 2000, Effect of Interfacial Phenomena on Evaporative Heat Transfer in Micro Heat Pipes, Int. J. Therm. Sci., 39, pp. 498504. Sato, T., and Matsumura, H., 1963, On the Conditions of Incipient Subcooled Boiling and Forced Convection, Bull. Jpn. Sci. Mech. Eng., 7, pp. 392398. Serizawa, A., Feng, Z., and Kawara, Z., 2002, TwoPhase Flow in Microchannels, Exp. Thermal Fluid Sci., 26, pp. 703714. Sobhan, C. B., and Garimella, S. V., 2001, A Comp arative Analysis of Studies on Heat Transfer and Fluid Flow in Microchannels Microscale Thermophys. Eng., 5, pp. 293311. Son, G., and Dhir, V. K., 1997, Numerical Simulation of Saturated Film Boiling on a Horizontal Surface, J. Heat Transfer, 119, pp. 525533. Son, G., and Dhir, V. K., 1998, Numerical Simulation of Film Boiling Near Critical Pressures with a Level Set Method, J. Heat Transfer, 120, pp. 183192. PAGE 211 211 Son, G., Dhir, V. K., and Ramanujapu, N. K., 1999, Dynamics and Heat Transfer Associated with a Single Bubble During Nucleate Boiling on a Horizontal Surface, J. Heat Transfer, 121, pp. 623631. Spedding, P. L., and Spence, D. R., 1993, Flow Regimes in TwoPhase GasLiquid Flow, Int. J. Multiphase Flow, 19, pp. 245280. Steinke, M., and Kandlikar, S. G., 2006, S inglePhase Liquid Friction Factors in Microchannels, Int. J. Therm. Sci., 45, pp. 10731083. Stenning, A. H., 1964, Instabilities in the Flow of a Boiling Liquid, J. Basic Eng. Trans. ASME Ser. D., 86, pp. 213228. Stephan, P., and Busse, C. A., 1992, Analysis of the Heat Transfer Coefficient of Grooved Heat Pipe Evaporator Walls, Int. J. Heat Mass Transfer, 35, pp. 383391. Sujanani, M., and Wayner, P. C., Jr., 1992, Tra nsport Processes And Inte rfacial Phenomena in an Evaporating Meniscus Chem. Eng. Comm., 118, pp. 89110. Suo, M., and Griffith, P., 1964, TwoPhase Flow in Capillary Tubes, J. Basic Eng. 86, pp. 576582. Swanson, L. W., and Herdt, G. C., 1992, Model of the Evaporating Meniscus in a Capillary Tube, J. Heat Transfer, 114, pp. 434441. Tabatabai, A., and Faghri, A., 2001, A New TwoPhase Flow Map and Transition Boundary Accounting for surface Tension Effects in Horiz ontal Miniature and Micro Tubes, J. Heat Transfer, 123, pp. 958 968. Tadrist, L., 2007, Review on TwoPhase Flow Instab ilities in Narrow Spaces, Int. J. Heat Fluid Flow, 28, pp. 5462. Taitel, Y., and Dukler, A. E., 1976, A Model for Predicting Flow Regime Transitions in Horizontal and Near Horizontal GasLiquid Flow, AIChE J., 22, pp. 4755. Thom, J. R. S., Walker, W. M., Fallon, T. A., a nd Reising, G. F. S., 1965, Boiling in Subcooled Water During Flow Up Heated Tubes or Annuli, Symposium on Boiling Heat Transfer in Steam Generating Units and Heat Exchangers, Institute of Mechanical Engineers, New York, pp. 226246. Thome, J. R., 2006, StateoftheArt Overview of Boiling and TwoPhase Flows in Microchannels, Heat Transfer Eng., 27, pp. 419. Thome, J. R., Dupont, V., and Jacobi, A. M., 2004, Heat Transfer Model for Evaporation in Microchannels. Part I: Pr esentation of the Model, Int. J. Heat Mass Transfer, 47, pp. 33753385. PAGE 212 212 Tran, T. N., Chyu, M. C., Wambsganss, M. W., and France, D. M., 2000, TwoPhase Pressure Drop of Refrigerants During Flow Boiling in Small Channels: An Experimental Investigation and Correlation Development, Int. J. Multiphase Flow, 26, pp. 17391754. Triplett, K. A., Ghiaasiaan, S. M., AbdelKhal ik, S. I., Sadowski, D. L., 1999a, GasLiquid TwoPhase in Microchannels Part I: TwoPhas e Flow Patterns, Int. J. Multiphase Flow, 25, pp. 377394. Triplett, K. A., Ghiaasiaan, S. M., AbdelKhalik, S. I., LeMouel, A., and McCord, B. N., 1999b, GasLiquid TwoPhase Flow in Microchannels Part II: Void Fraction and Pressure Drop, Int. J. Multiphase Flow, 25, pp. 395410. Tuckerman, D. B., and Pease, R. F. W., 1981, High Performance Heat Sink for VLSI, IEEE Electron Dev. Let., ELD2, pp. 126129. Wambsganss, M. W., Jendrzejczyk, J. A., and France, D. M., 1991, TwoPhase Flow Patterns and Transitions in Small, Horizontal, Recta ngular Channels, Int. J. Multiphase Flow, 17, pp. 327342. Wang, G., Cheng, P., and Bergles, A. E., 2008, Effects of Inlet/Outlet Configurations on Flow Boiling Instability in Parallel Microchannels, Int. J. Heat Mass Transfer, 51, pp. 22672281. Wang, G., Cheng, P., and Wu, H., 2007, Unsta ble and Stable Flow Boiling in Parallel Microchannels and in a Single Microchannel, Int. J. Heat Mass Transfer, 50, pp. 42974310. Wang, B. X., and Peng, X. F., 1994, Experiment al Investigation of Li quid Forced Convection Heat Transfer Through Microchannel s, Int. J. Heat Mass Transfer, 37, pp. 7382. Wayner, P. C., Jr., 1979, Effect of Thin Film H eat Transfer on Meniscus Profile and Capillary Pressure, AIAA J., 17, pp. 772776. Wayner, P. C., Jr., 1991, The Effect of Inte rfacial Mass Transport on Flow in Thin Liquid Films, Colloids and Surfaces, 52, pp. 7184. Wayner, P. C., Jr., 1994, Thermal and Mechanical Effects in the Spreading of a Liquid Film Due to a Change in the Apparent Fin ite Contact Angle, J. Heat Transfer, 116, pp. 938945. Wayner, P. C., Jr., 1997, Constrained Va por Bubble, Microscale Thermophys. Eng., 1, pp. 111118. Wayner, P. C., Jr., 1999, Intermolecular Forces in PhaseChange Heat Transfer: 1998 Kern Award Review, AIChE J., 45, pp. 20552068. Wayner, P. C., Jr., Kao, Y. K., and Lacroix, L. Y., 1976, The Interline HeatTransfer Coefficient of an Evaporating Wetting Film, Int. J. Heat Mass Transfer, 19, pp. 487492. PAGE 213 213 Wee, S. K., Kihm, K. D., and Hallinan, K. P., 2005, Effects of the Liquid Polarity and the Wall Slip on the Heat and Mass Transport Charac teristics of the MicroScale Evaporating Transition Film, Int. J. Heat Mass Transfer, 48, pp. 265278. Weisman, J., Duncan, D., Gibson, J., and Crawford, T., 1979, Effects of Liquid Properties and Pipe Diameter on TwoPhase Flow Patterns in Horizontal Lines, Int. J. Multiphase Flow, 5, pp. 437462. Weisman, J., and Kang, S. Y., 1981, Flow Pattern Transitions in Vertical and Upwardly Inclined Lines, Int. J. Multiphase Flow, 7, pp. 271291. Welch, S. W., 1998, Direct Simulation of Vapor Bubble Growth, Int. J. Heat Mass Transfer, 41, pp. 16551666. Welch, S. W., and Wilson, J., 2000, A Volume of Fluid Based Method for Fluid Flows with Phase Change, J. Comput. Phys., 160, pp. 662682. Wu, H. Y., and Cheng, P., 2003a, Fricti on Factors in Smooth Trapezoidal Silicon Microchannels with Different Aspect Ra tios, Int. J. Heat Mass Transfer, 46, pp. 25192525. Wu, H. Y., and Cheng, P., 2003b, An Experime ntal Study of Convective Heat Transfer in Silicon Microchannels with Different Surface Conditions, Int. J. Heat Mass Transfer, 46. pp. 25472556. Wu, H. Y., and Cheng, P., 2003c, Visualization and Measurements of Periodic Boiling in Silicon Microchannels, Int. J. Heat Mass Transfer, 46, pp. 26032614. Wu, H. Y., and Cheng, P., 2004, Boiling Instab ility in Parallel Silicon Microchannels at Different Heat Flux, Int. J. Heat Mass Transfer, 47, pp. 36313641. Wu, H. Y., Cheng, P., and Wang, H., 2006, P ressure Drop and Flow Boiling Instabilities in Silicon Microchannel Heat Sinks, J. Micromech. Microeng., 16, pp. 21382146. Xu, J., Gan, Y., Zhang, D., and Xiuhan, L., 2005, Microscale Boiling Heat Transfer in a MicroTimescale at High Heat Fluxes, J. Micormech. Microeng., 15, pp. 362376. Xu, B., Ooi, K. T., Wong, N. T., and Choi, W. K., 2000, Experimental Investigation of Flow Friction for Liquid Flow in Microchannels, Int. Co mm. Heat Mass Transfer, 27, pp. 11651176. Xu, J., Zhou, J., and Gan, Y., 2005, Static and Dynamic Flow Instability of a Parallel Microchannel Heat Sink at High He at Fluxes, Energy Convers. Manage., 46, pp. 313334. Yadigaroglu, G., 1981, TwoPhase Flow Instabilities and Propagation Phenomena, Thermohydraulics of TwoPhase Systems for Industrial Design and Nuclear Engineering, Delhaye, J. M., Giot, M., and Riethmuller, M. L., eds., Hemisphere, Washington DC, pp. 353403. PAGE 214 214 Ye, T., Shyy, W., and Chung, J. N., 2001, A FixedGrid, SharpInterface Method for Bubble Dynamics and Phase Change, J. Comput. Phys., 174, pp. 781815. Yu, W., France, D. M., Wambsganss, M. W., and Hu ll, J. R., 2002, TwoPhase Pressure Drop, Boiling Heat Transfer, and Criti cal Heat Flux to Water in SmallDiameter Horizontal Tube, Int. J. Multiphase Flow, 28, pp. 927941. Zhang, L, Koo, J. M., Jiang, L., Ashegi, M., Goods on, K. E., Santiago, J. G., and Kenny, T. W., 2002, Measurement and Modeling of TwoPhas e Flow in Microcha nnels with Nearly Constant Heat Flux Boundary Conditi ons, J. Microelectromech. Syst., 11, pp. 1219. Zhang, L, Wang, E. N., Goodson, K. E., and Ke nny, T. W., 2005, Phase Change Phenomena in Silicon Microchannels, Int. J. Heat Mass Transfer, 48, pp. 15721582. Zhao, T. S., and Bi, Q. C., 2001, CoCurrent Air Water TwoPhase Flow Patterns in Vertical Triangular Microchannels, Int. J. Multiphase Flow, 27, pp. 765782. Zivi, S. M., 1964, Estimation of SteadyState Stem VoidFraction by means of the Principle of Minimum Entropy Production, J. Heat Transfer, 86, pp. 247252. Zuber, N., 1963, Nucleate Boiling The Region of Isolated Bubbles Similarity With Natural Convection, Int. J. Heat Mass Transfer, 6, pp. 5365. PAGE 215 215 BIOGRAPHICAL SKETCH Yun W han Na was born in Seoul, South Korea. He grew up and sp ent most of his early life in the city. In 1995, he graduated with a bachel ors degree in mechan ical engineering from Hongik University, Seoul, South Korea. He c ontinued his study at Hongik University and received a Master of Science degree in mechan ical engineering in 1997. With his friends, he founded POSLINK, an IT venture company ba sed in Seoul, South Korea, in 1998. With enthusiasm for and the challenge of the IT ve nture, he managed the company until 2002. In August 2003, he enrolled in the Department of Mechanical and Aerospace Engineering at the University of Florida. While enrolling at UF, he has worked on the topic of this dissertation under the guidance of Dr. Chung and obtained a Doctor of Philosophy degree in fall 2008. 