UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 A FLOW BOILING MICROCHANNEL THERMO SYPHON FOR FUEL CELL THERMAL MANAGE MENT By PATRICK THOMAS GARRITY 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 2009 1 PAGE 2 2009 Patrick Thom as Garrity 2 PAGE 3 To m y parents, James and Nancy Garrity. Wit hout their continuous support and encouragement this work would not have been possible. 3 PAGE 4 ACKNOWL EDGMENTS The author would like to express his sinceres t gratitude to his academic advisor and PhD committee chairperson, Professor James Fred erick Klausner. His guidance, support, encouragement and insight contributed immensel y to the accomplishment of this work. The author would also like to express his appreciati on to Professor Renwei Mei for his invaluable input throughout the course of my studies. Th e author would like to thank Professor Jacob Chung, Professor Mark Sheplak and Professo r Ranga Narayanan for serving on his PhD committee. Their insight has helped the author to develop the necessary critical thinking skills needed to become a contributing member of the academic community. This research was supported by the Nati onal Aeronautics and Space Administration Glenn Research Center, through contract NA632750. Without its financial support this work would not be possible. 4 PAGE 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF TABLES................................................................................................................. ..........8 LIST OF FIGURES.........................................................................................................................9 LIST OF ABBREVIATIONS........................................................................................................12 ABSTRACT...................................................................................................................................17 CHAPTER 1 INTRODUCTION................................................................................................................. .17 PEM Fuel Cells................................................................................................................. ......18 TwoPhase Flow and Thermosyphons....................................................................................22 Microchannels.........................................................................................................................25 Compact Condensers............................................................................................................. .30 Instabilities and Subcooling................................................................................................... .35 Summary.................................................................................................................................37 2 FLOW BOILING MICROCHANNEL EVAPORATOR AND THERMOSYPHON PERFORMANCE...................................................................................................................3 9 Background and Objectives....................................................................................................39 Experimental Facility.......................................................................................................... ....42 Experimental Protocol............................................................................................................46 Results.....................................................................................................................................48 Flow Modeling........................................................................................................................51 3 MICROCHANNEL PRESSURE DROP CORRELATION..................................................59 Background and Objectives....................................................................................................59 Experimental Facility.......................................................................................................... ....59 Experimental Protocol............................................................................................................60 Results.....................................................................................................................................61 Pressure Drop Correlations.....................................................................................................64 Conclusion..............................................................................................................................69 5 PAGE 6 4 INSTABILITY PHENOMENA IN A TWOPHASE THERMOSYPHON..........................70 Background and Objectives....................................................................................................70 Experimental Considerations..................................................................................................71 Facility....................................................................................................................... ......71 Experimental Protocol.....................................................................................................72 System Subcooling and Manifold Pressure Drop............................................................72 Flow Modeling........................................................................................................................74 Experimental Results..............................................................................................................76 Instability Prediction......................................................................................................... ......82 5 PERFORMANCE OF ALUMINUM AND CA RBON FOAMS FOR AIR SIDE HEAT TRANSFER AUGMENTATION..........................................................................................86 Background and Objectives....................................................................................................86 Experimental Facility.......................................................................................................... ....88 Experimental Results..............................................................................................................91 Heat Transfer Analysis...........................................................................................................95 Results...................................................................................................................................101 Hypothetical Heat Exchanger Performance..........................................................................103 Numerical Model................................................................................................................ ..105 6 CONCLUSIONS.................................................................................................................. 111 APPENDIX A INSTABILITY IDENTIFICATION....................................................................................113 B UNCERTAINTY AND CALIBRATION............................................................................116 Data Acquisition (DAQ).......................................................................................................11 6 Pressure Drop.................................................................................................................. ......117 Flow Meter Calibration......................................................................................................... 120 Uncertainty in Power Measurement.....................................................................................122 Temperature.................................................................................................................... ......123 Various Parameters............................................................................................................. ..125 Wall Temperature..........................................................................................................125 Vapor Quality................................................................................................................125 Vapor Superficial Velocity............................................................................................125 Pitot Tube Velocity........................................................................................................126 C PRESSURE GRADIENT AND HEAT TR ANSFER DATA FOR FOAM SAMPLES......127 D NUMERICAL CODE FOR SOLVING POROUS MEDIA TRANSPORT EQUATIONS.......................................................................................................................130 LIST OF REFERENCES.............................................................................................................145 6 PAGE 7 BIOGRAPHICAL SKETCH .......................................................................................................152 7 PAGE 8 LIST OF TABLES Table page 11 Typical values of the convec tion heat transfer coefficient................................................24 12 Various types of twophase flow instability......................................................................37 31 Frictional pressure drop data.............................................................................................. 62 32 Error between measured and predicted frictional pressure drop for various models........69 41 Experimentally observed limiting heat flux.......................................................................78 51 Foam properties............................................................................................................ .....92 52 Relative error for pressure dr op and Nusselt number correlations..................................103 53 Geometric parameters for louvered fin configuration.....................................................107 B1 Diaphragm maximum pressure, sensitivity parameter, and standard error of the estimate for each pressure transducer..............................................................................120 C1 Pressure gradient for different foams...............................................................................128 C2 Volumetric heat transfer coefficient and upper fall temp erature for different foams......129 8 PAGE 9 LIST OF FI GURES Figure page 11 Operation of a PEM fuel cell.............................................................................................19 12 Thermosyphon............................................................................................................... ....23 13 Microchannel and conventional channel configurations...................................................26 14 Crossflow condenser........................................................................................................31 15 Thermal transport mechanisms inside the condenser........................................................32 16 Thermal resistance inside the condenser............................................................................32 17 Flow field and boundary la yer for louver type fins...........................................................34 21 Experimental twophase thermosyphon.............................................................................42 22 Exploded view of the cooling plate assembly....................................................................43 23 Section view of flow meter................................................................................................4 4 24 Thermocouple locatio ns on cooling plate..........................................................................45 25 Degree of inlet sub cooling variation with heat flux..........................................................47 26 Variation of flow rate with heat flux..................................................................................48 27 Wall temperature at various locations................................................................................49 28 Variation in vapor quali ty with channel height..................................................................51 29 Measured and predicted mi crochannel pressure drop va riation with flow rate.................54 210 Comparison between measured and predicted flow rate...................................................55 211 Comparison between measured and pred icted wall temperature at thermocouple location 5 based on measured flow rate and vapor quality................................................56 31 Thermosyphon experimental facility.................................................................................60 32 Experimental versus predicted frictiona l pressure drop using Friedal correlation............67 33 Experimental versus predicted frictio nal pressure drop using MuellerSteinhagen correlation..........................................................................................................................67 34 Experimental versus predicted frictiona l pressure drop using homogeneous model.........68 9 PAGE 10 35 Experimental versus predicted frictio nal pressure drop using Lee correlation..................68 41 Variation of inlet s ubcooling with heat flux......................................................................73 42 Pressure drop across inlet and outlet manifolds.................................................................74 43 Mass flow rate variation with increasing heat flux............................................................76 44 Pressure drop variation across the micr ochannel evaporator plate with respect to mass flow rate................................................................................................................. ...77 45 Upper and lower bounds for predicted mass flow rates for H = 1.33 m and H = 0.79 m........................................................................................................................................78 46 Measured probability density function for pressure drop prior to unstable flow..............79 47 Measured probability density function fo r mass flow rate prior to unstable flow.............79 48 Measured probability density function for mass flow rate following the onset of unstable flow......................................................................................................................80 49 Measured probability density function for pressure drop following the onset of unstable flow......................................................................................................................81 410 Power spectrum for mass flow rate....................................................................................82 411 Pressure drop across the microchannel pl ate for varying flow rate at H=1.33 m..............84 412 Comparison of measured and predicted heat flux at the onset of instability.....................85 51 Foam samples............................................................................................................... ......89 52 Porous media experimental facility...................................................................................91 53 Axial pressure upper wall temperature va riation for 10 PPI aluminum foam sample, = 1.13mu ms.....................................................................................................................93 54 Carbon foam pressure gradient va riation with mean air velocity......................................93 55 Aluminum foam pressure gradient variation with mean air velocity................................94 56 Foam upper wall temperature variation with mean air velocity........................................94 57 Dimensionless velocity profile..........................................................................................97 58 Dimensionless temperature profile....................................................................................98 59 Carbon foam volumetric heat transfer coeffi cient variation with mean fluid velocity......99 10 PAGE 11 510 Aluminum foam volumetric heat transfer coefficient variation with mean fluid velocity.............................................................................................................................100 511 Nusselt number correlation for aluminum foams............................................................102 512 Nusselt number correlation for carbon foams..................................................................103 513 Hypothetical heat ex changer in cross flow......................................................................104 514 Geometric comparison of louvered fin a nd foam heat exchange r configurations...........107 515 Comparison of coefficient of performance for louvered fin and foam configurations....108 516 Comparison of compactness factor for louvered fin and foam configurations................109 517 Comparison of power density for louv ered fin and foam configurations........................110 B1 Pressure drop calibration curve used for measurement of pressure drop........................119 B2 Flow rate versus pressure drop calibration curve used for measurement of flow rate within the venturi flow meter...........................................................................................121 B3 Temperature calibration curv e for a typeE thermocouple..............................................124 11 PAGE 12 LIST OF ABBRE VIATIONS A Area 2ma Surface area per unit volume 1m Bi Biot number C Thermal conductivity ratio F C Inertia coefficient C0 Distribution parameter Cp Constant pressure specific heat (J/kgK) pfC Constant pressure specific heat of fluid J kgK CF Compactness factor 3Wm COP Coefficient of performance d Hydraulic diameter (m) D Machined passage diameter (m) Da Darcy number kD Dimensionless parameter de scribing small scale friction input E Electrical Energy (W) h F Fin height (mm) p F Fin pitch (mm) f Friction factor s f Sampling frequency 1 s G Mass velocity (kg/m2s) g Gravitational constant (m/s2) 12 PAGE 13 H Condenser Height (m) s fh Interstitial heat transfer coefficient 2WmK vh Volumetric heat transfer coefficient 3WmK hfg Latent heat of vaporization (J/kg) K Permeability 2mk Thermal conductivity WmK L Axial length of cartridge heater m aL Louver angle (degrees) dL Louver depth (mm) H L Louver height (mm) pL Louver Pitch (mm) m Mass (kg) m Mass flow rate (kg/s) or (g/s) N Number of samples K Nu Nusselt number P Pressure (Pa) or (kPa) PD Power Density Wkg Pr Prandtl number q Heat rate (W) removedQ Heat removed from hypothetical condenser (W) wq Heat flux along heater wall 2Wm 13 PAGE 14 R Flow resistance Paskg Re K Reynolds number based on permeability Re L Reynolds number based on L Re Reynolds number St Stanton number T Temperature ( C) U Normalized fluid Uvelocity u Fluid uvelocity ms mu Mean fluid velocity ms vu Superficial vapor velocity (m/s) V Volume 3mVvj Drift velocity (m/s) W Width m x Vapor quality x Dimensionless xcoordinate y Dimensionless ycoordinate z Axial location (m) Greek Symbols Vapor volume fraction Porosity p Carbon foam porosity based on interstitial pores f r 2phase friction multiplier 14 PAGE 15 Dynamic viscosity 2Nsm Kinematic viscosity (m2/s) Inclination angle (radians) Density (kg/m3) Surface tension (N/m) P Pressure drop Pa s P Total system pressure drop (Pa) mP Measured pressure drop (Pa) pP Predicted pressure drop (Pa) Subscripts 2 Twophase mixture a Accelerational component amb Ambient air ,bf Bulk fluid e Exit f Frictional component fl Fluid fe Effective fluid property g Gravitational component i Inlet K Permeability based length scale l Liquid 15 PAGE 16 m Mixture mf Mean fluid velocity s Solid s e Effective solid property T Total v Vapor w Wall z zlocation 16 PAGE 17 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 A FLOW BOILING MICROCHANNEL THERMO SYPHON FOR FUEL CELL THERMAL MANAGEMENT By Patrick T. Garrity May 2009 Chair: James Frederick Klausner Major: Mechanical and Aerospace Engineering To provide a high power density thermal management system for proton exchange membrane (PEM) fuel cell applications, a passi vely driven thermal management system was assembled to operate in a closed loop twophase thermosyphon. The system has two major components; a microchannel evaporator plate and a condenser. The microchannel evaporator plate was fabricated with 56 squa re channels that have a 1 mm x 1 mm cross section and are 115 mm long. Experiments were conducted with a liqui d cooled condenser with heat flux as the control variable. Measurements of mass flow ra te, temperature field, a nd pressure drop have been made for the thermosyphon loop. A m odel is developed to predict the system characteristics such as the temperature and pressure fields, flow rate, flow regime, heat transfer coefficient, and maximum heat flux. When the system is subjected to a heat load that exceeds the maximum heat flux, an unstable flow regime is observed that cause s flow reversal and eventual dryout near the eva porator plate wall. This undesira ble phenomenon is modeled based on a quasisteady state assumption, and the model is capable of predicting the heat flux at the onset of instability for quasisteady twophase flow. Another focus of this work is the performa nce of the condenser portion of the loop, which will be air cooled in practice. The aim is to reduce air side thermal resistance and increase the 17 PAGE 18 18 condenser performance, which is accomplished with extended surfaces. A testing facility is assembled to observe the air side heat transfer performance of three aluminum foam samples and three modified carbon foam samples, used as extended surfaces. The aluminum foam samples have a bulk density of 216 kil ograms per cubic meter with pore sizes of 0.5, 1, and 2 mm. The modified carbon foam samples have bulk dens ities of 284, 317, and 400 kilograms per cubic meter and machined flow passages of 3.2 mm. in diameter. Each sample is observed under forced convection with air velocity as the control variable. Ther mocouples and pressure taps are distributed axially along the test section and measurements of pressure and temperature are recorded for air velocities ra nging from 16 meters per second Using the DarcyForcheimer equation, the porosity is determined for each sample The volumetric heat tr ansfer coefficient is extracted by means of solvi ng the coupled energy equations of both the solid and fluid respectively. Nusselt number is correlated with Reynolds number. The optimal foam configuration is explored base d on a Coefficient of Performance, (COP), Compactness Factor (CF) and Power Density (PD). The COP is the ra tio of total heat removed to electrical heat consumption of the blower, CF is the total heat removed per unit volume, and PD is the total heat removed per unit mass. These performance parameters are computed for a hypothetical heat exchanger using each foam sample at various flui d velocities. They are also compared against those for the hypothetical heat exchanger fitted with conventional louvered fins. Given a proper weighting function based on the importance of CF COP, and PD in the condenser design, an optimal configuration for an air cooled conde nser can be obtained for various operating conditions. PAGE 19 CHAP TER 1 INTRODUCTION Fuel cells are expected to play a major ro le in energy production within the foreseeable future. Increasing concerns about pollution a nd possible anthropogenic global warming, coupled with economic issues involving fossil fuels have accelerated research inte rest in alternative energy systems. Energy Systems that run on a clean, economical energy source such as hydrogen are of particular interest It is anticipated that the development and deployment of economical and reliable fuel cel ls could be the catalyst for hydrogen production and usage. Proton Exchange Membrane Fuel Cells (PEMFC) ar e of significant interest to the automobile, avionics and space industries, due to the potentia l for high power density, relatively quick start up, rapid response to varying loads, and low operating temperatures. In addition, there is effectively no pollution produced during operation of PEM fuel cells as water and heat essentially make up all of the exhaust to the environment (Vishnyakov, 2006). Thermal issues arise within a PEM fuel cell stack that can significantly limit its cells performance. By rapidly removing heat from the fu el cell stack as it is generated, it is possible to maintain the stack temperature within its optim al range. Such a heat removal technology can increase both the power density a nd energy density. The objective of this work is to understand and analyze the heat transfer pe rformance of a twophase microcha nnel evaporator plate that can be inserted within a PEMFC stack used in conjunction with a twophase thermosyphon. The system operates in a natural circulation mode which reduces the required pumping energy. Initial experiments are conducted with a liquid cooled condenser, with focus being on the evaporator plate. After taking extensive measur ements of the temperature and pressure fields within the microchannel evaporat or plate, a 1D model is deve loped to predict the thermal hydraulics of the cooling plate including the maximum heat removal capabilities, where the 17 PAGE 20 lim iting heat flux is governed by system instabilities. The ideal working fluid for the thermosyphon loop is chosen with a saturation temper ature that allows for the fuel cell stack to operate within its optimal range. Once the maximum heat removal of the cooling plate is known, along with the properties of the ideal working fluid, focus is shifte d to the condenser design, with a focus on maximizing air side heat transfer rates. Carbon and me tallic foams are configured as extended surfaces for heat transfer augmenta tion. The various metallic and carbon foam materials are placed in a forced convection arrangement and the volumetric heat transfer characteristics of each sample are analyzed and co rrelated. Each foam configuration is evaluated based on a Coefficient of Performance, (COP), Compactness Factor (C F) and Power Density (PD). These performance parameters are comput ed for a hypothetical heat exchanger using each foam sample at various fluid velocities. They are also compared against those for the hypothetical heat exchanger fitted with conventional louvered fi ns. Given a proper weighting function based on the importance of CF, COP, and PD in the condenser design, an optimal configuration for an air cooled condenser can be obtained fo r various operating conditions. PEM Fuel Cells The basic operation of a PEMFC converts hydro gen fuel to electricity by passing the fuel across a solid polymer membrane, a thin plastic film that is permeable to protons, but does not conduct electrons. Thus electrons are forced through an extern al electric circ uit for power production, see Figure 11. The solid polymer me mbrane, also known as the electrolyte, is located between the anode and the cathode side of the cell. At the anode, the hydrogen makes contact with a thin layer of platinum catalyst and the electr on is stripped away from the hydrogen. The liberated electrons c ontinue to flow through the exte rnal circuit while the protons travel through the electrolyte to the cathode side (Barbir, 2005). Th e anode side half reaction is 18 PAGE 21 222 HHe (11) On the cathode side, protons passing thr ough the membrane encounter oxygen and the electrons from the external circuit. They recomb ine at a catalyst electr ode layer to form liquid water, which is drained away so as to not block the fuel cell. The cathod e side half reaction is given in Eq. (12), while the total re action is provided in Eq. (13). 2442 HeOH2 O (12) 2221 2 HOH O (13) Figure 11. Operation of a PEM fuel cell. Prior to operation, in an open circuit, there exists a voltage potential between the anode and cathode sides of the cell due to the free elec trons generated at the electrolyte. Once the circuit is closed and the load is connected, external current is generated initiating the chemical reaction within the cell. As the external resistance is reduced the react ion rates increase allowing 19 PAGE 22 the cur rent to adjust, while simultaneously generati ng different types of inte rnal resistance in the process. In the case of PEM fuel cells, the polyme r membrane is solid; reducing corrosion and electrolyte management problems. However, in order for the hydrogen ions to conduct through to the cathode, sufficient hydration is required. The proton has to first bind with water to form hydronium ions, which are then drawn across the me mbrane by the existing electric field. This implies that every proton moving across the memb rane carries along with it a certain number of water molecules. This phenomenon is referred to as electroosmotic drag (Coppo et al., 2005). If there is not enough water present, a high resistance to proton c onduction takes place while on the other hand, if too much water is present, flooding will occur and block the transport of reactants. This process makes water management essential. Likewise, thermal problems can arise within the stacks that significantly lo wer the fuel cell efficiency. With the hopes of optimizing efficiency, PEM fu el cells have recently been subject to an everincreasing number of modeling efforts for th e past 15 years. Accurate models describing the complex coupling of the heat, mass and moment um transfer within the stack are a crucial step in optimizing fuel cell efficiency with regards to thermal management. Coppo et al. (2005) developed a 3D computational model to describe the effects of temperature on the operation of PEM fuel cells. It is shown that temperature affects the fuel cell differe ntly depending on current density. At low current density, (0 PAGE 23 overpoten tial, which increases current density. However this simultaneously raises electrical resistance, reducing the cell power ou tput. At slightly higher curre nt densities, classified as the ohmic region, electroosmotic drag becomes more significant causing the optimal temperature to be governed by membrane ionic conductivity and dissolved water diffusivity. Higher temperature decreases the ioni c conduction due to the reduce d hydration in the presence of electroosmotic drag. On the other hand, loweri ng the temperature decreases the dissolved water diffusivity. This causes the membrane hydratio n level to be less homogeneous, bringing about localized regions of low ionic conductivity. Finally, for high current density operation, also known as the mass transport limited regime, the optimal temperature is governed by the resistance to reactant transport that results fr om the high mass flow rates required to sustain the electrochemical reactions at a rapid rate (Coppo et al., 2005). By increasing the operating temperature, the oxygen diffusivity within the gas diffusion layer as well as the ionconducting polymer is enhanced. Simultaneously, the air/ water concentration rati o influencing reaction kinetics decreases with temperat ure. The higher water content in the gas stream reduces oxygen concentration at the cell inlet. So even though the oxygen diffu sivity is higher, less oxygen is available when the temperature is too high, reducing the cathode overpot ential and consequently cell voltage. Similarly, water diffusivity in creases with temperature due to the strong temperature dependency of kinematic viscosity, contact angle, and surf ace tension. It should also be noted that there is a strong coupling between the liquid water and oxygen diffusion, since higher liquid water diffusion increases the water c ontent in the gas stream The resulting optimal temperature within the mass tran sport limited regime is governed by a balance of the coupled oxygen and liquid water diffusivities as well as the resulting saturation level of the oxygen at the 21 PAGE 24 cell inle t. It can than be c oncluded that optimal operating temperatures vary depending on the fuel cell output or current density. You and Liu (2002) have developed a twopha se flow transport model for the cathode side of PEM fuel cells. They ha ve demonstrated that there exis ts an operating temperature for optimum performance that tends to be in the range of 7080 C. This low operating temperature, coupled with the large heat flux removal demand renders thermal management a difficult task and further elucidates the need for enhanced mode s of thermal transport. There are also several advantages to the low operating temperature that allow PEM fuel cells to be used as lightweight and mobile energy sources. For instance, at low ope rating temperature, the fuel cell is capable of warming up quickly which can be beneficial for m obile applications. Other advantages include the elimination of expensive containment stru ctures and less wear on system components, resulting in better durability. By implementing va rious heat removal methods such as twophase heat transfer in a thermosyphon, this low operati ng temperature can make PEM fuel cells more practical to energy production. TwoPhase Flow and Thermosyphons The concept of passively driven thermal ma nagement has been the focus of intense research in a wide variety of applications both small and large scale. Twophase thermosyphons are considered to be high efficiency heat transf er devices capable of tran sporting relatively large amounts of heat with a relati vely small temperature drop. This low thermal resistance encountered in thermosyphons has proven useful in small applications such as microelectronics cooling and solar water heaters, as well as larger scale applications such as nuclear reactors and heat exchangers used in the pe troleum industry. Faghr i (1995) and Peterson et al. (1993), among others, have reviewed the theory and the applications of the two phase thermosyphon 22 PAGE 25 technology. Aside from the low thermal resistance and the uniform temperature distribution in both the evaporator and condens er sections of the thermosyphon loop, thermosyphons have many advantages over other cooling techniques. Figure 12. Thermosyphon. In a twophase thermosyphon, see Figure 12, heat is conducted from the heat source, across an interface, to an evaporator plate. At this point the cool ant inside the plate is vaporized, taking with it the latent heat of vaporization. The vapor than travels up the riser and into the condenser where it condenses, tran sferring the latent heat to the t ube wall. Forced air blown over the condenser dissipates heat to the ambient, while the condensate travels back down to the evaporator plate and the cycle is repeated. The density difference between the vapor and the liquid portions of the lo op, create a pressure head that driv es the flow. The thermal resistance encountered within the evaporator plate and throughout the condenser will play a critical role in 23 PAGE 26 the over all system heat removal capabilities. This process has several ad vantages over single phase, forced convection techniques. By utilizing phase change heat transfer, heat rejection capabilities become an order of ma gnitude greater than that for singl ephase, (Table 11). On the other hand, there exists a limiting heat flux at which poi nt the system flow becomes unstable. Table 11. Typical values of the conv ection heat transfer coefficient. Process 2hWmK Free convection Gases 225 Liquids 50100 Forced convection Gases 25250 Liquids 10020,000 Convection with phase change Boiling or Condensation 2500100,000 (Incropera et al. 2007) The idea of a thermosyphon dates back to 1836 with the invention of the Perkins Tube (Palm & Tengblad, 1996). The idea was thought up for transferring heat from a furnace to a baking oven with the hopes of producing more uniformly baked bread. More recently, thermosyphons have been put into practice in demanding thermal management application. However, the specific purpose of a thermosyphon is not necessarily for coo ling, but to transfer heat away from a heat source, where cooling can be achieved more conveniently by means of a secondary coolant. By implementing a multiphase evaporator plate to tr ansfer the heat away from the source, it becomes beneficial to incl ude an aircooled condenser with an extended surface and relatively low thermal resistance to c ondense the working fluid. Palm and Tengblad (1996) classify this type of thermosyphon as an advanced thermosyphon loop in a detailed review of thermosyphon loops. In this paper Palm and Tengblad (1996) clarify the reasons for thermosyphon system failure. They show that the heat transfer capacity of a simple 24 PAGE 27 therm osyphon loop with adequate filling is rest ricted only by the pres sure drop limit and the boiling limit. In the case of an advanced ther mosyphon loop, it is possible to enhance the boiling surface in order to prevent any system failure that results near the boiling limit. Therefore, the heat removal capabilities of a thermosyphon with an advanced configuration are only limited by the pressure drop limit, which entails flow oscill ations and eventually system failure due to dryout at relatively high heat flux. In order to reduce the thermal resistance w ithin the thermosyphon considered in this work, a microchannel plate will be used in the evaporator. The purpose of the microchannels is to enhance the thermal transport from the plat e to the working fluid, which will help the efficiency of the condenser by increasing the temperature difference between the vapor and ambient. Microchannels Within the past twenty years, there have b een a growing number of experimental studies on twophase flow and evaporation heat transfer in microchannels that suggest significant improvement in heat transfer ove r conventional channels. Pioneeri ng research in this field began about 2 decades ago with Tuckerman and Pease (1981). Applications such as microheat exchangers, microcooling assemblies and microthermalmechanical systems (MTMS) are rapidly advancing. While microchannel cooling ha s contributed largely to the success of these smaller scale systems, a lack of fundamental understanding still exists that limits design methods. Thome (2004) has provided a literature review on microchannel boi ling. He explains that there is currently no distin ct definition for the transition from conventional channels to microchannels. Mehendal et al. (2003) r ecommends a size classification as follows: microchannels (1100m ), mesochannels (100m to 1mm), macrochannels (16mm), and 25 PAGE 28 conventional channels ( 6mm) while Kandlikar (1997) recom mends the following classification and size range s: microchannels (50600hDm ), minichannels (600m to 3mm) and conventional Figure 13. Microchann el and conventional ch annel configurations. channels ( 3mm). Although there is no clear understanding of the transition, Thom e does discuss several explanations for th e error that arises when using conventional models to describe the heat transfer and pressure drop present in mi crochannels. In both twophase heat transfer and pressure drop correlations found in the literature for large scale configurations, there is a large reliance on experiment. The experime nts used to fit these correlat ions are more often than not found to be in a turbulent type flow regime. Do to the small length scales associated with microchannels; laminar flow is more frequently obse rved. Therefore, it is not fair to assume that the flow characteristics such as flow transition, vapor volume fraction and drift velocity, among others found in a laminar flow will accurately reflect those encountered in experiments with hD 26 PAGE 29 flows of a turbulen t nature. Othe r effects such as capillary (sur face tension) forces will also strongly effect flow transitions as well as flow structure. The flow structure observed in microchannels is described by Thome (2004), see Figure 13. He suggests that under high heat flux conditions, there exist three important regions, a liquid slug, elongated bubble, and dry zone region. This alternate flow configuration that is observed in microchannels results from the combination of vapor bubble departure diameter to channel size ratio bhdd and surface tension effects. While Mehendal et al. (2000) and Kandlikar (2001) have made preliminary estimates of the channel size in which these alternate flow struct ures are observed, it has been accepted that a proper threshold value at which a channel can be classified as micro will depend on the fluid properties, system pressure, surface characteristics, flow properties, etc. Kandlikar (2005) provides a detailed review on experimental work and observation of microchannel flow boiling. He concludes that nucleating bubbles are pr esent in flow boiling under high shear conditions during the bubbly, sl ug, and annular flow regimes. The bubbles depart into the flow as individua l bubbles unless their size is sm aller than the channel dimension normal to the nucleating surface. Small bubbles con tinue to grow until they eventually become confined by the channel walls unde r the confined flow pattern which is similar to a plug flow seen in conventional channels followed by an annular flow pattern at higher vapor qualities. Throughout each stage, Kandlikar argues that a combination of nucleate boiling and convective heat transfer are the dominant m echanisms responsible for thermal transport, as is the case in conventional channels. Kasza et al. (1997) use a high speed camera to observe boiling phenomenon at high bubble nucleation frequencies. Through flow visualization, he provides evidence that nucleate boiling within the thin liqui d film is present within the slug flow regime, 27 PAGE 30 supporting Kandlikars argum ent. Relying heavil y on experimental measurements, Kandlikar (2004) suggests a correlation for heat transfer within microchannels and compares it against a large database. For experiments with Reynolds number ranging from 400 to 3000, (laminar to turbulent transition type flow), the average deviation between the predicted heat transfer coefficient and the data is 21.9%. For low Reynolds number flows (100 Re410D ), Kandlikars correlation deviates 16.3% from th e data and for very low Reynolds number flows ( ), the average deviation from experim ental data is 17.3%. Re100 Thome et al. (2004) claim that transient ev aporation of the thin liquid film surrounding the elongated bubbles is the dominant heat tran sfer mechanism as opposed to nucleate boiling. Due to the strong surface tension forces present in microc hannels, he suggests there are nearly no stratification effects or orient ation effects, making the most important flow regimes bubbly, elongated bubble, annular, mist, and flows with partial dryout. For ev aporating flow, Thome reasons that the lifespan of the bubbly flow regime in microchannels is short lived and the dominating regimes are elongate d bubble followed by annular flow. Here Thome presents a threezone flow boiling heat transfer m odel to describe the evaporat ion of the thin liquid film that is present in the elongated bubble flow regime that also incl udes the effects of partial dryout, see Figure 13. He argues that even though the evaporation process is heat flux dependent, bubble nucleation is not the dominant form of ther mal transport. The eva poration of the thin liquid film is the dominant heat transfer mech anism within microchannels. This dynamic system is modeled based on the time averaged heat transfer that occurs from the cyclic pattern of each of the three zones. In conclusion, th is model predicts the transient variation in local heat transfer coefficient during the cyclic passage of (1) a liquid slug, (2) an eva porating elongated bubble, and (3) a vapor slug. There is strong dependency on the bubble fr equency and the liquid film 28 PAGE 31 geom etry at formation and dryout. He shows that the heat transfer in the elongated bubble zone, where a thin liquid film exists, promotes heat tran sfer that is an order of magnitude greater than that in the liquid slug region. In contrast the vapor slug, or dryout zone, provides negligible heat transfer. The relative length of each zone is crucial for proper modeling as they influence the time period for each zone to pass the point of observation for each cycle and thus the value of the local timeaveraged heat transfer coefficient. With experiment al data taken from 9 sources covering tube diameters ra nging from 0.51 to 3.1 mm, mass velocities of 50 to 564 2kgms, pressures ranging from 124 to 5776 kPa, heat fluxes from 5 to 178 2kWm, and vapor qualities in the range of 0.01 to 0.99, Th ome (2004) compares against his 3zone model and predicts 67% of the points within 30%. Excluding multichannel experiments, the model was able to predict 77% of the data to within 30%. An experimental investigati on on pressure drop within microchannels has been provided by Kandlikar (2005). He shows that in the case of multiple parallel channels, there are pressure fluctuations that are proportional to surface te mperature, and occasionally flow reversal is observed. Ribatski (2006) provide s a twophase pressure drop co rrelation in microchannel flow. In the laminar region he uses a homogeneous model while using MuellerSteinhagen (1986) friction correlation for turbulent flow. Ribatski (2006) comp ares this model against 2210 experimental twophase frictional pressure drop data points, and w ith exception to the transition region, provides reasonable agr eement. More than 85.7% of the data were correlated within 20% and over 96% within 30%. Lee and Lee (2001) provide a pressure drop correlation for twophase flow through sm all channels. Usi ng a modified LockhartMartinelli (1949) type correlation, Lee and Lee (2001) corre late all of the data within 10%. 29 PAGE 32 During flow boiling experiments near the lim iting case of a microchannel configuration, conventional two phase frictional pressure dr op and flow boiling models are capable of reasonably predicting the pressu re and heat transfer characteristics. Gungor and Winterton (1986)(1987) provide twophase heat transfer correlations that are applicable to subcooled flow. The first of the two is based on the convective heat transfer co efficient and the boiling number, while the latter uses the method of superpositio n similar to the Chen (1966) correlation, which couples the convective boiling heat transfer and the nucleate boiling heat transfer to obtain an overall 2phase heat transfer coefficient. Both agree reasonably with experiments for conventional channels. Kandlikar (1991) also presents a corre lation based on the boiling number that does a reasonably good job in predicting 2phase heat transfer coefficients in large scale channels. Several correlations for fric tional pressure drop through conv entional size channels have been presented in the literatu re. MuellerSteinhagen (1996) uses a method of superposition which couples the singlephase pressure drop of the vapor and liquid resp ectively. The average deviation between the predicted and the measured values are 41%. Other twophase frictional pressure drop models include Lockhart a nd Martinelli (1949) a nd Friedel (1979) which incorporate a 2phase frictional multiplier to ac count for the increased pressure drop within a twophase flow. Compact Condensers In order to design a high performance thermo syphon loop, the design and fabrication of a compact efficient condenser is necessary. The co ndensation process will transfer the latent heat of the vapor to the ambient air. There are a limited number of conf igurations available for the air cooled condensation process. On the airside it is most common to ha ve a platefin type arrangement. The addition of fins increases th e surface area and reduces the thermal resistance 30 PAGE 33 significan tly. More recently however, materials such as carbon and metal foams have been used in place of fins and have shown low airside ther mal resistance. With the goal of minimizing the air side pumping power and overall condenser volume, difficulty ar ises on the liquid side when choosing an arrangement that allo ws for the proper drainage. The most common configurations include the Serpentine/Paralle l Cross Flow Condenser, Downflow Condenser, and the Reflux Condenser; (Figure 14). Each configuration has its own benefits and drawbacks. Figure 14. Crossflow condenser. A) Downflow. B) Reflux. C) Serpentine. In the case of plain fins, as is shown in Figures 14a through 14c, there exist several copper tubes, either in serpentine or parallel, surrounded by an array of fins, usually alum inum or copper, which add surface area and pull heat away from the copper tubing. The mechanisms responsible for the thermal transpor t are illustrated in Figure 15. A B C 31 PAGE 34 Figure 15. Thermal transport mechanisms inside the condenser. Figure 16. Thermal resistance inside the condenser. As the heat transfers from the condensing fl uid to the ambient air, several forms of resistance to heat flow are encountered. The vapor enters the condenser and condenses along the tube wall, The latent heat is than conducted through the wall, to the fins where it passes through any added contact resistances that may re sult from the m anufacturing "condensationR"pipewallR 32 PAGE 35 process, The forced air passing over both the fins, "contactR" f ins R and tube banks, "base R allows the heat to dissipate to the ambien t. The resulting thermal resist ance is shown in Figure 16. Several innovative techniques have been im plemented over the years to significantly enhance heat transfer by reducing the thermal resistance. Starting w ith the manufacturing process, techniques such as electroplating and alum inum braising have been put into practice that minimize interfacial gaps that minimize the localized thermal resistance. The liquid side heat transfer enhancement has been investigated by Zuo et al. (2000) Applicable to both the downflow and reflux conf igurations, Zuo shows th at by adding drainage disks within the tube, the film thickness is re duced, allowing considerab le improvement in the heat transfer. However, flooding within the ch annels becomes increasingly more significant because of the increased liquid entrainment. While proper draina ge is highly beneficial to enhancing heat transfer, the adde d pressure drop that results clos e to the flooding limit within a reflux type condenser can make for a relatively inefficient process. The highest thermal resistance found within an aircooled co ndenser comes from the airside heat transfer coefficient. Even with the extended fin surface, the heat transfer is still constrained by airside limitations. By alteri ng the geometry and gene rating turbulence, the thermal boundary layer can be interrupted, allo wing adequate mixing between the plates and providing a significant increas e in heat removal. Fin geometries such as louver, offset strip and wavy fins are among the geometri es classified as interrupt ed surfaces that are commonly encountered in industrial applications due to th e enhanced heat transfer performance, ease of fabrication, and low cost. Fi gure 17 shows the louvered fin geometry and the corresponding boundary layer thickness along the plate. This noticeable decreas e in boundary layer thickness, keeps the ambient air well mixed, although the pre ssure drop is increased Proper modeling tools 33 PAGE 36 therefore require a firm understand ing of the heat transfer capabilities as well as the pressure drop characteristics for a given system. Several papers in the open litera ture compare the heat transfer performance of plain fins versus louver or offset strip fins claiming a vast improvement in efficiency for the interrupted surfaces. In cont rast, the performance of extended surfaces in the form of carbon or metal foams is not as frequently encountered in the lite rature. Therefore this study will focus on evaluating the heat transf er performance of carbon and metal foams for comparison with louvered fins. Figure 17. Flow field and boundary layer for louver type fins. Due to the vast amount of experiments conduc ted on louver type fins, there have been several empirical correlations pr esented in the literature that provide reasonable accuracy. Wang et al. (1999), Chang et al (1997) and Chang et al. (2006) have provided gene ralized heat transfer and friction correlations that are widely used and each has been tested against a large database consisting of 74, 91, and 91 data samples. The co rrelations were constructed from experiments in which the fin efficiencies are determined us ing the Schmidt (1949) approximation and the average local heat transfer coefficients are computed. The idea of using foam materials to enhance he at transfer is a relatively recent concept. Lu et al. (2006) and Zhao et al. (2006) conducted a thermal an alysis on a metal foam tubein34 PAGE 37 tube heat exchanger and com pared it to that of a conventional uninterr upted finned tube heat exchanger. They conclude that the use of metal foams can significantly enhance the heat transfer performance when compared to conventional plai n fin heat exchangers due to the increased surface area and mixing capabilitie s; however no information on pressure drop is provided. Boomsma et al (2003) conducted experiments using 50% waterethylene glycol solutions and claim that for a given pressure drop, metal fo ams are capable of removing 2 or 3 times higher heat loads. Implementing carbon foam into heat exchange r design has recently b een investigated by Klett (2002). He demonstrates that the lig ament conductivity of the carbon exceeds 1700 WmK, compared to aluminums 237 WmK and copper 401 WmK. High thermal conductivity foam drastically increases the fin efficiency and heat transfer rates. He shows that for a solid foam configuration, heat transfer ra tes compared to aluminum foams are enhanced by a factor of 10, while the pressure drop increases by a factor of 40. Other configurations such as through holes or fin shapes provide similar heat transfer results wh ile detailed pressure drop data are not provided. Instabilities and Subcooling Several reviews on twophase flow instabil ity phenomena exist within the literature including Boure et al. (1973), Ishi (1970), Bergles (1977), Yadi garoglu et al. (1981), Fukuda & Kobori (1979), and Kakac & Bon (2008). A genera l classification for inst abilities in forced convection twophase flow systems is de scribed in Table 12 (Boure et al. 1973). A more detailed description of the various instabilities ca n be found in Boure at al. (1973). Numerous studies have been conducted to observe the mechanisms responsible for the onset of instability under natural circulati on conditions (Kyung et al. 4996)(Aritome & Chang 35 PAGE 38 1993)(V an Bragt, 1998)(Yang et al. 2005)(Durga Prasad et al. 2007)(Jiand et al. 2000). For natural circulation systems, Yang et al. (2005), Durga Prasad et al. ( 2007), and Jiang et al. (2000), describe natural circulation instability unique to passively driven flow. In this case the instability is first triggered by the static instability (Ledinegg) before becoming dynamic. After the flow excursion, dynamics are introduced in th e form of both pressure drop (1995) and natural circulation instability. Pressure drop instab ilities result from the added dynamics of a compressible volume within the system. Once the Ledinegg excursion takes place, the vapor or gas within the system is compressed, and the build up of pressure will cause the flow rate to increase again. The result is an oscillatory flow rate, or a dynamic stat e. Natural circulation instability results from added dynamics that ar ise when hydrostatic head replaces the pumping mechanism to drive the flow. Once a Ledinegg excu rsion occurs, the flow rate is decreased, causing an increase in void fraction. This in tu rn results in an increase in hydrostatic head causing the flow rate to increase and the cycle is repeated causing oscillations in flow rate. Both pressure drop oscillations a nd natural circulation instability occur as secondary phenomenon triggered by the static instabi lity (Ledinegg). Tadrist (2006) provides a review on instability phenomena found in narrow spaces. He shows that instability can arise in small channels due to the onset of critical heat flux, CHF. The dr yout and rewetting cycle of the heated surface can result in vapor recoil within the channel. Depending on the compressibility of the inlet zone of the channel and the magnitude of density oscillat ions, quasiperiodical pre ssure fluctuations are observed. For ideal thermal management design in the case of a thermosyphon, one must pay close attention to system instability. By controlling the system pressure and bulk fluid temperature it is 36 PAGE 39 possible to alter the point of instability or im prove the heat transfer for flow rates encountered during stable operation. Table 12. Various types of twophase flow instability. Static Instabilities Dynamic Instabilities Ledinegg instability Density wave oscillations Boiling crisis Acoustic oscillations Bumping, chugging, or geysering Pressure drop oscillations Flow pattern transition Thermal oscillation Multichannel instabilities BWR instabilities As is well known for twophase flow, th e pressure drop can be double valued for a specified flow rate. This characteristic of twophase flow can give rise to the Ledinegg instability. In a given system, there exists a cr itical vapor quality in which the friction in the system is significant enough to cau se fluctuations in the flow rate ultimately resulting in dryout and system failure. Joshi et al (2005) provide a numerical study on twophase flow in a thermosyphon. They conclude that an increase in power or a decrease in subcooling results in a destabilizing effect on the natural circulation within the loop. He also concludes that the margin of instability increases for reduced riser height. The problems that arise from this unstable fl ow structure makes reducing pressure drop on the liquid side critical for operation at high heat flux. By lowering th e bulk fluid temperature entering the evaporator plate, it is possible to keep the vapor qu ality down and delay the onset of instability. On the other hand, decreasing the bul k fluid temperature requi res a larger condenser or a more powerful fan. Summary Due to the demand for high efficiency energy co nversion devices, there is current interest in improving the overall power and energy densitie s of PEM fuel cells. Major challenges arise 37 PAGE 40 38 with respect to thermal management, and new te chniques are needed to handle the high thermal loads introduced in high power density systems. The goal of th is research is to construct a flow boiling microchannel thermosyphon system that can be integrated into a PE M fuel cell stack to improve the overall power/energy density of the cell. By successfully measuring and modeling the thermalhydraulics of the system, this rese arch should allow for future optimization of a microchannel thermosyphon facility for fuel cell thermal management applications. PAGE 41 39 CHAPTER 3 FLOW BOILING MICROCHANNEL EV APORATOR AND THERMOSYPHON PERFORMANCE Background and Objectives Protonexchange membrane (PEM) fuel cells have recently been the focus of intense research for mobility applications. When considering fuel cells for space applications Burke (2003) suggests that fuel cell development should focus on both power density and energy density. The added value potentially provided by fuel cells to space science missions takes the form of mass savings, volume savings or more power or energy at the same mass and volume when compared with battery energy storage. Th e basic operating principle of a PEM fuel cell is that hydrogen fuel is oxidized at the anode, thus liberating electrons a nd producing protons. The free electrons flow to the cathode through an external circuit a nd combine with the protons and dissolved O2 to produce water and heat. The electric ci rcuit is completed as the protons flow from the anode to the cathode through a solid el ectrolyte membrane. As elucidated by Djilali and Lu (2006) and Berning and Djilali (2003) se veral important coupled mass, momentum, and heat transport problems occur simultaneously w ithin the fuel cell stack. Three particularly crucial issues that strongly influence PEM perf ormance are: 1) thermal management, 2) water management, and 3) mass transport limitations During current PEM ope ration, humidification of the reactant gas stream is us ed to control the membrane moisture content. At high current densities, product water is removed by air stream convection, and is controlled by adjusting pressure drop, temperature, and moisture content. Up to 50% of the heat generated during high current density PEM fuel cell operation must be removed to prevent dry out and excessive heating of the membrane, otherwise local degradation of the membrane will occur. The small temperature difference available between the fu el stack and the ambient environment renders PAGE 42 40 thermal management a very difficult issue. NASA Glenn has been developing unitized regenerative fuel cells for unmanned vehicle a pplications (Burke 2003). Currently the limiting factor to increasing the power densit y of their test fuel cell is an in ability to remove the heat load at higher current densit ies (Bents, 2004). Burke (2003) is currently working to develop a loop heat pipe evaporator, where low temperature heat sinks are ava ilable. The more conventional approach to thermal management is to utilize fo rced convection of liquid coolant through the fuel cell stack (Zhang et al. 2004). A major concern in fuel cell operation is that the required pumping energy for the coolant flow reduces the net power output of th e fuel cell (Wetton, 2004). The objective of this work is to develop a flow boiling microchannel cooling plate that may be easily inserted into the stack of a PEM fuel cell for thermal control. The current thermal control system operates in the natural circulati on mode where the flow is driven by twophase density gradients. You and Liu (2002) have developed a twophase flow transport model for the cathode side of PEM fuel cells. They have demonstrated that ther e exists an operating temperature for optimum performance. At low operating temperatur e the reduction reaction requires higher activation overpoten tial to produce the same current density. However, at higher temperatures the membrane resistance is higher. The optimum temperature tends to be in the range of 7080 C. Evaporative cooling inherently involves twophase flow and changing flow regimes. The ideal configurat ion would utilize a cooling medi um whose saturation temperature is approximately 10 to 20 degrees less than the optimum PEM stack operating temperature. Thus, once the PEM stack ramps up to its designe d operating temperature, sufficient superheat will be available for boiling, and a naturally driven flow will be established. At steadystate, the flow rate is established by a bala nce of the pressure head and flow resistance. At higher heat PAGE 43 41 loads more evaporation will occur and establish a la rger pressure head and flow resistance. The advantages of using phasechange for therma l management include: 1) the wall superheat associated with flow boiling does not change appr eciably over a large change in heat flux, and thus the membrane temperature may be kept near the optimum, 2) natura l circulation flow does not consume additional electric power to drive a pump, and 3) high heat fluxes may be achieved at low superheat, thus allowing for high power density operation. A poten tial disadvantage of operating in the natural circulati on mode is that the flow is su sceptible to a disrupting Ledinegg instability if the operatin g heat flux is too high. This work describes the thermal hydraulic performance of a flow boiling microchannel cooling plate operating within a closed loop twophase thermosyphon. The cooling medium used is HFE7100 which has a saturation temperature of 61 C at atmospheric pressure. Conventional flow boiling heat transfer correlations are found to be applicable to predict the wall superheat for an imposed heat flux in such microchannels. The mass flow rate through the thermosyphon can also be well predicted using st andard macroscale thermal hydraulic models. Although the Ledinegg instability could potentially limit satisfactor y operation of the system at very high heat flux, the present flow boiling micr ochannel cooling plate is found to be suitable for high power density PEM fuel cell applications since it can already handle a significantly larger operating heat flux than currently encount ered in industrial practi ce. This work should provide the tools necessary to accurately model the thermal hydraulics of the microchannel evaporator plate and give insight into the heat removal capabiliti es when integrated into a PEM fuel cell stack. PAGE 44 42Experimental Facility The experimental twophase thermosyphon is show n in Figure 21. The main test section consists of a 56 channel aluminum cooling plate. Each channel has a cro ss section of 1 x 1 mm and has a length of 115 mm, (the 1 mm. x 1 mm. cr oss section of the channels are classified by Kandlikar (2001) as minichannels while Mehendal (2000) classifies these at the threshold of mesochannels and conventional channels). Figure 21. Experimental twophase thermosyphon. The cooling plate is covered with Lexan to allow for flow visualization. An imposed heat flux is placed on the cooling plate using a Minco foil heater, 127 mm X 12 7 mm, that covers the back of the cooling plate, and the maximum heater flux is 70 kW/m2, based on the surface area of the heater. The heated surface area in contact with the liquid is 193.2 During operation, 2cm PAGE 45 43 the cooling plate is thoroughly in sulated with polyurethane foam insulation rated at a thermal conductivity of 0.04 WmK. The temperature difference between the insulation and ambient is used to evaluate the heat loss, which is based on calibration. By taking direct measurements of voltage and current, the total heat input is determ ined. The heat transferred to the working fluid is taken as the difference between the to tal heat input and the heat loss. Figure 22. Exploded view of the cooling plate assembly. An exploded view of the cooling plate assembly is shown in Figure 22. The facility is filled with HFE7100 which enters and discharges the cooling plate through an inlet and outlet manifold, respectively. The twophase mixture di scharges the cooling plate and rises to a water cooled condenser operating at 15 C with a water flow rate of up to 43 g/s. The vapor is condensed, and the pure liquid phas e flows out of the condenser a nd through a flow meter. The liquid discharge from the flow meter returns to the cooling plate. The inner diameter of the PAGE 46 44 tubing carrying fluid between the condenser and cooling plate is 9.5 mm. The vertical height between the condenser outlet and manifold inlet is 114 cm, while the vertical height between the outlet manifold and condenser in let is 102 cm. The test section and piping throughout the system are insulated with elastomeric foam. 6 cm. 9 cm. Pressure Taps Screen Screen Figure 23. Section vi ew of flow meter. Since the facility operates in the natural circ ulation mode, it is important to have a flow meter that does not add a significa nt pressure drop. Therefore, a flow meter was fabricated using two 29 gauge screen meshes as depicted in Figur e 23. The pressure drop across the screens is measured with a Validyne differential pressure transducer with a diaphragm capable of measuring a pressure drop up to 0.55 kPa. The maximum pressure drop recorded was 0.125 kPa which corresponds to a flow rate of 43 g/s. The relationship between flow rate and pressure drop was determined through calibration over a range of 50 g/s. The uncertainty of the flow rate measurement is explained in detail in AppendixB. PAGE 47 45 Figure 24a. Thermocouple lo cations on cooling plate. A second Validyne differential pressure transduc er is used to measure the pressure drop across the cooling channel and ha s a range of 03.5 kPa, respectively. Type E thermocouples have been inserted into the inlet and outlet mani folds to measure the fluid temperature. Also, 8 thermocouples have been embedded into the back of the cooling plate to measure wall temperature as depicted in Figures 24 a and b. The aluminum wall thickness between each thermocouple and the fluid is 1.1 mm. Th e thermocouples were bonded securely using Omegabond200 high thermal conductivity adhesive. Each thermocouple has been labeled 18 for reference. All calibrations and uncertainti es are explained in detail in Appendix B. PAGE 48 46 Figure 24b. Thermocouple locations on cooling plate. All signals were measured us ing digital data acquisition c onsisting of a Measurement and Computing CIOEXP32 multiplexer and PCIMDAS 1602/16 16 bit analog to digital converter. A custom algorithm was developed to convert the measured voltage s to pressure, temperature, and flow rate, see Appendix B for more details. Experimental Protocol HFE7100 is added to the facili ty through the top of the conde nser until the height of the fluid is just below the condens er coil. Liquid water at 15 C is then circulated through the condenser coil. The differential pressure transdu cers are then purged of ai r. The heat flux on the cooling plate is gradually raised so that vigorous boiling is ob served through the channels. Once vigorous boiling is established for a half an hour, the condenser is purged of noncondensable gas. Following purging, the heat flux is reduced to zero so that the fl ow through the facility PAGE 49 47 ceases. The heat flux is then raised in increments, and at each increm ent the flow rate and temperatures are allowed to reach a steady stat e. Measurements of pressure, flow rate, and temperature are made at each heat flux interval. The heat flux is raised until the flow becomes unstable, at which point large s cale flow oscillations and flow reversal are easily visible. Heat Flux (kW/m2) 05101520253035 T sub (C) 2 4 6 8 10 12 Figure 25. Degree of inlet sub c ooling variation with heat flux. The fluid leaving the condenser is subcooled to some degree. The degree of subcooling increases for increasing heat flux, as shown in Figure 25. The condenser coils are located above the liquid level within the conde nser. When the vapor condenses it is highly subcooled and mixes with the bulk liquid in the storage area. Therefore, as more vapor flows through the condenser, the bulk fluid b ecomes more subcooled. PAGE 50 48Results Figure 26 shows the measured flow ra te circulating through the twophase thermosyphon at various heat fluxes, (where the h eat flux is based on the heated area in contact with the liquid). Initially, there is a rapid rise in flow rate with increasing heat flux. This rapid rise in flow rate is due to an increase in vapor volume fraction in the riser, and thus the density difference between the down flow line and riser increases, which results in an increased flow rate. As the heat flux increases further, the incr ease in vapor volume fraction is marginal, but the increase in vapor quality is significant which re sults in increased friction. Thus, the mass flow rate peaks at 43 g/s, corresponding to a heat flux of approximately 24 kW/m2. As the heat flux is further increased, the mass flow rate decrease s. The maximum achievable heat flux during these sets of experiments is 32 kW/m2. At that point, the flow become s unstable due to the onset of a Ledinegg instability. 05101520253035 0 10 20 30 40 50 Figure 26. Variation of fl ow rate with heat flux. PAGE 51 49 Heat Flux (kW/m2) 05101520253035 Wall Temp (C) 64 66 68 70 72 74 76 78 80 82 84 86 Figure 27. Wall temperatur e at various locations. Figure 27 shows the measured wall temper atures for thermocouple positions 18 at various heat fluxes. The wall temperatures n ear the entrance to the cooling plate, which correspond to locations 4 & 8 from Figure 4b, tend to be the highest. This is because the fluid enters the cooling plate as a pure liquid, and the h eat transfer coefficient is lower than that near the exit of the cooling plate, where the vapor quality is relatively high. Also, the wall temperatures are slightly less near the side of the plate than at the central region. This is due to the fact that there is finite h eat leakage at the edges of the plate. In general, the thermal performance of the twophase thermal management system is quite good. The wall temperatures of the cooling plate vary from approximately 66 to 82 C depending on the heat flux and location. This is very close to the optimum temperature range for PEM fuel cell membranes. PAGE 52 50 Furthermore, PEM fuel cells currentl y only require a heat flux of 5 kW/m2 (Lasbet, 2006) heat removal capacity. Therefore, the present system far exceeds the current need and has sufficient heat removal capacity for larger power density systems prior to th e onset of Ledinegg instability. It should be noted, that the thermal performance of the system can be tuned by raising or lowering the height between the microchannel co oling plate and the condenser. Increasing the height will tend to increase the flow rate, and the vapor quality and heat transfer will adjust accordingly. In order to elucidate how the variation of va por quality influences the thermal field along the cooling plate, the axial vari ation of thermodynamic vapor qua lity is computed. Since the flow is subcooled entering the cooling plate, the change in sensible heat must be considered when evaluating the thermodynamic equilibrium vapor quality, x z, at any zlocation along the cooling plate; therefore from an energy ba lance the axial vapor quality is computed, (()) () z pi fgqmCTzT xz mh (21) where z q is the rate of heat flow into the cooling plat e from the entrance to the zlocation, T(z) is the bulk fluid temperature at the zlocation, Tin is the inlet fluid temperature, is the mass flow rate, and m f gh is the latent heat of vaporization for the fluid. In order to evaluate T(z), the inlet and exit fluid temperatures are measured, and it is assumed that there is a linear change in the bulk fluid temperature along the zdirection. This assumption ma y not be exact, but since the maximum temperature change betwee n the inlet and outlet is only 10 C, it is acceptable for the estimation of the vapor quality x(z). The nonuniform heating is apparent from the measured wall temperatures shown in Figure 27. The late ral heat flux toward th e edges of the cooling plate is approximately 27% of the heat flux into the fluid. Figure 28 shows the variation of PAGE 53 51 vapor quality along the zdirecti on for various heat fluxes at th e central region of the cooling plate. While the quality appears to be small fo r low heat flux, this approximation is based on thermodynamic quality and does not account for th e localized subcooled boiling which produces net vapor generation at negative thermodynamic e quilibrium quality. Also, the quality increases rather slowly for low heat flux. At higher heat flux, the flow begins to decr ease in flow rate and a rapid increase in vapor quality is observed. z (mm) 020406080100120140 Vapor Quality 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 q"=32 kW/m2 q"=28 kW/m2 q"=24 kW/m2 Figure 28. Variation in vapor quality with channel height. Flow Modeling In a closed loop natural circul ation system, the pressure at the discharge of the condenser is the same as that at the return of the condenser Therefore, in order to predict the system flow rate it is necessary to comput e the pressure change around the flow loop and converge on a flow PAGE 54 52 rate that insures the inlet and di scharge pressure at the condenser match. The two phase pressure gradient may be expressed as, 2 g fadPdPdPdP dzdzdzdz (22) It consists of three components, gravitationa l, frictional, and accel erational, denoted by subscripts g, f, and a respectively. Each must be modeled accordingly to ensure proper representation of the flow behavi or. The gravitational pressure gr adient results from a body force acting on the fluid and is given by, sinm gdP g dz (23) where m is the mixture density, g is th e gravitational acceleration, and is an inclination angle. For upflow = /2 and =/2 for downflow. The mixture de nsity depends on the vapor volume fraction (1)mv l (24) where is the density and the subscripts l and refer to the liquid a nd vapor, respectively. In order to estim ate the vapor volume fraction, the standard ZuberFindlay (1965) drift flux model for slug flow is used, v 1 01 1vv vj lx CV xGx (25) 1/4 21.53lv vj lVg (26) where G is the mass flux and C0=1.2. The pressure gradient asso ciated with the acceleration of the fluid is derived from a momentum balance on a differential section of the tube in a constant PAGE 55 53 cross section. The resulting expression for the cha nge in velocity due to the change in vapor quality and vapor volume fraction is, 2 2 21 1lv ax dPdx G dzdz (27) The frictional pressure gradient is modeled us ing the correlation from MuellerSteinhagen and Heck (1986). This correlation co nsiders the singlephase fricti onal pressure gradient for both liquid and vapor and uses supe rposition along with an empirical fit to model the twophase frictional pressure gradient, 1/3 31fdP F xB dz x (28) 2 F ABAx (29) where the superficial velocity of the liquid a nd vapor phases is used to evaluate A & B as, 22l lG Af d (210) 22v vG Bf d (211) Here, d is the hydraulic diameter, and the Darcy friction factor is computed as, 1/464 ,Re1187 Re 0.3164 ,Re1187 Rel l l l lf (212) 1/464 ,Re1187 Re 0.3164 ,Re1187 Rev v v v vf (213) Rel lGd (214) PAGE 56 54 Rev vGd (215) Flow Rate (g/s) 01020304050 Pressure Drop (Pa) 200 0 200 400 600 800 1000 1200 1400 1600 1800 Measured Figure 29. Measured and predic ted microchannel pressure dr op variation with flow rate. At high heat flux the frictiona l pressure gradient within the microchannel plate is approximately 50 times the pressure gradient in the liquid section of the thermosyphon loop and 2.5 that in the twophase upper se ction. To test the correlation of MuellerSteinhagen and Heck (1986), the total measured pressure drop across the microchannel plate, which is dominated by friction, is compared against that predicted. As observed in Figure 29, the predicted pressure drop is in reasonable agreement with that measured. The model of MuellerSteinhagen (1987), commonly used to predict 2phase frictional pressu re drop within conventional sized channels, is used here and found to be satisfactory. The 1 mm. x 1 mm. cross sectio n of the channels under investigation are classified by Kandlikar ( 2001) as minichannels while Mehendal (2000) PAGE 57 55 classifies these at the threshol d of mesochannels and conventi onal channels. Therefore it is understandable that the MuellerSteinhagen mo del works well in pr edicting the frictional pressure drop in this configuration. Heat Flux (kW/m2) 0510152025303540 Flow Rate (g/s) 0 10 20 30 40 50 Figure 210. Comparison between meas ured and predicted flow rate. As is well known for twophase flow, the pressure drop can be double valued for a specified flow rate. This characteristic of twophase flow can give rise to the Ledinegg instability. In regionI, where th e slope of the pressure drop as a function of flow rate curve is positive, a small perturbation decreasing the flow rate will result in increasing vapor quality and larger gravitational pumping poten tial. Simultaneously, the frictional pressure drop decreases due to decreased flow rate, and thus the flow rate will adjust back to its original value. PAGE 58 56 Likewise, in regionI, the flow rate will return to its original value following a perturbation increasing the flow rate. The flow resistance, defined as P R m in this region is positive, Heat Flux (kW/m2) 05101520253035 Wall Temperature (C) 50 55 60 65 70 75 80 85 GungorWinterton (1986) MuellerSteinhagen (1995) Measured Figure 211a. Comparison between measured a nd predicted wall temper ature at thermocouple location 5 based on measured flow rate and vapor quality. thus the flow is stable in regionI. In contrast in region II the slope of the pressure drop as a function of flow rate curve, R, is negative and a small perturbation increasing the flow rate will result in decreasing vapor quality and gravit ational pumping potential. Simultaneously, the frictional pressure drop is reduced. Thus further increases in flow rate and decrease in vapor quality can result, until the flow swings back to region I. However, the flow state does not remain in region I for long due to the high heat flux and rapid increase in vapor quality. This behavior continues back and forth in an oscillatory mode until eventual dryout. In regionII, PAGE 59 57 where the resistance is negative, power is bei ng produced instead of dissipated; therefore the possibility exists for unstable flow. In the current experiment, this os cillatory behavior was observed at a heat flux of 32 kW/m2. Further explanation of the ons et of Ledinegg instability for microchannel systems will be discussed in Chapter 4. Heat Flux (kW/m2) 05101520253035 Wall Temp (C) 50 55 60 65 70 75 80 85 GungorWinterton (1986) MuellerSteinhagen (1995) Measured Figure 211b. Comparison between measured a nd predicted wall temper ature at thermocouple location 6 based on measured flow rate and vapor quality. In order to compute the mass flow rate for a specified heat flux, an algorithm has been developed that searches for the flow rate that yields an equal pre ssure at the inlet and exit of the condenser. Figure 210 shows a comparison between the measured and predicted mass flow rate at various heat fluxes assuming a linear bulk fl uid temperature profile along the microchannel plate and a nonlinear temp erature profile using the Saha & Z uber (1974) correlation. Very little difference is observed between the two predictions. In general, th e flow rate is slightly under PAGE 60 58 predicted, but considering the complexity of the sy stem, the predicted flow rate is satisfactory for design purposes. Figures 211 a and b compare the measured wall temperatures at thermocouple locations 5 and 6, respectively, with those predicte d using the GungorWinterton (1987) and MuellerSteinhagen (1996) subcooled flow boiling correlations. The computations for these comparisons use the measured flow rate and vapor quality. The GungorWinterton (19 87) correlation gives better agreement at lower heat flux and vapor qual ity, and in general th e data lie between both predictions. The wall temperatures were also computed using the predicted mass flow rates, and the results are virtually identical to those shown in Figures 211 a and b. This indicates that the predicted wall temperatures are not sensitive to small variations in vapor quality and mass flow rate. It is also noteworthy that flow boiling correlations devel oped for large diameter tubes are also satisfactory for the current microchannel plate design. The thermalhydraulic models presented in this work allow the thermal field in the microchannel plate to be predicted from fundamental models without the need for experime ntal inputs. While thes e correlations provide satisfactory results for the current configuration, it is unlikely smaller sized channels will be correlated accurately using the models discussed here. PAGE 61 59 CHAPTER 3 MICROCHANNEL PRESSURE DROP CORRELATION Background and Objectives A lot of attention has been paid to the twopha se flow pressure drop in small channels as the ability to fabricate flow channels at incr easingly smaller sizes continues to improve. New technologies allowing fabrication of microchannels on the order of a few micrometers (Hwang) have found to be helpful in a variety of advanced applications including microelectromechanical systems (MEMS), fuel cells, electronics cooling, and others. In such applications, it is beneficial from a thermal management perspective to fabri cate small channels due to an ability to yield high heat transfer coefficients and allow for a compact design. Since the reduction in channel size often leads to an undesira ble increase in pressure drop, it is important to have an understanding of the transport phenomena governing pressure drop within th ese small channels. A limited number of work exists in the literature describing pressure drop in small channels (Lee & Lee, 2001) (Lowry, 1988) (Id e, 1990) (Wambsganss, 1991) (Wambsganss, 1992) (Mishima, 1993) (Fujita, 1995) (Mishima, 1996) (Kandlikar 2005) (Ribatski, 2006), while others exist for conventional channels (MullerSte inhagen, 1986) (Friedal, 1979) (Lockhart & Martinelli, 1949). The goal of this work is to experimentally measure both the single phase and multiphase pressure drop across a microchannel plate over a wide parameter space and compare against various pressure drop correlations av ailable in the literature. This work is expected to contribute to the overall understanding of transp ort phenomenon within microchannels. Experimental Facility The experimental facility is s hown in Fig. 31. The main test section consists of a 56 channel aluminum microchannel plate. Each cha nnel has a cross section of 1 x 1 mm. and has a length of 115 mm, (the 1 mm. x 1 mm. cross section of the channels are classified by Kandlikar PAGE 62 60 (2001) as minichannels while Mehe ndal (2000) classifies these at the threshold of mesochannels and conventional channels). The microchannel pl ate is covered with Lexan to allow for flow visualization. A Minco foil heater, 127 mm. x 12 7 mm, is in direct contact with the bottom of the aluminum plate that provi des a maximum heat flux of 70 2kWm. A micropump is integrated into the loop to allow for flow rate control. During ope ration HFE7100 circulates between the microchannel plate and the condens er. A more detailed description of the experimental facility is given in Chapter 2. Figure 31. Thermosyphon experimental facility. Experimental Protocol During operation HFE7100 is first added to the system until the liquid height is just below the condenser coils. All pressure lines are than purged of air to allow fo r accurate measurements. PAGE 63 61 The pump is than increased in increments to allow for flow through the system. At each increment of pumping power, the heat flux is also raised in increments to allow for varying vapor quality. By controlling both the pumping power and the heat flux, a large parameter space is obtained with respect to flow rate and vapor quality. At each incr ement, the pressure drop, flow rate, heat input, inlet temperature, and outlet temp erature are recorded after the system is allowed to reach a quasisteady state condition. Results The important parameters governing pressure drop aside from the geometric dimensions are the flow rate and vapor quality. During operation, the heat input, flow rate, m, inlet temperature, and outlet temperature, are directly measured. Using these direct m easurements the vapor quality can be extracted at each axial lo cation. First the temperature field is computed using the SahaZuber (1974) model for subcooled boiling. Once this is obtained, the variation in vapor quality in the axial direction is obtained from Eq. 21. qinToutTThe next ste p is to extract the frictional component of pressure drop from the total pressure drop measurements that were measured dire ctly. This is done by subtracting both the accelerational and gravitational pressure drops from the measured pressure drops as shown in Eq. 22. The details on computing both the accelera tional and gravitational pr essure drop are given in Eqs. 24 through 27. Table 31 summarizes th e experimentally measured frictional pressure drop data for different heat inputs and mass flow rates. PAGE 64 62 Table 31. Frictional pressure drop data. q m f P inT outT q m f P inT outT W kgs Pa C CW kgs Pa CC 2.6 1.0 111 22.0 22.2 890.3 12.3 1069 47.8 60.8 2.6 1.0 116 21.9 22.1 0.5 0.4 34 21.3 21.6 57.3 4.7 123 47.4 56.8 0.5 0.5 35 21.4 21.6 57.3 4.9 36 47.4 56.8 11.3 1.7 19 30.8 36.7 147.4 21.7 370 56.7 61.7 11.2 1.7 19 30.9 36.7 208.4 25.8 627 56.2 61.5 28.6 2.5 46 40.3 49.6 208.4 25.7 630 56.3 61.5 54.9 6.8 20 52.3 58.9 276.5 27.1 801 55.1 61.0 54.9 6.7 20 52.3 58.8 276.4 27.1 800 55.2 61.1 91.2 15.3 285 56.7 61.0 350.8 27.2 922 54.1 61.0 91.2 15.3 286 56.6 60.9 350.8 27.3 918 54.2 61.0 148.9 20.8 627 56.9 60.9 441.0 26.5 1034 53.0 61.0 148.9 20.8 634 57.0 60.9 441.0 26.6 1035 53.0 61.1 0.3 0.6 26 21.7 21.9 533.9 24.7 1150 51.5 60.9 0.3 0.6 24 21.7 21.9 533.9 24.7 1152 51.5 60.9 208.6 21.9 756 56.0 60.6 645.5 21.6 1308 50.1 61.2 208.5 21.8 752 55.9 60.6 645.6 21.5 1309 50.1 61.2 278.0 22.0 856 54.9 60.6 738.8 18.2 1354 49.3 61.2 277.9 22.0 858 55.0 60.6 738.7 18.3 1354 49.3 61.2 345.0 21.5 930 54.0 60.7 861.2 14.7 1268 48.8 61.2 345.0 21.5 931 54.0 60.7 861.2 14.7 1266 48.9 61.2 435.8 20.5 1019 52.1 60.5 998.5 12.3 985 49.2 61.3 435.7 20.5 1018 52.1 60.5 998.5 12.2 982 49.1 61.3 523.9 18.5 1096 50.7 60.5 0.6 0.2 46 23.3 23.4 523.9 18.4 1100 50.7 60.5 0.6 0.2 45 23.3 23.4 626.6 15.8 1130 49.2 60.8 12.2 2.1 21 30.5 35.6 626.6 15.8 1132 49.2 60.8 12.2 2.1 21 30.5 35.6 749.4 12.5 1044 48.4 60.8 29.5 3.0 17 39.6 47.8 749.4 12.5 1043 48.5 60.8 29.5 3.0 17 39.5 47.7 773.8 12.0 1008 48.6 61.0 12.3 1.9 20 30.5 35.8 773.8 12.0 1009 48.5 61.0 12.3 1.9 20 30.5 35.8 12.1 1.6 18 30.8 37.2 26.9 2.6 14 39.4 47.9 12.0 1.6 18 30.8 37.3 26.9 2.7 14 39.4 47.8 23.9 2.9 14 40.9 50.8 57.5 6.6 24 50.7 58.1 23.9 2.1 9 40.9 50.5 57.5 6.4 19 50.7 58.1 48.4 6.8 21 52.1 58.5 95.9 15.6 178 55.8 60.5 48.4 6.8 22 52.2 58.5 95.9 15.6 158 55.8 60.6 4.1 0.6 197 22.6 22.4 147.8 22.0 501 56.4 60.8 4.2 0.5 197 22.6 22.4 PAGE 65 63 Table 31 Continued. q m f P inT outT q m f P inT outT W kgs Pa C CW kgs Pa CC 203.2 24.2 692 55.8 60.7 94.0 16.5 429 56.2 60.0 203.2 24.2 691 55.8 60.7 139.2 18.1 636 55.7 59.8 271.6 24.8 836 54.8 60.6 139.2 18.0 642 55.9 59.9 271.5 24.8 836 54.9 60.6 193.0 18.9 743 55.2 59.8 352.4 24.8 944 53.8 60.8 193.0 18.9 741 55.2 59.7 352.3 24.8 937 53.9 60.8 3.9 0.4 192 22.7 22.6 432.8 23.9 1042 52.4 60.5 3.9 0.4 191 22.7 22.6 432.8 24.0 1044 52.4 60.5 337.6 18.1 874 52.3 59.5 540.6 21.7 1171 50.9 60.6 337.5 18.0 875 52.5 59.7 540.6 21.7 1167 50.8 60.6 431.4 16.6 940 50.6 59.8 656.2 18.4 1264 49.4 60.8 535.9 14.2 972 48.2 59.5 656.2 18.4 1262 49.4 60.8 535.9 14.1 973 48.2 59.5 780.1 14.7 1225 48.1 60.7 627.1 11.8 929 47.3 60.0 780.0 14.6 1225 48.1 60.7 627.0 11.8 931 47.3 60.0 890.4 12.3 1065 47.9 60.8 PAGE 66 64Pressure Drop Correlations Several correlations exist that describe twophase pressure gradients for flow through channels. The conventional twophase frictional pressu re drop correlations considered here are those by Freidal (1979), MullerS teinhagen (1986), and the ho mogenous approximation, while that of Lee (2001) is specific to small channe ls. The MullerSteinhagen correlation is described in Eqs. 2.82.15, while those for Friedal, and Lee are shown below in Eqns. 31 through 323. Friedal 2 f r fldPdP dzdz (31) 2 22 1l l hlf dP mx dzd (32) 0.2516 Re2000 Re 0.079 Re2000 Red l df (33) Redmd (34) 2 0.0450.0353.24fr HlFH E FrWe (35) 2 21lv vl f Exx f (36) 0.224 0.781 Fxx (37) 0.910.19 0.71lvv vllH (38) l H md We (39) 11H vl x x (310) For 0 PAGE 67 65 2 12 2 255 1900 2400 5001900 4.8500 Gkgm G s B Gkgm G Gkgm s s (311) For 9.5 PAGE 68 66 22l l lf dP G dzd (322) 22v v vf dP G dzd (323) After computing the frictional pressure drops us ing each correlation based on the experimentally measured flow rates and vapor qualities, predic tions using each model are plotted against the experimentally measured values in Figs. 32 through 35. The m ean relative error, defined as, 2Pmp m estPP P N (324) along with the maximum relative error are reported for each correlation in Table 32. As shown the Muller Steinhagen (1986) and homogeneous m odels have the lowest mean relative error while the homogeneous model pr ovides the lowest maximum rela tive error. The microchannel frictional pressure drop correlation of Lee (2001) significantly over predic ts pressure drop at high vapor quality while those for conventional channels provide better overall agreement. This observation however may not be applicable to smaller channe ls, as the 1 mm x 1 mm cross section of these channels is large enough where it may be mo deled as a conventional channel under certain flow conditions, (the 1 mm. x 1 mm. cross section of the channels under investigation are classified by Kandlikar ( 2001) as minichannels while Mehendal (2000) classifies these at the th reshold of mesochannels a nd conventional channels). PAGE 69 67 0 500 1000 1500 0 500 1000 1500 2000 2500 3000 Measured Pressure Drop (Pa)Predicted Pressure Drop (Pa) 40% 40% Figure 32. Experim ental versus predicted frictional pressure drop using Friedal correlation. 0 500 1000 1500 0 500 1000 1500 2000 2500 Measured Pressure Drop (Pa)Predicted Pressure Drop (Pa) 40% 40% Figure 33. Experim ental versus predicted frictional pressure drop using MuellerSteinhagen correlation. PAGE 70 68 0 500 1000 1500 0 500 1000 1500 2000 2500 Measured Pressure Drop (Pa)Predicted Pressure Drop (Pa) Figure 34. Experimental versus predicted frictional pressure drop using homogeneous model. 0 500 1000 1500 0 500 1000 1500 2000 2500 3000 3500 4000 Measured Pressure Drop (Pa)Predicted Pressure Drop (Pa) 40% 40% 40% 40% Figure 35. Experim ental versus predicted frictional pressure drop using Lee correlation. PAGE 71 69 Table 32. Error between measur ed and predicted frictional pre ssure drop for various models. Model Max Relative Error Mean Relative Error Homogeneous 53% 21% Freidal 318% 64% MullerSteinhagen 153% 27% Lee 867% 82% Conclusion The frictional pressure drop has been expe rimentally measured for twophase flow through a microchannel plate and co mpared against various correlati ons within the literature. The largest error occurs in the correlations obtained from experiment on pressure drop through small channels while the best correlations are those used for conventiona l size channels. Both the homogenous model as well as the frictional pr essure gradient model of MuellerSteinhagen provide the best overall agreemen t. Due to the chaotic nature of twophase flow, along with the uncertainties that arise within small channels, th e relative error is still quite large and further work should be done to better co rrelate twophase frictional pressure gradient within small channels. PAGE 72 CHAP TER 4 INSTABILITY PHENOMENA IN A TWOPHASE THERMOSYPHON Background and Objectives Twophase flow system instabilities and ther mohydraulic oscillations have long been of interest in the energy field since many steam genera tors operate in the natural circulation mode. The twophase instability phenomenon was first reported by Ledinegg (1938) when investigating twophase flow in steam generators. The onset of instabilities can be f ound in other applications involving passively driven twophase flow systems, such as thermosyphons for electronic cooling and the present application, fuel cell th ermal management. A number of reviews have considered twophase flow instability phenomena, including those of Boure et al. (1973), Ishi (1970), Bergles (1977), Yadigaro glu et al. (1981), F ukuda & Kobori (1979), and Kakac & Bon (2008). Instabilities in forced c onvection twophase flow systems can be categorized as static or dynamic. Static instabilities include Ledinegg, bo iling crisis, flow pattern transitions, bumping, chugging, or geysering. Dynamic instabilities include density wa ve, acoustic, pressure drop, and thermal oscillations, in addition to multichannel and BWR. Boure et al. (1973) provides a more detail ed description on the various instabilities reported. Numerous studies have been conducted to observe the mechanisms responsible for the onset of instability under natural circulation conditions (Kyung et al. 1996) (Aritome & Chang 1993) (Van Bragt, 1998) (Yang et al. 2005) (Durga Prasad et al. 2007) (Jiang et al. 2000). For natural circulation systems, Yang et al. (2005), Durga Prasad et al. ( 2007), and Jiang et al. (2000), describe instabili ties unique to passively driven flow which are initially triggered by the static Ledinegg instability prior to transitioning to a pressure dr op (Lui et al. 1995) and natural circulation instability. Inst ability phenomena associated w ith twophase flow through narrow spaces is described by Tadrist (2006). It is shown that instability can arise in small channels due 70 PAGE 73 to the onset of local dryout. A more detailed explanation of th e instability identification for the current study can be found in AppendixA. An experimental investigati on that focuses on understandin g and predicting the onset of instability within a twophase microchannel thermosyphon facility is presented. Since it is a passive flow system, the emphasis is on the static Ledinegg instability. The experimental facility incorporates a fixed geometry with variable condenser height. The applied heat flux is the control variable, while the system flow rate is a dependent variable Under these operating conditions, small scale flow perturbations that occur under quasisteady state conditions are measured and used to predict the onset of inst ability. To validate these predictions, both the pressure drop as well as the limiting heat flux are compared against those determined experimentally. The limiting heat flux is satisfactorily predicted to within 10% uncertainty. The uncertainty in the model prediction arises from uncertainties in the frictional pressure gradient and vapor vol ume fraction models. Experimental Considerations Facility The experimental twophase thermosyphon is shown in Figure 21. The cooling plate assembly includes 56 microchannels, 115 mm in length, and having a 1 mm x 1 mm cross section. A Minco foil heater is positioned on the backside of the microchannel plate and provides a maximum uniform heat flux of 70 2kWm. A lexan cover plate is stacked on top of the microchannels and allows for flow visual ization. During operation, HFE7100 circulates between the cooling plate and conde nser. The condenser may be ra ised or lowered in order to change the gravitational head driving the flow through the facility. A mo re detailed description of the experimental facility is reported in Chapter 2. 71 PAGE 74 Experimental Protocol The condenser is positioned above the evaporator plate at one of four different heights considered for this investigation, 1.33, 1.15, 0.97, & 0.79 0.02 m. HFE7100 is first added to the system until the liquid height is just below the condenser coils All pressure lines are then carefully inspected for air bubbles followed by the purging of all pressure transducers. A heat load is applied to the evaporator plate to generate vi gorous boiling. The pressure release valve, located at the top of the condenser, is opened to purge all nonc ondensable gas from the system and improve condenser performance. The heat load is then reduced to ze ro so that the flow through the facility ceases. The h eat flux is raised in small increments. At each increment the flow rate and temperatures are allowed to reach a quasisteady state. Measurements of pressure, flow rate, and temperature are made at each heat fl ux interval. The heat flux is raised until large scale flow oscillations and flow re versal are easily visible; it is th is flow condition that is deemed as unstable. The condenser he ight is then repositioned and the experiment is repeated. System Subcooling and Manifold Pressure Drop The temperature of the working fluid leaving th e condenser is subcooled to some degree. At relatively low heat flux where no boiling take s place, the fluid is significantly below the saturation temperature, Tsub=1040 With a heat load of approximately 5 C 2kWm, vapor bubbles are observed in the riser. The vapor enters the condenser a nd back into the bulk liquid. Fig. 2 shows the degree of subcooling measur ed for each heat flux and condenser height considered. At higher heat flux more vapor comes in contact with the coils, and as the heat flux is increased, the degree of subc ooling becomes larger. Since th e degree of subcooling affects vapor generation, the onset of instability depends on the inlet temperature to the evaporator. 72 PAGE 75 Measurem ents of subcooling have been made at each heat flux interval, and are shown in Figure 41. 0 100 200 300 400 500 600 700 800 900 1000 0 5 10 15 20 25 30 35 40 Heat Flux (kW/m2)Tsub (C) H=1.33 m H=1.15 m H=0.97 m H=0.79 m Figure 41. Variation of inle t subcooling with heat flux. The pressure drops across the evaporator plate manifolds contribute significantly to the total pressure drop in the riser portion of the thermosyphon. Since the manifold design is unique to this experiment, direct measurements of pressure drop have been made at each heat load interval and condenser height and are shown in Figure 42. The pressure drop through the inlet manifold correlates with liquid velocity and th at for the exit manifold correlates with vapor superficial velocity and are shown in Figure 42. The vapor supe rficial velocity is, v vGxz u, (41) where G is the mass flux, x z is the vapor quality at an axia l location along the evaporator plate and v is the vapor density. The vapor quality is dete rmined from Eq. (21). The vapor quality is determined at the exit of the evaporator plate using the measured exit temperature of the bulk fluid, in order to compute the exit manifold superficial velocity. As shown the inlet manifold pressur e drop is larger than the outlet mani fold where twophase flow is present. For eT 73 PAGE 76 sim plification, the pressure taps for the inlet ma nifold were located so that the pressure drop across the flow meter could be included in the m easurement. It should also be noted that the pressure loss across the inlet manifold shown in Fig. 42 includes a small amount of tubing and a 90 degree bend. Inlet Manifold Liqui d Velocity (m/s) 0.000.050.100.150.200.250.300.350.40Pressure Drop (Pa) 0 1000 2000 3000 4000 5000 6000 Outlet Manifold Superficial Vapor Velocity (m/s) 024681012141 6 Outlet Manifold Inlet Manifold Measured Fit Figure 42. Pressure drop acro ss inlet and outlet manifolds. Flow Modeling To predict the flow rate circulating th rough the thermosyphon loop, the momentum equation is used to compute the pressure gradient The twophase pressure gradient is comprised of three components: that due to gravity, friction, a nd acceleration denoted by the respective subscripts g, f, and a in Eq. 2.2. In the case of a closed circulation loop, th e difference in pressure between the discharge and return of the condenser is considered negligibly small. Therefore computing the pressure field throughout the thermosyphon loop allows determination of the mass flow rate, m. The 74 PAGE 77 flow rate is com puted as that which gives th e same pressure at the inlet and outlet of the condenser. The gravitational pressure gradie nt component depends on the liquidvapor mixture density and the gravity vector is defined in Eq. 2.3 where m is the mixture density, g is gravity, and is the inclination angle with respect to the hori zontal. The mixture density is given in Eq. 24. The vapor volume fraction, is estimated using the ZuberFindlay (1965) drift flux model for slug flow is shown in Eqs. 25 through 26. Here =1.2 is the distribution parameter. A mom entum balance on a differential element is used to determine the accelerational pressure gradient component as described in Eq. 2.7. As suggested by Ga rrity et al. (2007), the frictional pressure gradient is computed using the MuellerSteinhagen a nd Heck (1986) correlation found in Eqs. 2.8.11. Here d is the hydraulic diameter and the Darcy friction factors are computed as shown in Eqs. 2.12.15. 0CT o compute the local vapor quality, x z, from Eq. 2.1, the temperature field, Tz is computed using the SahaZuber (1974) model for subcooled boiling. Using Eqs. 2.12.15 and the empirically determined pressure drop across the manifolds, the change in pressure from the condenser discharge to the inlet is computed. The mass flow rate is taken as that which yields no net change in total system pressure, s P This is expressed as 2 sdP Pd dz z=0. (42) 75 PAGE 78 Experimental Results Heat Flux (kW/m2) 01020304050607 0 Flow Rate (kg/s) 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 1.33 1.15 0.97 0.79 Measured Predicted H(m) Figure 43. Mass flow rate varia tion with increasing heat flux. The variation of the steady st ate mass flow rate of HFE7100 circulating through the twophase thermosyphon system with various applied h eat fluxes is shown in Figure 43 for different condenser elevations. At low heat loads the flow rate increases rapidly w ith increasing heat flux due to the increasing difference in pressure he ad between the riser and downcomer. This increase in pressure head results from the increasing vapor volume fraction. At slightly higher heat flux, the vapor quality becomes more signifi cant and an increase in frictional pressure drop is observed. As the frictional pressure drop beco mes more substantial, the flow rate begins to decrease. Although the mixture de nsity in the riser is low, th e high velocity and frictional dissipation results in a decreasing flow rate. The largest flow ra te occurs when the condenser elevation above the microchannel evaporator plate is highest. This is because a larger difference in pressure head between the ri ser and downcomer is achieved wh en the condenser elevation is largest. The pressure drop variation across the microchannel plate with resp ect to mass flow rate 76 PAGE 79 is shown in Figure 44. For H = 0.97 m and H = 0.79 m, the predicted pressure drop decreases with increasing flow rate under relatively low heat fl ux. This is a result of operating at low vapor quality, where the frictional pressu re drop approaches that of si ngle phase flow. Also at low vapor quality, the vapor volume fr action increases significantly w ith small changes in quality, and the gravitational cha nge in pressure drop is smaller at higher volume fraction. Thus a decrease in pressure drop is pr edicted. Given the complex nature of twophase flow, the pressure drop prediction across the microchannel plate is cons idered to be in reas onable agreement with the measured data. Flow Rate (kg/s) 0.0000.0050.0100.0150.0200.0250.030 Evaporator Pressure Drop (Pa) 1000 1500 2000 2500 3000 3500 4000 1.33 1.15 0.97 0.79 Measured Predicted H (m) Figure 44. Pressure drop vari ation across the microchannel eva porator plate with respect to mass flow rate. The reported uncertainty associated with the twophase frictional pressure drop computed with the MuellerSteinhagen and Heck (1986) correlation is approximately while the uncertainty associated w ith the void fraction com puted with the ZuberFi ndlay (1965) correlation is approximately. In order to determine the uncerta inty bounds of the m ass flow rate predictions, the mass flow rate is computed using th e upper and lower bounds of the frictional 41% 2.5% 77 PAGE 80 pressure drop and vapor volum e fraction. Figure 45 shows the upper and lower bounds of the predicted flow rates for H = 1.33 m and H = 1.15 m. The measured flow rates are compared with those predicted bounds. Only two c ondenser elevations are presen ted for discussion purposes. As shown, the measured flow rates generally fall within the predicte d bounds, thus providing further confirmation that the agreement between the measured flow rates and those predicted is reasonable. Heat Flux (kW/m2) 01020304050607 0 Flow Rate (kg/s) 0.00 0.01 0.02 0.03 0.04 1.33 0.79 Measured Error Bounds H(m) Figure 45. Upper and lower bounds for predicte d mass flow rates for H = 1.33 m and H = 0.79 m. Table 41. Experimentally observed limiting heat flux. Condenser Height, H (m) Limiting Heat Flux, 2 limqkWm 0.79 61.91 0.97 47.98 1.15 55.20 1.33 61.91 The heat fluxes shown in Figure 43 range fr om zero to a maximum value, at which point the system goes unstable. The system limiting heat flux for each condenser height is tabulated in 78 PAGE 81 T able 41. As observed, the limiting heat fl ux increases with increasing elevation of the condenser above the microcha nnel evaporator plate. 1000 1500 2000 2500 3000 3500 4000 0 0.01 0.02 0.03 0.04 0.05 Mean = 2510 Pa Standard Deviation = 287 Pa Skewness = 0.018 Kurtosis = 3.1 Percent Fluctuation = 11 % Standard Deviation = 282 Pa Skewness = 0.025 Kurtosis = 3.1 Percent Fluctuation = 12 %Probability Density FunctionPressure Drop (Pa) Mean = 2333 Pa H = 1.33 m (q" = 58 kW/m2) H = 1.15 m (q" = 50 kW/m2) Figure 46. Measured probability density functi on for pressure drop prior to unstable flow. 0.01 0.012 0.014 0.016 0.018 0.02 0 0.01 0.02 0.03 0.04 0.05 0.06 Mean = 0.016 kg/s Standard Deviation = 6.8x104 kg/s Skewness = 0.039 Kurtosis = 2.9 Percent Fluctuation = 4.1 %Probab ility Density FunctionFlow Rate (kg/s) Mean = 0.014 kg/s Standard Deviation = 8.3x104 kg/s Skewness = 0.099 Kurtosis = 3 Percent Fluctuation = 6 % H = 1.33 m (q = 58 kW/m2) H = 1.15 m (q = 58 kW/m2) Figure 47. Measured probability density functi on for mass flow rate prior to unstable flow. 79 PAGE 82 0.03 0.02 0.01 0 0.01 0.02 0.03 0.04 0.05 0 0.01 0.02 0.03 0.04 0.05 0.06 Mean = 0.011 kg/s Standard Deviation = 0.011 kg/s Skewness = 0.093 Kurtosis = 1.6 Percent Fluctuation = 93 %Probability Density FunctionFlow Rate (kg/s) Mean = 0.011 kg/s Standard Deviation = 0.012 kg/s Skewness = 0.19 Kurtosis = 1.5 Percent Fluctuation = 113 % H = 1.33 m (q = 62 kW/m2) H = 1.15 m (q = 56 kW/m2) Figure 48. Measured probability density functi on for mass flow rate following the onset of unstable flow. While the quasisteady mass flow rate and pres sure drop are very important variables for understanding the response of the system to st ep increases in heat flux, apparent random variations about the quasi steady variables are eq ually important. It is the variation from the quasisteady state that can send the twophase thermosyphon into an uns table operating mode. Therefore, the stochastic variations in pressure drop and flow rate are examined. Figures 36 & 47 show the probability density functions for bo th the pressure drop and flow rate immediately prior to encountering the flow instability for both H = 1.33 and H = 1.15. In each case the skewness is almost zero and the kur tosis is near to three, sugges ting the distribution is close to normal. The percent fluctuation, defined here as the ratio of the standard deviation to the mean, is 11% for the pressure drop and 4.1% for the flow rate for H = 1.33 m. For H = 1.15 m., the percent fluctuation is 12% fo r the pressure drop and 6% for the flow rate. The percent fluctuations are not statistically different be tween the two condenser elevations. At each condenser elevation considered in this study, the flow regime within the riser is observed to be annular immediately prior to instability. Knowing that bubble coalescence and void fraction play 80 PAGE 83 a critical role in the fluctuati on, it is quite reasonable that the fluctuations should not vary significan tly between condenser elevations as observed in Figures 46 and 47. Once the heat flux is increased beyond the lim iting heat flux, large s cale fluctuations are observed in the flow field. Figur es 48 & 49 show the probability density functions for pressure drop and flow rate for H = 1.33 m and H = 1.15 m slightly beyond the limiting heat flux. In each case the distributions resemble those for a sinusoi dal signal in the time domain. There are two distinct peaks for both th e pressure drop and flow rate. The flow rate is observed to go negative which is indicative of flow reversal. The averag e measured flow rate is 0.011 kg/s for each case. The flow rate has a positive skewness while the pr essure drop has a negative skewness. This is most likely a result of the dependency of pressure drop on flow rate within the system. As the flow rate is increased the vapor quality decreases and vice versa. Since pressure drop is strongly dependent on vapor quality, the pressure drop is la rgest at low flow rates where vapor quality is largest. Therefore, it is expected that a pos itive skewness in the flow rate would result in a negative skewness in pressure drop as observed in Figures 48 & 49. 1000 0 1000 2000 3000 4000 5000 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Mean = 2675 Pa Standard Deviation = 958 Pa Skewness = 0.5 Kurtosis = 2.4 Percent Fluctuation = 36 % Mean = 2475 Pa Standard Deviation = 1075 Pa Skewness = 0.37 Kurtosis = 2.1 Percent Fluctuation = 43 %Probability Density FunctionPressure Drop (Pa) H = 1.33 m (q = 62 kW/m2) H = 1.15 m (q = 56 kW/m2) Figure 49. Measured probability density func tion for pressure drop following the onset of unstable flow. 81 PAGE 84 A Fast Fourier Transform (FFT) wa s applied to the pressure drop and flow rate data taken at H = 1.33 m and H = 1.15 m. The data were taken for N = 512 data points with a sampling frequency, s f of 14 Hz, resulting in a period, T, of 37 s. Fig. 11 shows the power spectrum for mass flow rate, and the dominant frequency is approximately 0.8 HZ and 1 Hz, respectively, for H=1.33m and 1.15 m. The measured power sp ectrum indicates a very strong periodic low frequency fluctuation when the flow becomes unstable. 0 1 2 3 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 Power Spectrum (kg/s)2frequency (Hz) 0 1 2 3 0 0.005 0.01 0.015 0.02 Power Spectrum (kg/s)2frequency (Hz) H = 1.33 (q" = 62 kW/m2) H = 1.15 (q" = 56 kW/m2) Figure 410. Power spectrum for mass flow rate; (a) H = 1.33 m, (b) H = 1.15 m. Instability Prediction In order to reliably predict the onset of instability, it is necessary to identify the controlling mechanism leading to th e instability excursion. After giving consideration to the instability mechanisms described by Boure et al. (1973) and Tadris t (2006) it is concluded that the onset of the heat flux limiting instability obser ved in the thermosyphon facility is due to a static Ledinegg instability. Such an instability involves a sudden change in the flow rate to a lower value followed by oscillations and ultimately system failure. In forced convection, the phenomenon is governed by the slope of pressure gr adient with respect to mass flux for both the internal system characteristics (friction, acceleration, and gravity) and the external system 82 PAGE 85 characteristics (the pum ping mechanism). In the case of natural circulation, the pumping mechanism itself is the gravitational pressure he ad exerted by the downcomer. The criterion for unconditionally stable flow is simply, 0sP G (43) For the current facility, s Pis the change in pressure between the condenser outlet and inlet and is zero for quasisteady state operation. In orde r to investigate whether the system operation is stable, s P Gis first computed at each heat flux interval using Eq. 2.2. The quasisteady state mass flux is that which gives zero pressure drop between the cond enser inlet and exit and is determined for each heat flux interval. A pe rturbation in mass flux from the quasisteady state,, of 5% is used based on the percent fluctuati ons revealed experimentally in Figure 47. The change in system pressure drop resul ting from a 5% perturbation in mass flux is approximated as follows, ss s GGGPP P GG (44) The system heat flux for which s P G is negative is associated with unstable flow, and thus the limiting heat flux is that for which 0sP G For the flow to transition to the unstable regi me, the sum of the frictional, gravitational, and accelerational pressure change with respect to mass flux must be less than zero. Throughout the thermosyphon loop, the accelerational component of pressure drop is an order of magnitude smaller than the frictional and gravitational comp onents and can be considered negligible here for discussion purposes. At high heat flux, wher e the instability occurs, the frictional pressure gradient within the microchannel plate is approximately 50 times the pressure gradient in the 83 PAGE 86 liquid downcom er section of the thermosyphon loop and 2.5 times greater than that in the twophase riser section. Therefor e the pressure drop across the microchannel plate, which is dominated by friction, gives a reasonable approxi mation of the shape of the frictional pressure drop with respect to mass flux throughout the syst em and is shown in Figure 411 for illustration purposes. Flow Rate (kg/s) 0.0000.0050.0100.0150.0200.0250.030 Evaporator Pressure Drop (Pa) 1000 1500 2000 2500 3000 3500 4000 Measured Predicted H = 1.33 m Unstable (Region II) Stable (Region I) Figure 411. Pressure drop across the microchanne l plate for varying flow rate at H=1.33 m. In regionI, where the slope of the pressure drop as a function of flow rate curve is positive, a small perturbation decreasing the flow rate will result in increasing vapor quality and larger gravitational pumping poten tial. Simultaneously, the frictional pressure drop decreases due to decreased flow rate, and thus the flow rate w ill adjust back to its original value. Similarly, in regionI, the flow rate will return to its original value following a perturbation increasing the flow rate. Thus the flow is stable in regionI. In contrast, in region II the slope of the pressure drop as a function of flow rate curve is negative, and a small perturbation increasing the flow rate will result in decreasing vapor quality and gr avitational pumping potential. Simultaneously, the 84 PAGE 87 frictional pressure drop is reduced, and further in creas es in flow rate and decrease in vapor quality can result, until the flow swings back to region I. However, the flow state does not remain in region I for long due to the high heat flux and rapid increase in vapor quality. This behavior continues back and forth in an oscillato ry mode until eventual dryout. It is noted that the flow in region II may be stable as long Equation 43 is satisfied. 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 35 40 45 50 55 60 65 70 Condenser Height (m)Limiting Heat Flux (kW/m2) Experimental Predicted Upper & Lower Bounds Figure 412. Comparison of measured and predicted heat flux at the onset of instability. Figure 412 shows the computed limiting heat fluxes compared with those determined experimentally. Also shown are the upper a nd lower bounds of the computed limiting heat fluxes based on the uncertainties associated with the frictional pressure drop and vapor volume fraction predictions. To compute the upper and lower bounds, the model was run using the upper and lower bounds of both the frictional pressure drop and the vapor volume fraction, based on the associated error in the correlations, and the highest and lowest predicted limiting heat flux were taken as the bounds. The experimental m easurements tend to fall on the lower bound of uncertainty. Given the complex nature of twopha se flow, the agreement is considered to be satisfactory. 85 PAGE 88 CHAP TER 5 PERFORMANCE OF ALUMINUM AND CA RBON FOAMS FOR AIR SIDE HEAT TRANSFER AUGMENTATION Background and Objectives Improving condenser and evaporator performanc e has been a widely studied topic in the context of compact and efficient he at exchanger design. It has been a subject of interest in the process, cryogenic, power, refrigeration and ai r conditioning industries am ong others. While air flow heat exchangers have been in use by i ndustry for many years, new challenges arise and opportunities for performance enhancement persist. Since most air cooled heat exchangers are limited in performance by the thermal resistance on the air side, focus has been on advancing the air side capabilities. This is achieved by increasing the heat transfer surface area in contact with the flowing air, or by increasing th e heat transfer coefficient. In many applications it is common to use a louvered fin that act s to break up the thermal boundary layer and improve net heat transfer with an increase in bot h surface area and heat transf er coefficient (Chang and Wang, 1997). While densely packed louvered fins will greatly enhance heat tr ansfer, the increased pressure drop can make such a configuration undesirable due to incr eased pumping power. Typically there is an optimum sp acing that represents a tradeoff between enhanced heat transfer and increased pumping power requirements. Another approach to enhancing the air side heat transfer rate is to replace the fins with a high thermal conductivity porous structure. These structures are capable of increasing the heat exchange surface area while producing a tortuous flow path that breaks up the thermal boundary layer and enhances the heat transfer coefficien t. However, the pressu re drop across the porous media is amplified as a result, and increased pumping power is required. Numerous experiments have considered forced convective air heat transfer enhancement associated with porous metal foams. In general the experiments can be categorized as those that 86 PAGE 89 m easure a bulk heat transfer coefficient and thos e that measure the inters titial heat transfer coefficient. The bulk heat transf er coefficient is defined as the ratio of the surface heat flux to the difference in temperature between the heater wall and inlet air temperature. Convective air bulk heat transfer coefficient measurements for porous metal foams have been reported by Hsieh et al. (2004), Haak et al. (2001) Kim et al. (2000), and Tzeng ( 2007). While such measurements may be useful to test a specialized geometry and flow configuration, th eir utility is generally limited. The other approach is to trea t the metallic foam as flow through a porous media and solve the solid and fluid coupled energy eq uations to extract the interstitial heat transfer coefficients from the measured thermal fields. This appr oach has been pursued by Calmidi and Mahajan (2000), Jeng et al. (2005) and Hwang et al. (2002). The advantage to this approach is that the interstitial heat transfer coefficient only depends on the local interstitial geometry of the metallic foam. Thus, the results may be combined with a solution of the solid and fluid coupled energy equation to analyze porous foam heat exchangers with any arbitrary geometry. Analyses of this type are reported by Lu et al (2006) and Zhao et al. (2006). This paper considers the heat transfer augmentation provided by porous metal foams and compares them against conventional louvered fins. Three types of aluminum and three types of carbon foam are considered. An experimental f acility has been set up to measure the pressure drop and volumetric heat transfer coefficient with air flow veloci ties ranging from 1 to 6 m/s. The DarcyForcheimer equation is used to corr elate the pressure drop data. The Nusselt number is correlated with Reynolds number based on appropr iate length scales. Due to very significant differences in the porous structures, different correlations are used for the aluminum and carbon foams. For comparison against conventional louv ered fins, a hypothetical heat exchanger is 87 PAGE 90 analyzed an d the performance is determined based on a Coefficient of Performance (COP), Compactness Factor (CF), a nd Power Density (PD). Experimental Facility Three different aluminum foam samples manufactured by Energy Research and Generation Inc. are studied here along with th ree carbon foam samples manufactured by Koppers Inc. The samples are shown below in Figs. 51 a & b and the properties are given in Table 51. The aluminum foams are identified by the manufactur er according to pore density, with values of 10, 20, and 40 pores per inch (PPI). The bul k density of all aluminum foam samples is 216 3kgm. The pore sizes for the 10, 20, and 40 PPI aluminum foams are 0.5, 1, and 2 respectively. The three carbon foam samples are identified by the manufacturer as L1, D1, and L1A. The respective pore sizes according to the manufacturer are 600, 650, and 500mmm Due to the small pore size and large pressure drop required to convey air through the foams, they were modified by machining cylindrical air passages in the axial directi on. In each case, 80 passages, 6.7 apart, each with a diameter of 3.2 were machined into the foam providing for a uniform homogeneous geometry. The modified L1 D1 and L1A, foams have bulk densities of 317, 400, and 284mm mm 3kgm, respectively. Due to the large ratio of passage diameter to pore diameter, there is not a large pressure drop variation between carbon foam samples. It is assumed throughout the analysis that convection h eat transfer is dominated by air flow through the machined passages. Therefore the same surface area per unit volume, a, for each carbon foam sample is 5.24 and the porosity based on th e flow passages is 0.166. 1m 88 PAGE 91 (A) (B) Figure 51. A) Carbon foam samples, from t op to bottom: L1A, D1, & L1. B) Metal foam samples, from top to bottom: 10PPI, 20PPI, & 40 PPI. The experimental facility displayed in Fig. 52 is used to measure the pressure drop and heat transfer coefficient with ai r blowing through the foam. As shown, a variable speed blower forces air through the foam sample in which a 15.24 cm. x 15.24 cm. aluminum plate is bonded to the foam using Sbond materi al. A 15.24 cm. x 15.24 cm. foil heater is in contact with the bottom side of the aluminum plate. Five type E thermocouples are embedded into the aluminum 89 PAGE 92 plate to m easure the lower wall temperature whil e five more type E thermocouples are located axially along the top surface of the foam to m easure the upper wall temperature. The upper wall temperature is used in this analysis instead of the lower wall temperature in order to eliminate the error associated with contact resistance between the heater and the foam. The foam sample is 15.24 cm. x 15.24 cm. x 2.54 cm. in volume. Vali dyne differential pressure transducers measure the pressure drop across the foam. Seven pressure taps are locate d axially along the foam sample to measure the pressure varia tion along the channel. Another Validyne differential pressure transducer measures the pressure drop from a pito t tube to obtain the mean air velocity. All analog signals were captured using a digital data acquisition system cons isting of a Measurement and Computing CIOEXP32 multiplexer and PCIM DAS1602 16 bit analog to digital converter. A custom algorithm was developed to convert measured voltages to pressure, temperature, and velocity. The experimental facility is used to measur e the heat transfer coefficient and pressure drop for each foam sample at different velocity increments. The pressure drop is measured adiabatically. Throughout the heat transfer experiments the heat flux remains fixed at 9.77 2kWm, and the air velocity is increased in incr ements. The mean velocity ranges from approximately 1 to 6 ms for the pressure drop and heat transfer experiments, beyond which point there is a substantial increase in pressure drop and only a marginal increase in the heat transfer rate. At each velocity increment, the upper foam temperature is observed until steady state is reached, at which poi nt 500 data points are sampled at 1 Hz. The measurement is repeated 5 times to insure repeatability. The ai r velocity is then increas ed and the procedure is repeated. The uncertainty in temperature and pre ssure drop is explained in detail in Appendix B. 90 PAGE 93 Figure 52. Porous m edia experimental facility. Experimental Results For each foam sample considered, the meas ured pressure drop and foam upper wall temperature are tabulated in Tables C1 and C2 in AppendixC. The axial pressure variation is shown in Fig. 53 for the 10 PPI sample operating at =1.13mu ms. The observed linear variation is indicative of fully developed flow. dP dx, provided for each foam sample in TableC1, gives the pressure gradient. The measured pressure gradient is shown in Fig. 54 & 55 with varying air velocity. As observed, the pressure gradient increases parabol ically with increasing air velocity. The foam upper wall axial temperat ure profile is shown in Fig. 53 for the 10 PPI sample operating at =1.13mu ms. The wall temperature gradient sh own in Fig. 53 is used to determine the sensible heat transfer to the flui d, which is compared against the measured heat input confirm an energy balance. The upper wall temperature of the upper foam at x = 7.62 cm. is plotted versus mean fluid velocity in Fig. 56. As observed, the temperature drop decreases with increasing air velocity. The upper wall temp erature profiles of the three aluminum foams appear to be similar. At low velocity an increm ental increase in velocity results in a large drop 91 PAGE 94 in foam wall temperature and a moderate increase in pressure drop. Howeve r, at relatively high air velocities, an incremental change in velocity results in a large increase in pressure drop and only moderate decrease in foam temperature. Table 51. Foam properties. Foam Sample K 2m F C p p d m s ek W mK a 1m 40 PPI 6.98 9100.020 0.918 NA 5.08410 9.78 2760 20 PPI 1.21 8100.021 0.918 NA 1.02310 9.78 1770 10 PPI 1.98 8100.027 0.918 NA 2.03310 9.78 804 L1A 1.66 8100.034 0.166 .806 5410 48.6 5.24 D1 1.66 8100.034 0.166 .752 6.5410 97.2 5.24 L1 1.66 8100.034 0.166 .735 6410 61.8 5.24 The pressure gradient for carbon and aluminum foams may be correlated using the DarcyForcheimer equation (Vafai and Kim, 1989), 2 flflF mC dP uu dxK Km (51) where K is the permeability, is a dimensionless parameter that accounts for inertia effects, and is the mean velocity. and K are determined from th e tabulated pressure drop data shown in Table 41 by fitting the data to Eqn. (51) and extracting and K The measured K and are listed in Table 51 for each foam. In th e case of carb on foam, the difference in the pressure drop between data samples is not larg e because the majority of air flows through the machined passages. The small difference is attrib uted to flow through th e interstitial pores. FCFCmuFCFC 92 PAGE 95 Axial Location (m) 0.000.020.040.060.080.100.120.140.16 Foam Upper Wall Temperature (C) 50 55 60 65 70 75 Pressure (Pa) 10 20 30 40 50 60 70 80 Measured Temperature Measured Pressure Figure 53. Axial pressure upper wall temperature variation for 10 PPI aluminum foam sample, = 1.13mu ms. Mean Air Velocity (m/s) 012345 dP/dx (Pa/m) 0 1000 2000 3000 4000 5000 6000 D1 L1 L1A Figure 54. Carbon foam pressure gradie nt variation with mean air velocity. 93 PAGE 96 Mean Air Velocity (m/s) 012345 6 dP/dx (Pa/m) 0 2000 4000 6000 8000 10 PPI 20 PPI 40 PPI Figure 55. Aluminum foam pressure grad ient variation with mean air velocity. Mean Air Velocity (m/s) 01234567 Upper Wall Temperature at z = 7.62 cm (C) 20 30 40 50 60 70 80 90 100 10 PPI 20 PPI 40 PPI D1 L1 L1A Figure 56. Foam upper wall temperature variation with mean air velocity. 94 PAGE 97 Heat Transfer Analysis In order to extract the local heat transfer co efficient from the measured foam temperatures, the steady volumeaveraged momentum equatio n and a two equation nonequilibrium heat transfer model proposed by Calmidi and Mahajan (2000) are employed. In this analysis it is assu med that all thermophysical properties of the solid and fluid are independent of temperature. The foams are modeled as homogenous porous media, and the effects of radiation, natural c onvection, and thermal dispersion are negligible (Jeng, 2005). The governing conservation equations can be expressed as (Lu et al, 2006), Momentum 2 2 flFflC PVVVVVV K K (52) Energy (Solid) 21sssfflskThaTT 0 (53) Energy (Fluid) 2 f lpflflflflsfsflCVTkThaTT (54) where f l is the dynamic viscosity of the fluid, is the porosity, K is the permeability, s fhis the interstitial heat tran sfer coefficient, and is the surface area per unit volume for the porous m aterial. The product of a s fh and is the volumetric heat transfer coefficient, Here avh f l and are the respective density and specific heat of the fluid, while ,pflC f lkand s k are the thermal conductivity of the fluid and solid, respectively. During the experiments, the upper wall of the fo am test section is in sulated. The boundary conditions are, 00uy0uyH 0 0f s s ef e y yT T kk yy w q 95 PAGE 98 0f s sefe yH yHT T kk yy 0fTxT 2 2 00f xT x 2 2 00s xT x and 2 20s xLT x The following dimensionless variables are introduced (Calmidi and Mahajan, 2000)(Lu et al., 2006) ws eTT qLk x x L 'y y L mu U u 24KK Da H, 322m F f luH C 2 sf shaL Bi k ReL f luL Rem K fluK, Pr f l f l and sf f lpflmhaL St Cu where and are the respective local and mean velocities, is the applied heat flux along the lower wall, and umuwq f ekkfl and 1 s eksk are the effective thermal conductivities of the respective fluid and solid. Conservation equations (52 through 54) are nondimensionalized as 2 2 2Re 11 'K F K KU UCU yDa Da (55) 22 220 ''ss flsBi xy (56) 22 2211 'RePr''fl flfl s fl LSt U xxy (57) with boundary conditions, '00 Uy '0 UyHL '0 '01 ''fef s se y yk yky '0'0flx''fefl s se yHL yHLk yky 0 2 2 '00 'fl xx 2 2 '10 'fl xx and 2 2 '10 's xx The exact solution to the momentum e quation has been provided by Vafai and Kim (1989) and is given below, 21'c mu AB UsechDy uA C (58) 122 3KA Da (59) 96 PAGE 99 1214 3KKB DaDa (510) 11 1/2 A Csech DAB (511) 2 AB D (512) 0.5 0.4 0.3 0.2 0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 y/HU uc/um = 1.007 = 1113.7 DaK = 1.22 x 106 Figure 57. Dimensionless velocity profile. 97 PAGE 100 (a) (b) Figure 58. Dimensionless temperatur e profile. a) solid. b) fluid. 98 PAGE 101 In order to com pute the veloci ty field, a guess is made for cmuu The velocity profile in then computed from Eq. (58), and a check is made to insure mass conservation, 1 1AUdA A If mass conservation is satisfied the solution is complete, otherwise a new guess is made for cmuu and the procedure is repeated until the integral is satisfied. Mean Air Velocity (m/s) 123456 Volumetric Heat Transfer Coefficient (W/m3K) 0 10000 20000 30000 40000 50000 60000 D1 L1 L1A Figure 59. Carbon foam volumetric heat transfer coefficient variat ion with mean fluid velocity. Solutions to the energy conser vation equations were obtained using the Thomas algorithm. For 1.007cmuu 1113.7 5 KDa=1.2210 Re7,3048 Pr0.71 0.54St 0.92 and the computed dimensionless velocity prof ile is shown in Fig. 57 while that for tem perature is shown in Figs. 58 a & b. 7.56Bi As observed, the velocity is uniform along th e centerline with a valu e close to the mean velocity. The solid foam temperature decrease s monotonically in the ydirection and increases 99 PAGE 102 linearly in th e xdirection. The fluid temperat ure initially decreases in the ydirection while increasing in the xdirection. Mean Air Velocity (m/s) 0123456 Volumetric Heat Transfer Coefficient (W/m3K) 0 5e+4 1e+5 2e+5 2e+5 3e+5 3e+5 10 PPI 20 PPI 40 PPI Figure 510. Aluminum foam volumetric heat tran sfer coefficient variation with mean fluid velocity. When solving Eqs. 56 and 57 the heat tr ansfer coefficient is not known a priori. Therefore, a guess is made for the heat transfer coefficient, s fh and solutions for Eqs. (56 & 57) yield s and f l The upper foam temperature measured at 2 x L is compared to that computed using the guessed value for s fh If agreement is achieved, the guessed heat transfer coefficient is correct. If they do not agree, a new guess is ma de for the heat transfer coefficient and the computation is repeated until agreement is achieved. In cases where the mean air velocity is large, the solid foam and air are in near ther mal equilibrium. Thus large errors arise when extracting the volumetric heat tr ansfer coefficients since most thermal transport is due to conduction through the porous medium. In such cases, a small change in the measured foam 100 PAGE 103 upper wall tem perature can result in a large change in the computed volumetric heat transfer coefficients. Results The volumetric heat transfer coefficients, vsfhh a are computed for each of the six foams and are summarized in Table C2 in the Appendix. Their variations with mean air velocity are shown in Figs. 59 & 510. The aluminum foam s have the largest volumetric heat transfer coefficient. The aluminum foams are capable of handling a large heat load because they have large specific surface areas and large effective ther mal conductivities. In order to correlate the volumetric heat transfer data, the perm eability based Nusselt number is introduced, v K flhKa Nu k. Fig. 511 shows the Nusselt number as a function of the Reynolds number, Rem K fluK, for the three aluminum foams. It is observed that the data are well collapsed for relatively low air velocity. At high air veloci ty, the solid and fluid are in near thermal equilibrium, and a large change in heat transfer coefficient results in only a small change in the solid foam temperature. Therefore, the error in the correlation at high velocities has a minimal effect on the computed solid temperature profile. Overall, the correlation allows for an accurate prediction of wall temperature for the full range of air velocities investigated here. The aluminum foam Nusselt number may be expressed as, Re 0.01 4887p K Kd Nu K (513) For the carbon foam samples, if it is assumed th at all of the air passes through the machined channels, and ideally each foam sample shoul d have the same volumetric heat transfer coefficient at the same Reynolds number. However, variations in the data are observed. These 101 PAGE 104 are explained by the fact that so m e air passes through the smaller in terstitial pores of the foam. Thus, it is reasonable to expect th at the data can be correlated us ing the pore diameter. Fig. 512 shows the product of the Nusselt number, 2 v D f lhD Nu k and 3/2 pD d as a function of Reynolds number, Rem D f luD for the three carbon foams. As observe d, the data are also well correlated, and the carbon foam Nusselt number may be expressed as 3/2 1.4Re 511108 90D D pD Nu d (514) The standard error between the measured and predicted pressure drops and volumetric heat transfer coefficient for the aluminum and car bon foams are summarized in Table 52. As observed, the volumetric heat transfer coeffici ent is predicted with reasonable accuracy. ReK 050100150200250300350400 NuK(dp/K1/2)1 0.02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 40 PPI 20 PPI 10 PPI Eq.(13) Thermal Equilibrium Figure 511. Nusselt number corre lation for aluminum foams. 102 PAGE 105 ReD 20004000600080001000012000 NuD(dp/D)1.5 0 100 200 300 400 500 600 D1 L1 L1A Eq.(14) Figure 512. Nusselt number co rrelation for carbon foams. Table 52. Relative error for pressu re drop and Nusselt number correlations. Foam Sample Relative Error P Relative Error kNu 10 PPI 4% 5% 20 PPI 3% 19% 40 PPI 5% 5% L1A 2% 4% D1 9% 5% L1 9% 9% Hypothetical Heat Exchanger Performance For comparison purposes, a hypothetical heat ex changer having specified dimensions is analyzed for each foam sample as well as for louvered fins. An illustration of the geometry is provided in Fig. 513. The louvered fin geometry shown in Fig. 513a is common in compact heat exchangers and can be further classified as multilouvered. The flat tube configuration has been chosen because it allows for enhanced performance. Fiebig et al. (1994) shows that the heat exchanger element with flat tubes gives nearly twice as much heat transfer and only half as much pressure 103 PAGE 106 loss as the corresponding heat ex changer elem ent with round tubes allowing for a more compact, efficient design. Figure 513. Hypothetical heat ex changer in cross flow. a) louve red fin. b) carbon/metal foam configurations. The foam geometry shown in Fig. 413b consti tutes the same flat tube arrangement as the louvered fin case and is oriented su ch that the heat exchanger worki ng fluid is in cross flow with the ambient air, as is the case with the louvered fins. Since the only difference in the two designs is the extended surface area be tween tubes, the configurati on allows for a quantitative comparison between each foam sample and the louvered fins. 104 PAGE 107 To determ ine which configuration or enhanced surface yields the best, it is important to first have a clear definition of what constitutes good performance. The ideal heat exchanger will be compact in size and capable of removing the required heat load, while consuming a minimal amount of electric energy from the blower. Sinc e a more compact heat exchanger must often consume more energy than a less compact heat exchanger under the same heat load, it is important to base the performance on both efficien cy and compactness so that the designer can choose an optimal configuration that best suits each specific application. Th is paper evaluates the performance of each sample based on a Coefficien t of Performance (COP; defined as the ratio of the total heat removed to the electrical input of the blower), Compactness F actor (CF; defined as the total heat removed per unit volume), and Po wer Density (PD; defined as the total heat removed per unit mass). Numerical Model In order to predict the heat removal capabili ties of the hypothetical heat exchanger shown in Fig. 14b, the momentum equation for the fluid as well as the coupled energy equations for the solid and fluid are solved using Eqs. 52 thr ough 54. The equations can be normalized by the following scales, wTT TT x x L y y L mu U u and reduced to the following for twodimensional rectangul ar coordinates, 2 2 2Re 11 'K F KU UCU yDa Da (515) 22 220 ''ss flsBi xy (516) 22 2211 'RePr''fl flfl s flSt U xxy (517) 105 PAGE 108 Whe re 24KK Da H 322m F f luH C 2 sf shaL Bi k ReL f luL Pr f l f l and sf f pflmhaL St Cu The hypothetical heat exchanger described here will have a constant wall temperature of 100 and an ambient air stream entering at 25 By implementing a constant wall tem perature, the hypothetical heat exchanger wi ll emulate a configura tion with zero thermal resistance between the tubular inner working fluid and the tube wall. Assuming a noslip condition at the wall, the boundary conditions along th e tube wall and the fluid at the inlet are as follows, Uy, CC'00Uy'/HL 0 '0','01fl sxyxy', '/syHL','/1xyHL ',xfl '0,'0flxy 2 2 '10 'fl xx and 2 2 '10 's xx To simplify this analysis, two tubes are to be considered, with heat transf er taking place along both boundaries of the foam. Using the volumetric heat transfer coeffici ent taken from Eqs. 513 through 514, along with the permeability, K, and the inertia coefficient, taken from Table 51, the governing equation s are solved and the temperature field is obtained. FC The total heat removed, is computed from removedQ ,,, 00 0 021LL f s removedpffoutfinsf y yT T QmCTTWkdzWkd yy zP, (518) while the pressure gradient is evaluated from Eq (51). The electrical in put to the blower is estimated as the product of th e volumetric flow rate of air and the pressure drop: inputmcEuA (519) where is the frontal area of the heat exchanger s een by the ambient air stream. Finally the COP, CF, and PD are evaluated: cA 106 PAGE 109 removed inputQ COP E (520) removedQ CF V (521) removedQ PD m (522) Where V is simply the volume of the hypotheti cal heat exchanger and m is the mass. Figure 514. Geometric comparis on of louvered fin and foam h eat exchanger configurations. For proper comparison, the heat removal capab ilities of the louvered fins are to be obtained as well. Since there ha ve been numerous experiments c onducted on this subject, several correlations exist in the literature that provide ac curate approximations for the interstitial heat transfer coefficient and pre ssure drop. Dong et al. (2007) have provided a correlation for interstitial heat transfer coefficient as well as pressure drop based on a large data base for multilouvered fins. Table 53. Geometric parameters for louvered fin configuration Fh (mm) 60 Fp (mm) 4 Lh (mm) 58 Lp (mm) 1.2 Ld (mm) 100 La (degrees) 22 107 PAGE 110 Using the correlations provided by Dong et al (2007), the COP, CF, and PD are computed for the louvered fin heat exchangers and comp ared against those for each foam sample at different air velocities. The geometric parame ters chosen to describe the louvered fin configuration are pr ovided in Table 53 and illustrated in Fig. 514. The bulk density of the louvered fins is estimated based on these geometric parameters to be 128 3kgm. The results of the comparison are illustrated belo w in Figs. 515 through 517. Mean Air Velocity (m/s) 01234567 COP 0 200 400 600 800 1000 1200 1400 10 PPI 20 PPI 40 PPI D1 L1 L1A Louveredfin Figure 515. Comparison of coefficient of perfor mance for louvered fin and foam configurations Fig. 515, shows that the multilouvered fin configuration is favorable over both the carbon and aluminum foam samples when COP is th e main consideration. This is due to the added pressure drop resulting from the tortuous flow path characteristic of the foam. In the case of a louvered fin, several indepe ndent length scales ex ist that have an effect on COP including the fin pitch and louver pitch. A small fin pitch will increase surface area, while a small louver pitch will break up the thermal boundary layer more frequently. By reducing both of these length 108 PAGE 111 scales, the heat rem oval should increase due to th e increase in both surface area and heat transfer coefficient. However, this will also constrict th e flow path and induce mixing, giving rise to an increase in pressure drop. Mean Air Velocity (m/s) 01234567 CF (W/m3) 0.0 200.0x103400.0x103600.0x103800.0x1031.0x1061.2x1061.4x1061.6x1061.8x106 Aluminum Foams D1 L1 L1A Louvered Fin Figure 516. Comparison of comp actness factor for louvered fi n and foam configurations In the case of flow through a homogeneous por ous structure such as aluminum foam, the length scales of importance are the pore diam eter, and the ligament diameter, or more importantly the square root of permeability. A more fundamental study on the flow structure within porous media is needed in order deve lop a designer foam, capable of producing COP comparable to conventional louvered fin conf igurations. A comprehensive study using the LatticeBoltzmann method may provide insight into the optimal geometric length scales needed to improve overall foam performance, as the cu rrent foam structures do not provide adequate COP when compared to louvered fins. 109 PAGE 112 Mean Air Velocity (m/s) 01234567 PD (W/kg) 0 2000 4000 6000 8000 10000 Aluminum Foams D1 L1 L1A Louvered Fin Figure 517. Comparison of power density fo r louvered fin and foam configurations Fig. 516 shows that the CF is highest for the D1 carbon foam sample. This is a result of several factors. One is the high thermal conductivity of the carbon fo am that helps to increase the fin efficiency. Another is due to finite leakage through the pores, whic h provides a large surface area in contact with the air and a tortuous flow path that induces mixing. As seen in Fig. 517 however, the power density is largest for the louvered fins. This is due to the relatively low density of the louvered fin configurat ion. It is interesti ng to note that the CF for the louvered fin configuration drops off more rapidly than the fo am samples at low air velocity. This could be because at low fluid velocity, or low Reynolds number, the thermal boundary layer will be largest, at which point it is advantageous to break up the boundary layer frequently. At high Reynolds number however, the boundary layer will be smaller, and the effects of the added mixing will be marginal. Therefore it should be e xpected that a foam structure with a tortuous flow path will provide better heat transfer ra tes at low Reynolds number than a louvered fin configuration that is less effective in breaking up the thermal boundary layer. 110 PAGE 113 CHAP TER 6 CONCLUSIONS The work presented here encompasses the development of an advanced high power density thermal management system for PEM fuel cells that significantly increases the heat removal capacity while decreasing the required pumping power compared with other state of the art techniques. Analytical design tools have been developed for twophase thermosyphon thermal management systems, which should prove useful to designers to optimize for power and energy density. An extensive investigation into the pe rformance of a microchannel evaporator plate has contributed significantly to th e overall understanding of the th ermal hydraulics encountered in small channels. It is envisioned that this work will prove very useful to NASA for the development and implementation of innovative thermal management systems suitable for space exploration. A microchannel evaporative c ooling plate thermosyphon facil ity has been developed to remove waste heat from PEM fuel cell stacks. The overall performance of th e evaporator plate is excellent. The plate operates in a natural circulation twophase thermosyphon, and the thermal capacity for the current microchannel plate varies between 39 and 62 kW/m2, depending on the condenser height, which far exceeds the heat flux currently being rejected from PEM fuel cell stacks. The plate wall temp eratures range from 6682 C, depending on heat flux and position on the plate. A thermal hydraulic model has been pr esented which gives a satisfactory prediction of the mass flow rate, pressure drop, and microchann el plate thermal field. Similarly, the quasisteady behavior of a twophase microchannel th ermosyphon has been investigated to obtain a quantitative understanding of th e onset of unstable flow phenomenon observed at large heat fluxes. Static instabilities are largely influen ced by the height of the condenser and stochastic variations in flow rate. A onedimensional momentum equationbased model has been presented 111 PAGE 114 to provide a fra mework for pred icting the onset of the instability. Reasonable agreement is obtained between the observed and predicted inst ability limiting heat flux. The overall thermal hydraulic model describing the microchannel evaporat or plate is found to be suitable for design and analysis applications. In order to improve overall power /energy density of a PEM fuel cell, it is beneficial to integrate an air cooled condenser with a large COP, CF, and PD within the thermosyphon. This work evaluated the air side heat transfer performance of thr ee aluminum foam and three carbon foam samples to gain insight into the heat remo val capabilities and pressure drop requirements of the various foam samples and to compare against conventional louvered fin configurations. Correlations for computing the pressure drop and volumetric heat transfer coefficient based on the geometric properties of the fo am have been provided. To be tter quantify the heat transfer performance of these foam samples, comparison with conventional louvere d fin geometries has been made using a COP, CF, and PD as the ba sis for comparison. Resu lts show the carbon foam samples provide significant improve ment in CF but the COP and PD are considerably lower. The highest COP and PD are achieved with the mu ltilouvered fin configuration; however the resulting CF is substantially lower than th at for carbon foam. In conclusion, carbon foam performs better than conventional platefin heat exchangers in a pplications where a reduced heat exchanger volume is favored over reduced mass and/or COP, but in cases where reduced mass and/or enhanced COP are important, a multilouver ed fin type configuration proves to be more advantageous. In conclusion, these correlations along with the thermal hydraulic model of the evaporator plate are expected to contribute to further opt imization of the microchannel thermosyphon with respect to bo th power and energy density. 112 PAGE 115 APPENDIX A INSTABILITY IDENTIFICATION In order to a ccurately predict the onset of inst ability, it is first important to identify the mechanism causing the phenomenon. By rulin g out various instab ilities one by one, identification can be achieved. Using the cl assification of Boure et al. (1973), several mechanisms can be ruled out. For instance th e boiling crisis has been observed to occur simultaneously with flow oscillations (Boure et al. 1973) and does not occur prior to any instability. Since oscillations in flow rate them selves result in oscillati ons in vapor quality, the critical heat flux (CHF), may occur at a lower heat flux than would o ccur under constant flow rate conditions. In general, the boiling crisis is an unfavorable consequence of flow instability and therefore is not in itself a mechanism for the onset of instability. However, special consideration must be made when consideri ng twophase flow in narrow spaces. Since the channels under investigation have a 1 mm x 1 mm cross section, it possible that the CHF occurs at a lower heat flux than found in larger channels. Brutin & Tadrist (2004) show that the instability does not occur when the slope of pressu re drop versus flow rate is negative as is the case in the current investigation. Flow pattern transition occurs near the transition from bubbl y to slug flow. Since the pressure drop in annular flow is less than that for bubbly flow, the tr ansition can result in a sudden increase in flow rate. This in turn will cause a decrease in void fraction and the flow will revert back to a bubbly regime and the cycle is repeated. This phenomenon is not observed in the current experiments. The rise r section of the system is obser ved to be in the annular flow regime prior to instability, and does not osci llate between regimes until after the flow goes unstable. Another possible mechanism for instabili ty is bumping, chugging, or geysering. These phenomenons are a result of a lack of nucleation s ites. Under these conditions, large quantities 113 PAGE 116 of superheat build up. At this point, nucleation sites appear a nd the wall tem perature is reduced. If the working fluid has a good wettability, the nucleation sites will diminish and the cycle repeats. This frequent change in the number of nucleation si tes causes a fluctuation in vapor quality, and in turn causes oscillations in flow rate. This can be ruled out in the current experiment since vigorous boiling is observed at h eat fluxes relatively far below the heat flux at the onset of instability. Acoustic oscillations are a re sult of the propagation of pres sure disturbances along the heated channel. The frequencies of these osci llations are generally large, and the amplitude small (Boure et al. 1973) compared to those found in the current experiment. These low amplitude high frequency fluctuations may be occu rring prior to the large scale instabilities but are of no real interest here. This paper focuses on large scale, system in stabilities that can cause a boiling crisis and eventually system failure. Density wave oscillations are explained by propagation time lags describing the flow rate, vapor generation rate, and pressure drop within the system. A more de tailed description of density wave oscillations is described by Yadi garoglu and Bergles (1981). Kakac and Bon (2008) provide an instability map, and show that density wave oscillations occur when the slope of pressure drop versus flow rate is positive. Th e current investigation shows that the instability arises when the slope of pressure drop versus flow rate is negative. Therefore, the density wave oscillations can be ruled out. Thermal oscillations have been identified by Stenning and Veziroglu (1967) to be associated with the thermal re sponse of the heated wall afte r dryout much like the CHF type instabilities found in sm all channels. Boure (1973), Kakac & Bon (2008), have shown that this type of oscillation requires an interaction with high frequency density wave oscillations that act 114 PAGE 117 115 as disturbances to destabilize film boiling. Sin ce density wave oscillations are not present in the current experiments, thermal oscillations can be ruled out as well. BWR type instabilities are specific to nuclear reactors and involve the coupling of fuel dynami cs. It is obvious that these phenomena do not occur in the thermosyphon facility used here for testing. Finally, parallel channel instabilities occur when th e test section is comprised of tw o or more parallel channels. The test section used in this experiment is not susceptible to this phenomenon since each channel is connected. Due to the fabrication process us ed to manufacture the test section to allow for flow visualization, a small amount of fluid is able to flow between channels, eliminating the possibility of any parallel channel instability. After ruling out several mechanisms for the onset of instability, the remaining possibilities include Ledinegg in stability, pressure drop oscillations, and natural circulation instability. Both pressure drop oscillations and natural circulation instab ility are first triggered by the static instability (Ledinegg) before becoming dynamic. Since both pressure drop oscillations and natural circulation instability are a direct result of Ledinegg instability, the onset of these types of instability can be determined by successfully predicti ng the onset of Ledinegg instability. However, it is reas onable to assume that the observe d dynamics are not largely due to pressure drop oscillations since the compressible volume is ra ther small. Throughout the loop there exists a compressible volume in both the condenser and the riser. Therefore, in order for the vapor to compress within the condenser an equal amount of vapor will have to expand within the riser in order to obtain a mass balance. Therefore even for the case where the density fluctuations within the riser a nd condenser are perfectly out of phase, the overall compressibility of the system is limited by the compressibility within the riser. PAGE 118 APPENDIX B UNCE RTAINTY AND CALIBRATION In order to gage the accuracy of the various measurements, such as temperature, pressure drop, flow rate, heat flux, and vapor quality, the uncertainty in these parameters must be estimated. Uncertainties exist in the various devices, such as thermocouples, pressure transducers, volt meters and amp meters, along with uncertainties from the data acquisition board and calibration curves. This section focuses on both the calibrati on procedures as well as the associated errors relating to each measurement. Data Acquisition (DAQ) The data acquisition signal is read using a Measuremen t and Computing CIOEXP32 multiplexer and a PCIMDAS1602/16 bit analog to digital (AD) converter. For each measurement 500 samples are taken at 1 Hz, the data are averaged, and voltage is recorded. The measurement is repeated 3 times to insure repeatability. When measuring voltage for thermocouples, the voltage input is on the or der of mV. The multiple xer takes the millivolt signal and applies an appr opriate gain before it is sent to the AD board. The absolute error associated with this process is estimated by the manufacturer to be 0.2%MPV The maximum voltage signal recorded after the gain is applied for the range of temperatures measured in these experiments is 3 V. Ther efore the maximum error associated with the thermocouples from the multiplexer is .006MPV After the signal passes through the multiplexer it is sent to the AD board where quantiz ation error arises. The quantization error of the AD board is defined as, maxmin1 22qe NVV (B1) 116 PAGE 119 where is the voltage range given to that specif ic channel and N is the number of bits which in this case is 16. The highest voltage range assigned to any channel is 20 volts, and the resulting quantization er ror is approxim ately, Due to the high resolution of the AD board, the quantization error is negligible. Usi ng the KlineMcclintock m ethod, the total voltage error from the DAQ is shown to be, maxminVV 41.5310qeV 220.006VqeMPV. (B2) For all channels that do not requ ire gain, the error associated with the DAQ system is estimated to come only from the quantization error. Since th e quantization error is so small, it is neglected for simplicity. Pressure Drop The pressure drop is measured using a Validyne differential pressure transducer model # DP1532. For each pressure measurement a differe nt diaphragm was used, all having a different range. The manufacturer specified rela tive uncertainty of the transducer is .25%PTP of full scale. The pressure transducer is directly calib rated to a monometer using R827 oil as the fluid. The smallest increment on the monometer is 1 mm, therefore the error associated with reading the monometer, H is estimated to be 0.5 mm. Since PgH the KlineMcclintock method shows that the uncertainty in the monometer is, mHg The calibration process involves adding pressure to the monometer, recording the voltage signal and pressure drop, increa sing the monometer pressure in crementally and repeating the process until sufficient data is recorded. The resulting pressure drop ve rsus voltage calibration curves are shown for each pressure tran sducer below in Figs. B1a B1c. 117 PAGE 120 Voltage (V) 024681012 Pressure (Pa) 0 2000 4000 6000 8000 10000 12000 14000 16000 Experimental Measurements Curve Fit (y=ax) A Voltage (V) 02468 Pressure (Pa) 0 50 100 150 200 250 Experimental Measurements Curve fit (y=ax) B Figure B1. Pressure drop calibration curve used for measurement of pressure drop. A) Venturi flow meter (PT1). B) Pitot tube (PT2). C) Pressure drop across both microchannel plate and foam samples (PT3). 118 PAGE 121 Voltage (V) 024681 0 Pressure (Pa) 500 0 500 1000 1500 2000 2500 3000 3500 Experimental Measurements Curve fit (y=ax) C Figure B1 Continued. The standard error of the estimate, P est is computed from the following equation, 2Pmc estPP N (B3) where is the measured pressure and is the pressure estimated by the calibration curve. The tota l error is taken as the root sum of square s of the standard error of the estimate, the error of the monometer and the error of the pressure transducer, mPcP 222mP P estPT (B4) The sensitivity, a, of each pressu re transducer is listed below in table B1, where PT1 is the pressure transducer used for measurement of pre ssure drop across the venturi flow meter, PT2 is the pressure transducer used for measurement of pressure drop across the pitot tube, and PT3 is 119 PAGE 122 the pres sure transducer used for measurement of pressure drop across both the microchannel plate and the foam samples. Table B1. Diaphragm maximum pressure, sensitiv ity parameter, and standard error of the estimate for each pressure transducer. max(PPa )a PaV PestPa PT1 14000 1435.0 54.5 PT2 860 26.3 1.1 PT3 3500 319.5 32 Flow Meter Calibration The flow meter used here is a venturi fl ow meter manufactured by Gerand Engineering Company. The calibration procedure involves pumping the working fl uid, HFE7100, through the venturi flow meter and reco rding the pressure drop directly Once the fluid passed through the flow meter it is discharged into a graduated cylinder. Using a stopwatch, the time required to fill the graduated cylinder to the level of 800 ml is recorded along with the pressure drop simultaneously. The mass flow rate is than computed from the following equation, Vol m t (B5) Where is the fluid density, Vol is the volume of the cylinder filled up by the fluid, and t is the time is takes to fill the cylinder. The pressure drop is measured using a Validyne differential pressure transducer. The error associated with the pressure transducer is gi ven in Eq. (B2). The error associated with the stopwatch,SW is estimated to be 0.2 s. The smallest increment on the graduated cylinder is 10 ml, therefore the estimated uncertainty in the volume reading from the graduated cylinder, g c is 5 ml. Using the KlineMcclintock m ethod, the error in the mass flow rate is calculated from the following equa tion assuming zero error in the fluid density, 22 2 mcgcSWVol tt (B6) 120 PAGE 123 The m easured flow rate is than calibrated directly to the measur ed pressure drops and is shown below in Fig. (B2). The equation for flow rate versus pressure drop is given in En. (B7). b cmAP (B7) where A = and b = 0.516. 43.13810 Pressure Drop (Pa) 0200040006000800010000 Flow rate (kg/s) 0.000 0.005 0.010 0.015 0.020 0.025 0.030 0.035 0.040 Experimental Measurements curve fit (y=axb) Figure B2. Flow rate versus pr essure drop calibration curve used for measurement of flow rate within the venturi flow meter. The standard error of the estimate is computed from, 2 42.58110mmc estmm N (B8) where is the measured mass flow rate and is the mass flow rate estimated by the calib ration curve. During the actual experime nts, the pressure drop across the venturi is measured and the resulting flow rate is computed from the calibration curve. Since the error in the pressure drop device shows up in both the original calibration as well as in the actual mmcm 121 PAGE 124 experim ents, this error is account ed for twice in the overall uncerta inty of the flow meter. The total uncertainty is given as, 2 22 ,2mmmcestPm P (B9) where m P is computed from the calibration curve. After some manipulation, Eq. (B9) can be shown in the following form, 2 1 1 2 2 2 2mb b mgcSWestmm m Ab VolVol A P (B10) After plugging in the quantitative va lues in Eq. (B10), the error in flow rate can be determined as a function of flow rate only, 22 2 24.726.0625.17862.58110.3474mmm m 2. (B11) It is necessary to mention that flow rate s extracted during the actual experiments are not taken directly from the flow rate versus pressu re drop calibration found in Fig. B2. Instead a Discharge Coefficient is extracted from the calibration curve to account for the frictional pressure drop across the venturi. By correla ting the Discharge Coefficient with Reynolds number, the calibration curve can be corrected for temperature changes during experiment. Since the true flow rate versus pressure drop cu rve requires an iterative technique to solve for flow rate, the direct flow rate ve rsus pressure drop calibration from Fig. B2 was used here in the uncertainty analysis for simplification. Uncertainty in Power Measurement The heat input is determined from direct measurements of voltage and current. The voltmeter is a Keithley model 177 microvolt d mm with a manufacturers estimated relative 122 PAGE 125 uncertainty of 0.3%VMV The amp meter is a MPJA m odel 9903 having a manufacturers estimated relative uncertainty of 2%AMI Using the KlineMcclintock method, the relative uncertainty in power is taken to be, 2 2 Q VM AMVIVI (B12) Temperature The temperature sensors used are typeE ther mocouples from Omega. The manufacturers estimated uncertainty is given as 1.7TCC Each thermocouple is calibrated in a bath with a mercury thermometer. The smallest increment on the thermometer is 0.1 therefore the estimated uncertainty in the thermometer measurement, C M C is 0.05 C The thermocouples are placed in the bath within 3 cm. distance of the thermometer. The heat flux is raised in increments and the bath is allo wed to reach a steady state temper ature. At each increment, the voltage signal is taken directly from the data acquisition and the temperature is simultaneously recorded from the thermometer. Since the voltage read from the data acquisition board is read as the difference in voltage between the junction on the board and the thermocouple probe, the temperature reading is read as the difference in temperature between th e thermocouple probe and the data acquisition board. The temperature of the data acquisition board is measured at the junction using a thermister. The manufacturers estim ated relative uncertainty in the thermister, TM is When calibrating the board, the juncti on temperature reading is adju sted to match the ambient air temperature, which is measured using the mercury thermometer. Therefore the error in the junction temperature is defined as the root mean sum of squares of the error in the thermister, the error in the thermomete r, and the error in temp erature due to the error in the DAQ voltage reading, 0.025 C 123 PAGE 126 2 22 JTTHTMVT V (B13) The following calibration curve is plotted relating the temperature directly to voltage. V 1.01.21.41.61.82.02.2 TTamb (C) 60 65 70 75 80 85 90 95 Experimental Measurements Curve Fit (9th order polynomial) Figure B3. Temperature calibration curve for a typeE thermocouple. The standard error of the estimate is computed from, 2.05Tmc estTT C N (B14) where is the measured temperature difference between the bath and the junction and mT cT is that estimated by the calibration curve. Us ing the KlineMcclintock method, the total uncertainty in the thermocouples is estimated as, 2 2222TTTHestJTTCVT V (B15) 124 PAGE 127 Various Parameters Once the uncertainty in all direct meas urements, pressure, flow rate, power, and temperature, are known, the error from other various parameters can be determined. Wall Temperature Since the distance between the microchannel wall and the thermocouple probe has some uncertainty, the error in all wa ll temperature measurements, wT is computed as follows, wm cqt TT kA (B16) 22 2wTTqt acactq kAkA (B17) Vapor Quality The vapor quality is computed from an energy balance, poutin fgqmCTT x mh (B18) Since the vapor quality, x, is a function of heat input, temperatur e, and flow rate, the error in vapor quality must be accounted for accord ingly. The equation is shown here as, 22 21 2 2p xqm fgfgfgC q mhhmh 2 T (B19) Vapor Superficial Velocity Similiarly, the error is computed for the vapor superficial velocity that results from both the flow rate and vapor quality as well as the erro r in the pitot tube veloc ity that results from the pressure drop uncertainty. v cvmx u A (B20) 125 PAGE 128 22vum cvcvxm AA x (B21) Pitot Tube Velocity 12 1.2m airP u (B22) 11 1.22muP airP (B23) 126 PAGE 129 127 APPENDIX C PRESSURE GRADIENT AND HEAT TRANSFER DATA FOR FOAM SAMPLES PAGE 130 Table C1. Pressure gradient for different foams. D1 L1 L1A 10 PPI 20 PPI 40 PPI mu ms dPdz Pam mu ms dPdz Pam mu ms dPdz Pam mu ms dPdz Pam mu ms dPdz Pam mu ms dPdz Pam 0.6 282 0.5 226 0.9 338 1.1 339 1.0 316 1.0 442 0.9 382 0.7 298 1.0 435 1.5 525 1.2 457 1.2 659 1.4 689 1.1 503 1.4 705 1.9 910 1.7 804 1.6 1080 1.8 1271 1.8 1083 1.7 1103 2.4 1520 2.2 1343 2.2 1930 2.2 1921 2.4 1819 2.1 1558 2.9 2199 2.6 1985 2.6 2670 2.6 2570 2.7 2475 2.3 1973 3.5 3035 3.1 2708 3.1 3652 3.1 3324 3 3051 2.6 2309 3.9 3763 3.6 3386 3.6 4596 3.4 4040 3.4 3680 3.0 3064 4.3 4562 4.0 4149 4.0 5426 3.6 4517 3.6 4101 3.3 3644 4.6 5152 4.3 4640 4.2 6090 3.8 4947 3.7 4463 3.4 3953 4.9 5653 4.5 5077 4.4 6670 3.8 5191 3.8 4695 3.8 4849 5.0 5934 4.6 5408 4.6 7033 128 PAGE 131 Table C2. Volum etric heat tran sfer coefficient and upper fall te mperature for different foams. L1A D1 L1 mu ms vh 3WmK ,wuT C mu ms vh 3WmK,wuT C mu ms vh 3WmK,wuT C 1.6 12364 73.0 1.5 10148 79.6 1.5 10497 77.6 2.1 20458 55.6 2.1 13586 65.2 2.0 15153 61.4 2.6 32539 46.2 2.6 19374 54.2 2.5 21689 51.0 3.2 39497 42.1 3.2 24076 48.4 3.1 27058 45.7 3.7 46448 39.3 3.6 27813 45.1 3.5 31775 42.3 4.2 50438 37.8 4.1 29818 43.2 3.9 33972 40.7 4.6 54321 36.6 4.5 32247 42.0 4.4 36414 39.4 4.9 54428 36.2 4.8 33654 40.8 4.6 37285 38.6 10 PPI 20 PPI 40 PPI mu ms vh 3WmK ,wuT C mu ms vh 3WmK,wuT C mu ms vh 3WmK,wuT C 1.2 12923 59.4 1.1 15552 70.2 0.7 23016 88.7 1.8 26779 44.7 1.7 35774 47.8 1.2 42191 57.4 2.4 45288 37.8 2.2 58228 39.8 1.7 66727 45.3 3.3 70038 33.1 3.0 104675 33.9 2.3 103577 38.3 3.9 89111 31.0 3.5 155029 31.3 2.8 137817 35.0 4.5 101196 29.7 4.2 215088 29.7 3.3 171021 33.2 5.0 116760 29.1 4.8 286133 28.6 3.8 179321 31.7 5.7 122192 28.6 5.2 372559 27.9 4.1 186646 31.1 P P P P P P P P P P 129 PAGE 132 APPENDIX D NUMERIC AL CODE FOR SOLVING POROUS MEDIA TRANSPORT EQUATIONS #include PAGE 133 double cc[512]; double d[512]; double dd[512]; double Pr[512]; double ST[512]; double nu; double C; double qs; double qf; double qint[512]; double qwp; double dy; double y[512]; double Tend[512]; double ErrorT; double Tb[512]; double ErrorTb; double hr; double hl; double Ac; double As; double LL; double Tinf[512]; double L; double x[512]; double xp[512]; double yp[512]; double dxp; double dyp; double resolutionx; double dx; double H; double ErrorTT; double TTend[51200]; double Tin; double Theta0; const double pi = 3.14159265358979323846; double Nu[512]; double Errorh; double EB; double LaA; double LaB; double LaC; double LaD; double uu[512]; double yy[512]; 131 PAGE 134 double yyp[512]; double alphas[512]; double gamm as[512]; double zs[512]; double alphaff[512]; double gammaf[512]; double zf[512]; double resolutionu; double umax; double du; double Ewall[512]; double Eint[512]; double Eint1[512][512]; double Eint2[512][512]; double Eint3[512]; double Eint4[512]; double Es[512]; double Eh[512]; double Eh1[512]; double Eh2[512]; double Eh3[512]; double Tbo[512]; double Tbo1[512]; double Tbo2[512]; double Tbi[512]; double Tbi1[512]; double Tbi2[512]; double rhof = {1.16}; double Cpf = {1030}; double mu = {1.85E6}; double dyy; //10 PPI //double K = {1.98E8}; //double epsilon = {0.918}; //double kse = {9.778651}; //double asf = {1}; //double Cf = .027; //double Texp; //double h[512]; //double Um[512]; //20 PPI //double K = {1.21E8}; //double epsilon = {0.918}; //double kse = {9.778651}; //double asf = {1}; //double Cf = .021; 132 PAGE 135 //double Texp; //double h[512]; //double Um[512]; //40 PPI double K = {6.98E9}; double epsilon = {0.918}; double kse = {9.778651}; double asf = {1}; double Cf = .02; double Texp; double h[512]; double Um[512]; //L1A //double K = {1.66E8}; //double epsilon = {0.166}; //double kse = {48.62}; //double asf = {1}; //double Cf = .034; //double Texp; //double h[512]; //double Um[512]; //double dp = {.0005}; //D1 //double K = {1.66E8}; //double epsilon = {0.166}; //double kse = {97.24}; //double asf = {1}; //double Cf = .034; //double Texp; //double h[512]; //double Um[512]; //double dp = {.00065}; //L1 //double K = {1.66E8}; //double epsilon = {0.166}; //double kse = {61.88}; //double asf = {1}; //double Cf = .034; //double Texp; //double h[512]; //double Um[512]; //double dp = {.0006}; //void WriteU(void) { ofstream file; 133 PAGE 136 file.open("Velocity.m "); file << "u = ["; for (int n = 0; n <= (( int)resolutiony); n++) { file << u[1][n] << "; if (n == ((int)resolutiony)) file << "];"; else file << std::endl; } } //void WriteY(void) { ofstream file; file.open("YY.m"); file << "yp = ["; for (int n = 0; n <= (( int)resolutiony); n++) { file << y[n] << "; if (n == ((int)resolutiony)) file << "];"; else file << std::endl; } } //void WriteX(void) { ofstream file; file.open("XX.m"); file << "xp = ["; for (int n = 0; n <= (( int)resolutionx); n++) { file << x[n] << "; if (n == ((int)resolutionx)) file << "];"; else file << std::endl; } } //void WriteTs(void) { ofstream file; file.open("Temps.m"); 134 PAGE 137 file << "TTs = ["; for (int n = 0; n <= resolutionx; n++) { for (int m = 0; m <= resolutiony; m++) { file << Ts[n][m] << "; } if (n == resolutionx) file << "];"; else file << std::endl; } } //void WriteTf(void) { ofstream file; file.open("Tempf.m"); file << "TTf = ["; for (int n = 0; n <= resolutionx; n++) { for (int m = 0; m <= resolutiony; m++) { file << Tf[n][m] << "; } if (n == resolutionx) file << "];"; else file << std::endl; } } //void Writeh(void) { ofstream file; file.open("Heat.m"); file << "h40PPI = ["; for (int k = 0; k <=resolutionu; k++) { file << h[k] << "; if (k == resolutionu) file << "];"; else file << std::endl; } } 135 PAGE 138 //void WriteUm(void) { ofstream file; file.open("MeanV.m"); file << "Um40PPI = ["; for (int k = 0; k <=resolutionu; k++) { file << Um[k] << "; if (k == resolutionu) file << "];"; else file << std::endl; } } //void Writejunk(void) { ofstream file; file.open("junk.m"); file << "c = ["; for (int k = 0; k <= (int)resolutiony; k++) { file << c[k] << "; if (k == (int)resolutiony) file << "];"; else file << std::endl; } } //void Writejunk2(void) { ofstream file; file.open("junk2.m"); file << "a = ["; for (int k = 0; k <= (int)resolutiony; k++) { file << a[k] << "; if (k == (int)resolutiony) file << "];"; else file << std::endl; } } void Writejunk3(void) 136 PAGE 139 { of stream file; file.open("junk3.m"); file << "b = ["; for (int k = 0; k <= (int)resolutiony; k++) { file << b[k] << "; if (k == (int)resolutiony) file << "];"; else file << std::endl; } } void Writejunk4(void) { ofstream file; file.open("junk4.m"); file << "d = ["; for (int k = 0; k <= (int)resolutionu; k++) { file << Eh[k] << "; if (k == (int)resolutionu) file << "];"; else file << std::endl; } } //void GenerateMesh(void) { resolutiony=64; H=0.06; L=.1; dyp=(H/L)/resolutiony; dy=H/resolutiony; y[0]=0; yp[0]=0; for ( n = 1; n <= resolutiony; n++) { y[n]=y[n1]+dy; yp[n]=yp[n1]+dyp; } resolutionx=64; dx=L/resolutionx; dxp=1/resolutionx; xp[0]=0; 137 PAGE 140 x[0]=0; for ( n = 1; n <= resolutionx; n++) { x[n]=x[n1]+dx; xp[n]=xp[n1]+dxp; } dyy = H/(resolutiony); yy[0] = 0; for (n = 1; n <= floor(r esolutiony/2); n++) { yy[n] = yy[n1]+dyy; yyp[n] = yy[n]/(H/2); } resolutionu = 64; umax = 6; du = (umax1)/resolutionu; Um[0] = 1; for (n = 1; n<= resolutionu; n++) { Um[n] = Um[n1]+du; } } //void Momentum(void) { ErrorU=10; Lal=Um[k]; Lar=2*Um[k]; while (ErrorU>1E3) { La=(Lal+Lar)/2; Da[k]=K/pow(H/2,2)*1/epsilon; eta[k]=pow(epsilon,1.5)*Cf/epsilon*La*H/2*rhof/mu; LaA = 2*eta[k]/3*pow(Da[k],0.5); LaB = 1/Da[k]+4*eta[k]/3*pow(Da[k],0.5); LaC = 1/LaD*acosh(sqrt((LaA+LaB)/LaA))1; LaD = sqrt(LaA+LaB)/2; for ( n = 0; n <= floor(resolutiony/2); n++) { uu[n]=1(LaA+LaB)/LaA*pow(1/cosh(LaD*(yyp[n]+LaC)),2); } Int[k]=(1(LaA+LaB)/LaA*(tanh(La D*(1+LaC))tanh(LaD*LaC))/LaD)*La; ErrorU=fabs(Um[k]Int[k]); if ((Um[k]Int[k])<0) { Lar=La; 138 PAGE 141 } else { Lal=La; } } for (int it = 0; it <= resolutionx; it++) { for (int n = int(floor(resoluti ony/2)); n <= int(resolutiony); n++) { u[it][n] = uu[int(nfloor(res olutiony/2))]*La/Um[k]; } for (int n = 0; n <= floor(int (floor(resolutiony/2))); n++) { u[it][n] = uu[int(floor(resoluti ony/2)n)]*La/Um[k]; } } std::cout<<" Momentum Complete"< PAGE 142 } } } return 0 ; } //int TemperatureMY(void) { { for ( it = 1; it <= (int)resolutionx1; it++) { for ( n = 0; n <= (int)resolutiony; n++) { b[n]=(2/(pow(dxp,2))+2/(pow(dyp,2))+Bi[k]); bb[n]=(2/(Re[k]*Pr[k]*pow(dxp,2))+2/(Re[k]*Pr[ k]*pow(dyp,2))+Stint[k]+u[it][n]/(dxp)); } for ( n = 0; n <= (int)resolutiony1; n++) { a[n]=1/pow(dyp,2); c[n]=1/pow(dyp,2); aa[n]=(1/(Re[ k]*Pr[k]*pow(dyp,2))+v[it][n]/(2*dyp)); cc[n]=(v[it ][n]/(2*dyp)1/(Re[k]* Pr[k]*pow(dyp,2))); } a[0] = 0; aa[0] = 0; b[0] = 1; b[(int)resolutiony] = 1; a[(int)resolutiony] = 0; c[0] = 0; bb[0] = 1; bb[(int)resolutiony] = 1; aa[(int)resolutiony] = 0; cc[0] = 0; alphas[0] = b[0]; gammas[0] = c[0]/alphas[0]; for ( i = 1; i <= (int)resolutiony1; i++) { alphas[i] = b[i]a[i]*gammas[i1]; gammas[i] = c[i]/alphas[i]; } alphas[(int)resoluti ony] = b[(int)resolutiony ]a[(int)resolutiony]*gammas[(int)resolutiony1]; alphaff[0] = bb[0]; gammaf[0] = cc[0]/alphaff[0]; for ( i = 1; i <= (int)resolutiony1; i++) 140 PAGE 143 { al phaff[i] = bb[i]aa[i]*gammaf[i1]; gammaf[i] = cc[i]/alphaff[i]; } alphaff[(int)reso lutiony] = bb[(int)resolutiony]aa[(int)resolutiony]*gammaf[(int)resolutiony1]; Tend[0]=10; ErrorT=10; { for ( n = 1; n <= (int)res olutiony1; n++) { d[n]=(Bi[k]*Thetaf[it] [n]+1/pow(dxp,2)*(Thetas[it+1][n]+Thetas[it1][n])); dd[n]=Stint[k]*Thetas[it ][n]1/(Re[k]*Pr[k]*pow(dxp,2))*(Thetaf[it+1][n]+Thetaf[it1][n])u[it][n]/(dxp)*(Thetaf[it1][n]); } d[0]=1; d[(int)resolutiony]=1; dd[0]=1; dd[(int)resolutiony]=1; zs[0] = d[0]/alphas[0]; for ( i = 1; i <= (int)resolutiony; i++) { zs[i] = (d[i ]a[i]*zs[i1])/alphas[i]; } Thetas[it][(int)resolutiony] = zs[(int)resolutiony]; for ( i = (int )resolutiony1; i >= 0; i) { Thetas[it][i] = zs[i]gammas[i]*Thetas[it][i+1]; } zf[0] = dd[0]/alphaff[0]; for ( i = 1; i <= (int)resolutiony; i++) { zf[i] = (dd[i]aa[i]*zf[i1])/alphaff[i]; } Thetaf[it][(int)resolu tiony] = zf[(int)resolutiony]; for ( i = (int )resolutiony1; i >= 0; i) { Thetaf[it][i] = zf[i]gammaf[i]*Thetaf[it][i+1]; } Tend[iter]=Thetas[(i nt)resolutionx/2][(int)resolutiony]; ErrorT=fab s(Tend[iter]Tend[iter1]); } fo r (i = 0; i <= (int )resolutiony; i++) { Ts[it][i]=Thetas[ it][i]*(10025)+25; 141 PAGE 144 Tf[it][i]=Thetaf[i t][i]*(10025)+25; //cou t< PAGE 145 Pr[k]= nu/alphaf; Momentum(); TemperatureInit(); { iterTTT=0; TTend[0]=100; ErrorTT=10; EB=1E12; while (fabs(ErrorTT)>EB) { iterTTT=iterTTT+1; nu=mu/rhof; kf=.0212; kfe=epsilon*kf; C=kse/kfe; alphaf=kfe/(rhof*Cpf); Tinf[k]=25; Stint[k]=h[k]*asf*L/(rhof*Cpf*Um[k]); Bi[k]=h[k]*asf*pow(L,2)/kse; Re[k]=Um[k]*L/nu; Pr[k]=nu/alphaf; TemperatureLBY(); TemperatureMY(); TemperatureRBY(); TTend[iterTTT]=Thetas[(int)resolutionx/2][(int)re solutiony/2]+Thetas[1][1]+Thetas[(int)resoluti onx1][(int)resolutiony1]; ErrorTT=(TTend[iterTTT]TTend[iterTTT1]); } std::cout<<"for k = "< PAGE 150 Kandlikar, S .G., Steinke, M. E., 2003, "Predicting Heat Transfer During Flow Boiling In Minichannels and Microchannels," ASHRAE Transactions 109 (1), pp. 19. Kandlikar, S. G., Balasubramanian, 2004, P., An Extension of the Flow Boiling Correlation to Transition, Laminar, and Deep Laminar Fl ows in Minichannels and Microchannels, Heat Transfer Engineering 25 (3), 8693. Kandlikar, S.G., Balasubramanian, P., 2005, Experimental Study of Flow Patterns, Pressure Drop, and Flow Instabilities in Pa rallel Rectangular Minichannels, Heat Transfer Engineering 26 (3), pp. 2027. Kasza, K. E., Didascalou, T., Wambsganss, M. W., 1997, Microscale Flow Visualization of Nucleate Boiling in Small Channels: Mechanisms Influencing Heat Transfer, Shah, R. K., (Ed.), Proceedings of International Conference on Compact Heat Exchangers for Process Industries Begell house, New York, pp. 343352. Kim, S. Y., Paek, J. W., Kang, B. H., 2000, Flo w and Heat Transfer Correlations for Porous Fin in a PlateFin Heat Exchanger, Journal of Heat Transfer Vol. 122, pp. 572578. Klett, J. W., McMillan, A., Gallego, N., 2002, Carbon Foam for Electronics Cooling, National Laboratory Fuel Cell Annual Report Klimenko, V. V., 1990, A generalized correlati on for twophase forced flow heat transfer, Internat. J. Heat Mass Transfer 31, pp. 541. Kyung. I.S., Lee, S. Y., 1996, Periodic flow excu rsion in an open twophase natural circulation loop, Nuclear Engineering Design 162, 233244. Lasbet, Y., Auvity, B., Castelain, C., Peerho ssaini, 2006, H., A Chaotic HeatExchanger for PEMFC Cooling Applications, Journal of Power Sources Journal of Power Sources, Vol. 156, Issue 1, pp. 114118. Ledinegg, M., 1938, Instability of flow during natural and forced circulation, Die Warme 61, 881. Lee, H.J, Lee, S.Y., 2001, Pressure drop correl ations for twophase flow within horizontal rectangular channels with small heights, International Journal of Multiphase Flow, 27 pp. 783796. Lockhart, R. W., and Martinelli, R. C., 1949, Pr oposed correlation of data for isothermal twophase, twocomponent flow in pipes, Chemical Engineering Progress Symposium Series 45, pp. 39. Lowry, B. and Kawaji, M., 1988, Adiabatic verti cal twophase flow in narrow flow channels, AIChE Symposium Series 84 263, pp. 133. 148 PAGE 151 Lu, W ., Zhao, C., Tassou, S., 2006, Thermal analys is on metalfoam filled heat exchangers. Part I: Metalfoam filled pipes, International Journal of Heat and Mass Transfer, Volume 49, Issue 1516, pp. 27512761 Liu, Z., and Winterton, R.H.S., 1991, A genera l correlation for satura ted and subcooled flow boiling in tube and annuli, based on a nucleate pool boiling equation, Internat. J. Heat Mass Transfer 34, pp. 2759. Lui, H. T., Kodak, G., Kakac, S., 1995, Dynamical Analysis of Pressure Drop Type Oscillations with a Planar Model, International journal of Multiphase Flow Volume 21, Issue 5, Pages 851859. Mehendal, S. S., Jacobi, A. M., Shah, R. K., 2000, Fluid Flow and Heat Tr ansfer at Microand Mesoscales with Application to Heat Exchanger Design, Appl. Mech. Rev. 53 (7), 175193. Mishima, K. and Hibiki, T., 1996, Some character istics of airwa ter twophase flow in small diameter vertical tubes, Int. J. Multiphase Flow 22, pp. 703. Mishima, K., Hibiki, T. and Nishihara, H., 1993, Some characteristics of gasliquid flow in narrow rectangular ducts, Int. J. Multiphase Flow 19 1, pp. 115. MullerSteinhagen, H. Heck, K., 1986, A Simp le Friction Pressure Drop Correlation for TwoPhase Flow in Pipes, Chemical Engineering Processes Vol. 20, pp. 297308. MullerSteinhagen, H., Jamialahmadi, M. 1996, Subcooled Flow Boiling Heat Transfer to Mixtures and Solutions, Convective Flow Boiling Editor John C. Chen, Taylor & Francis, Washington D.C. Palm, B. and Tengblad, N, 1996, Cooling of el ectronics by heat pipes and thermosyphonsA review of methods and possibilities, Proc. Nat. Heat Transfer Conf. vol. 7, HTDVol. 329, pp. 97. Peterson, G. P., Duncan, A. B., and Weichold, H. M., 1993, Experimental Investigation of Micro Heat Pipes Fabricat ed in Silicon Wafers, Journal of Heat Transfer vol. 115, pp. 751756. Ramaswamy, C., Joshi, Y., Nakayama, W., John son, W. I., 1998, Com bined effects of subcooling and operating pressure on the perf ormance of a twochamber thermosyphon, Thermal and Thermomechanical Phenomena in Electronic Systems ITH ERM '98. The Sixth Intersociety Conference on, Seattle, WA, May 1998. Ribatski, G., Wojtan, L., Thome, J. R., 2006, An analysis of experiment al data and prediction methods for twophase frictional pressure dr op and flow boiling heat transfer in microscale channels, Experimental Thermal and Fluid Science Vol. 31 (1), pp. 119. Ross, H. D., and Radermacher, R., 1987, Suppression of nucleate boi ling of pure and mixed refrigerants in turbulent annular flow, Internat. J. Multiphase Flow 13, pp. 759. 149 PAGE 152 Saha, P, Zuber, N, 1974, Point of Net Vapor Generation and Vapor Void Fraction in Subcooled Boiling, Pr oceedings of 5th International Heat Transfer Conference Vol. 4, pp 175179. Schmidt, ThE., 1949, Heat Transfer Calculations for Extended Surfaces, Refrig. Eng April, pp. 351357. Steiner, D., and Taborek, J., 1992, Flow boiling heat transfer in vertical tubes correlated by an asymptotic model Heat Transfer Eng 13, pp. 43. Stenning, A. H., Veziroglu, T. N., 1967, Osc illations in TwoComponent TwoPhase Flow, NASA CR72121, Volume 1. Tadrist, L., 2006, Review on twophase fl ow instabilities in narrow spaces, International Journal of Heat and Fluid Flow Volume 28, pp. 5462. Thome, J. R., 2004, Boiling in Microchannels : A Review of Experiment and Theory, International Journal of Heat and Mass Transfer Vol. 25, pp. 128139. Thome, J. R., Dupont, V., Jacobi, A. M., 2004, Heat Transfer Model for Evaporation in Microchannels Part I: Presentation of the Model. Tran, T. N., Wambsganss, M. W., and France, D. M., 1996, Small circularand rectangularchannel boiling with two refrigerants, Internat. J. Multiphase Flow, 22 3, pp. 485. Tuckerman. D and Pease, R. E, 1981, High Performance Heal Sinking for VLS, I. IEEE Electronic Device Letters EDL2 pp.126129. Tzeng, S. C., 2007, Spatial Thermal Regulation of Aluminum Foam Heat Sink Using Sintered Porous Conductive Pipe, International Journal of Heat and Mass Transfer Volume 50, pp. 117126 Vafai, K., Kim, S.J, 1989, "Forced convection in a channel filled with a porous medium: an exact solution", ASME Journal of Heat Transfer Vol. 111 pp.11036. Van Bragt, D.D.B., 1998, Analytical Modeli ng of Boiling Water React or Dynamics, Delft University Press. Vishnyakov, V. M., 2006, Proton exchange membrane fuel cells, Vacuum Volume 80, Issue 10, The World Energy Crisis: Some Vacuumbased Solutions, 3 August 2006, Pages 10531065. Wambsganss, M.W., Jendrzejczyk, J.A. and France, D.M., 1991, Twophase flow patterns and transitions in a small, horiz ontal, rectangular channel, Int. J. Multiphase Flow 5, pp. 40 56. Wambsganss, M.W., Jendrzejczyk, J.A., France, D. M. and Obot, N.T., 1992, Frictional pressure gradients in twophase flow in a sm all horizontal rectangular channel, Exp. Thermal Fluid Sci. 5, pp. 40. 150 PAGE 153 Wang, C. C. Lee, C. J., Chang, C. T., Lin, S. P., 1999, Heat Transfer and Friction Correlation for Compact Louvered Fin and Tube Heat Exchangers, International Journal of Heat and Mass Transfer Vol. 42, pp. 19451956. Wetton, B., Promislow, K., Caglar, 2004, A., A simple thermal model of PEM fuel cell stacks, Proceedings of the Second International C onference on Fuel Cell Science, Engineering and Technology, Rochester NY, June. Yadiaroglu, G., 1981, Twophase flow in stabilities and propagation phenomena, Thermohydraulics of Twophase Systems for I ndustrial Design and Nuclear Engineering (Edited by J. M. Delhaye, M. Giot and M. L. Riethmuller), Chap. 17 McGrawHill, New York. Yadigaroglu, G., Bergles, A. E., 1969, An Expe rimental and Theoretical Study of DensityWave Oscillations in TwoPhase Flow MIT rept. DSR 746293 (HTL 7462967). Yang, X.T., Jiang, S.Y., Zhang, Y.J., 2005, M echanism Analysis on Flow Excursion of a Natural Circulation with Low Steam Quality, Nuclear Engineering and Design 235, pp. 23912406. You L., Liu H., 2002, A Twophase Flow and Tr ansport Model for the Cathode of PEM Fuel Cells, International Journal of Heat and Mass Transfer 45 (11), pp. 22772287. Zhang, Y., Ouyang, M., Qingchun, L.,Luo, J., Li X., 2004, A Model Predicting Performance of Proton Exchange Membrane Fuel Cell Stack Thermal Systems, Applied Thermal Engineering Vol. 24, pp. 501513. Zhao, C., Lu, W., Tassou, S., 2006, Thermal analys is on metalfoam filled heat exchangers. Part II: Tube heat exchangers International Journal of Heat and Mass Transfer, Volume 49, Issue 1516, Pages 27622770. Zuber, N., and Findlay, J. A., 1965, Average Volumetric Concentration in TwoPhase Flow Systems, Journal of Heat Transfer Vol. 87, pp. 453468. Zuo, Z. J., Fale, J. E., Gernert, N. J., 2000, Study of Drainage Disk Thermosyphon Condenser Design, Journal of Enhanced Heat Transfer Volume 7, Issue 3. 151 PAGE 154 152 BIOGRAPHICAL SKETCH Patrick Garrity is a Ph.D. student focusing on thermal science and fluid dynamics at the University of Florida, United States. He received his bachelors degree from the University of Massachusetts in 2004. Currently he is working on thermal management of PEM fuel cells with a focus on microchannel flow boiling. 