Numerical Simulations of Boiling in Dielectric Fluid Immersion Cooling Scenarios of High Power Electronics by Seth Nathaniel Fincher A thesis submitted to the Graduate Faculty of Auburn University in partial fulfillment of the requirements for the Degree of Master of Science Auburn, Alabama May 4, 2014 Keywords: Numerical Simulation, Two-phase, Dielectric Fluid, Electronics Cooling, Immersion Cooling, Computational Fluid Dynamics Copyright 2014 by Seth Nathaniel Fincher Approved by Roy W. Knight, Assistant Professor of Mechanical Engineering Sushil H. Bhavnani, Professor of Mechanical Engineering Daniel Harris, Associate Professor of Mechanical Engineering ii ABSTRACT As server technologies have progressed, the traditional methods of blade level cooling have become very costly. Up to this point, most server units have been cooled through forced air convection; a process that requires large air handling systems that consume a great deal of power. In this study, free convective single phase, pool boiling, and flow boiling immersion cooling of high power electronics is modeled numerically using the commercially available computational fluid dynamics (CFD) codes ANSYS Fluent. Single phase simulations of free convective, dielectric fluid immersion cooling in FC-72 were performed using Icepak, and validated against experimental results. Five 1.8 cm x 1.8 cm equally spaced die were mounted on a PCB in a cross layout. Two die were vertically aligned and two die were horizontally aligned with the center die and spaced 2.5 cm apart. Excellent agreement was found between the numerical results and experimental data. Using a five die model and a similar four die model, a numerical study was then performed to see the effect die size and spacing had on die temperature. After determining the effects of the geometry, the same setup was also modeled with using Novec 649 and HFE 7100 as a working fluid. Two-phase simulations of pool and flow boiling in subcooled Novec 649 were performed using the Eulerian two-fluid boiling model available in ANSYS Fluent known as the RPI Boiling Model. Modeled using ANSYS Icepak and imported into Fluent, four 2.4 cm x 2.4 cm die were mounted on a PCB measuring 10 cm x 10 cm in a square layout to simulate an experimental apparatus. Using various combinations of the built-in boiling model parameters in ANSYS Fluent, simulations were performed and the results were compared to experimental results. In two different cases, the combinations of boiling model parameters used showed good agreement between simulations and experimental results. After iii establishing the creation of a model in Icepak, importing it into Fluent, and modifying the boiling model through its boiling parameters as the appropriate methodologies to model boiling in Fluent, flow boiling was investigated. Using a model of a high performance small form factor, liquid immersion cooled server, flow boiling was simulated. Using subcooled Novec 649 as a working fluid, the effect of the inlet flow rate was evaluated. Results from these simulations were compared to experimental data and good agreement was observed. iv ACKNOWLEDGMENTS I would like to thank my advisory committee, Dr. Roy Knight, Dr. Sushil H. Bhavnani, and Dr. Daniel Harris for their guidance in my research. Special thanks to Dr. Knight for my undergraduate heat transfer class that got me excited about the thermal sciences and for helping me to develop the tools necessary to accomplish my personal and professional goals. I would like to thank Dr. Bhavnani for his support in senior design and throughout my graduate career. Thanks to all my colleagues in my research group, you brought camaraderie to our group that will be missed. I want to extend my gratitude to my family, Aaron, Emily, Aimee, Lydia, and especially my parents, Jerry and Brenda, for encouraging me to pursue a higher education, and giving me guidance along the way. Most importantly I want to thank my wife, Emily for her patience and support through all the long nights spent working on my thesis. She is the rock on which I am found. v TABLE OF CONTENTS Abstract ................................................................................................................................................ ii Acknowledgments ............................................................................................................................... iv List of Tables ..................................................................................................................................... viii List of Figures...................................................................................................................................... ix Nomenclature ..................................................................................................................................... xii CHAPTER 1: INTRODUCTION .......................................................................................................... 1 1.1 NUMERICAL MODELING OF BOILING ................................................................................................. 1 CHAPTER 2: LITERATURE REVIEW ............................................................................................... 3 2.1 BOILING MODEL OF KURUL AND PODOWSKI (RPI BOILING MODEL) ........................................... 3 2.2 NUMERICAL SIMULATIONS USING THE RPI BOILING MODEL ....................................................... 5 2.2.1 Numerical simulations of boiling in water .............................................................................................. 6 2.2.2 Numerical simulations of boiling in other fluids ..................................................................................... 9 2.3 NUMERICAL SIMULATIONS USING OTHER METHODS AND MODELS ......................................... 15 2.4 CONCLUSIONS FROM LITERATURE ................................................................................................... 16 CHAPTER 3: SINGLE PHASE SIMULATIONS IN ANSYS ICEPAK .............................................. 17 3.1 SINGLE PHASE MODEL ........................................................................................................................ 17 3.1.1 Five die model ..................................................................................................................................... 18 3.1.1.1 Five die model description ......................................................................................................... 18 3.1.1.2 Icepak solver settings for the five die model ............................................................................... 19 3.1.1.2 Convergence criteria and mesh independence ............................................................................. 19 3.1.1.4 Compact conduction models ....................................................................................................... 20 3.1.2 Four die model .................................................................................................................................... 21 vi 3.1.2.1 Four die model description ......................................................................................................... 21 3.1.2.2 Icepak solver settings for the four die model ............................................................................... 22 3.1.2.2 Convergence criteria and mesh independence ............................................................................. 23 3.2 SINGLE PHASE MODEL VALIDATION ................................................................................................ 24 3.3 EFFECT OF DIE SPACING ..................................................................................................................... 25 3.3.1 Normalized die spacing ....................................................................................................................... 25 3.4 NORMALIZATION OF RESULTS .......................................................................................................... 27 3.4.1 Nusselt ? Rayleigh relation .................................................................................................................. 27 3.4.2 Inclusion of Prandtl number ................................................................................................................. 30 3.4.3 Normalization applied to the four die model ......................................................................................... 33 3.4.4 Issue of model scaling ......................................................................................................................... 34 3.5 EFFECT OF BOARD COPPER ................................................................................................................ 35 CHAPTER 4: TWO-PHASE SIMULATIONS ................................................................................... 37 4.1 OVERVIEW OF THE TWO-PHASE BOILING MODEL ......................................................................... 37 4.1.1 RPI Boiling Model as implemented in Fluent ....................................................................................... 37 4.1.2 Boiling Model Parameters ................................................................................................................... 39 4.1.2.1 Bubble departure diameter........................................................................................................... 39 4.1.2.2 Bubble departure frequency and fractional area of influence ........................................................ 40 4.1.2.3 Nucleation site density ................................................................................................................ 40 4.2 POOL BOILING SIMULATIONS ............................................................................................................ 41 4.2.1 Flat plate model ................................................................................................................................... 41 4.2.1.1 Flat plate model description......................................................................................................... 41 4.2.1.2 Fluent solver settings for the two-phase flat plate model .............................................................. 43 4.2.1.3 Convergence criteria and mesh independence .............................................................................. 43 4.2.1.4 Boiling parameter functions used................................................................................................. 44 4.2.1.5 Comparison to Rohsenow?s correlation ....................................................................................... 44 4.2.2 Multi-chip module model..................................................................................................................... 45 4.2.2.1 Multi-chip module model description .......................................................................................... 46 vii 4.2.2.2 Fluent solver settings for the two-phase multi-chip module model................................................ 48 4.2.2.3 Convergence criteria and mesh independence .............................................................................. 48 4.2.2.4 Multi-chip model die modeling ................................................................................................... 50 4.2.2.5 Boiling model parameters used .................................................................................................... 50 4.2.2.6 Comparison to experimental results ............................................................................................. 51 4.3 FLOW BOILING SIMULATIONS ........................................................................................................... 57 4.3.1 Cartridge model ................................................................................................................................... 57 4.3.1.1. Cartridge model description ....................................................................................................... 58 4.3.1.2 Fluent solver settings for the two-phase cartridge model .............................................................. 60 4.3.1.3 Convergence criteria and mesh independence .............................................................................. 61 4.3.1.4 Boiling model parameters used .................................................................................................... 62 4.3.1.5 Comparison to experimental results ............................................................................................. 62 CHAPTER 5: CONCLUSIONS .......................................................................................................... 70 5.1 SINGLE PHASE SIMULATIONS ............................................................................................................ 70 5.2 TWO-PHASE POOL BOILING SIMULATIONS ..................................................................................... 70 5.3 TWO-PHASE FLOW BOILING SIMULATIONS .................................................................................... 71 5.4 RECOMMENDATIONS FOR FUTURE WORK ...................................................................................... 71 REFERENCES ................................................................................................................................... 73 APPENDIX A ? USING ICEPAK FOR CREATING A BOILING MODEL IN FLUENT: A GUIDE .................................................................................................................... 77 APPENDIX B ? SINGLE PHASE RESULTS ..................................................................................... 94 APPENDIX C ? TWO-PHASE RESULTS ....................................................................................... 100 APPENDIX D ? FLUID PROPERTIES ............................................................................................ 101 APPENDIX E ? INTERFACIAL MOMENTUM TRANSFER MODELS FOR BOILING ................ 102 viii LIST OF TABLES Table 3.1 Constants for use with the five die array .............................................................................. 32 Table 3.2 Constants for use with the four die array .............................................................................. 33 Table B.1 Single heated die baseline for five die model ....................................................................... 94 Table B.2 The effect of die spacing on temperature rise, five die model, L = 18 mm, FC-72 ................................................................................................................................. 94 Table B.3 The effect of die spacing on temperature rise, five die model, L = 24 mm, FC-72 ................................................................................................................................. 95 Table B.4 Temperature, Nusselt, and Rayleigh data, five die model, FC-72 ......................................... 96 Table B.5 Temperature, Nusselt, and Rayleigh data, five die model, Novec 649 .................................. 97 Table B.6 Temperature, Nusselt, and Rayleigh data, five die model, HFE 7100 ................................... 98 Table B.7 Temperature, Nusselt, and Rayleigh data, four die model, FC-72 ......................................... 98 Table B.8 Temperature, Nusselt, and Rayleigh data, four die model, Novec 649 .................................. 99 Table C.1 Pool boiling simulation results from the multi-chip module model .................................... 100 Table C.2 Flow boiling simulation results from the cartridge model .................................................. 100 Table D.1: Fluid properties used for simulations................................................................................ 101 ix LIST OF FIGURES Figure 2.1 Wall superheat along the channel length as predicted by: Experiment data, RPI Boiling Model, and the RPI model as modified by H.-T. Kim et al. (Kim Model) [8] .................................................................................................................. 8 Figure 2.2 Boiling curve from experimental data compared to boiling curve from Narumanchi et al. [9] .................................................................................................. 10 Figure 2.3 Comparison of Narumanchi model [9] to experimental results for: (a) jet velocity of 0.41 m/s and (b) jet velocity of 11.36 m/s....................................................................... 11 Figure 2.4 Effect of variable liquid properties of R-113 on (a) predicted void fraction and (b) predicted liquid temperature at the M.P. (Measurement plane, located 0.5m down the heated section of the concentric channel) [13] ............................................ 13 Figure 3.1 Five die model configuration .............................................................................................. 18 Figure 3.2 Cut plane of one quarter of the non-conformal mesh used with the five die model............... 20 Figure 3.3 General layout of the four die model ................................................................................... 22 Figure 3.4 Cut plane of one quarter of the mesh used with the four die model ...................................... 23 Figure 3.5 Thermal performance of five die board, experiment and simulation .................................... 24 Figure 3.6 Heat flux as a function of temperature rise, two die sizes, three die spacings, five die model, FC-72 ......................................................................................................... 26 Figure 3.7 The effect of die spacing on temperature rise, five die model, FC-72 .................................. 26 Figure 3.8 Nusselt dependence on Rayleigh number, five die model, FC-72 ........................................ 29 Figure 3.9 Nusselt dependence of Rayleigh number, five die model, S/L=1, three fluids ...................... 30 Figure 3.10 Non-dimensionalized thermal performance for the five die model, two die sizes, three die spacings, and three fluids .......................................................................... 32 Figure 3.11 Non-dimensionalized thermal performance of four die model, two die sizes, three die spacings, and two fluids ........................................................................................ 33 Figure 3.12 Effect of scaling system dimensions on predicted thermal performance, five die model, S/L = 1, FC-72 .......................................................................................... 34 x Figure 3.13 Effect of twice the copper in the PCB on thermal performance, five die model, 18 mm die, S/L = 1, FC-72 .................................................................................... 35 Figure 4.1 Two-phase flat plate model, flat plate (grey), symmetry plane (green) ................................ 42 Figure 4.2 Two-phase flat plate model temperature contours, 96,000 W/m2, water, 10 ?C subcool .................................................................................................................... 45 Figure 4.3 The two-phase multi-chip module model ............................................................................ 46 Figure 4.4 Updated multi-chip module model with single heated die ................................................... 47 Figure 4.5 Cut plane of the mesh used with the Multi-chip module ...................................................... 49 Figure 4.6 Comparison of experimental data to numerical simulations of the single heated die multi-chip module model, CASE A, Novec 649, 15?C subcool .......................... 52 Figure 4.7 Comparison of experimental data to numerical simulations of the single heated die multi-chip module model, CASE B, Novec 649, 15?C subcool ........................... 53 Figure 4.8 Comparison of experimental data to numerical simulations of the single heated die multi-chip module model, CASE C, Novec 649, 15?C subcool ........................... 54 Figure 4.9 Temperature contour, CASE C, 4 W/cm2, Novec 649, 15?C subcool ................................... 55 Figure 4.10 Contour of vapor volume fraction on the die surface, CASE C, 4 W/cm2, Novec 649, 15?C subcool................................................................................................... 56 Figure 4.11 Velocity vector plot near heated die, CASE C, 4 W/cm2, Novec 649, 15?C subcool .......... 57 Figure 4.12 Cartridge boiling model, four die (light grey), PCB (green), flow distributor (dark grey), CWH (pink), Novec 649, 15 ?C subcool......................................................... 59 Figure 4.13 Cartridge boiling model, four die (light grey), PCB (green), flow distributor (dark grey), no CWH, Novec 649, 15 ?C subcool .............................................................. 60 Figure 4.14 Cut plane of non-conformal mesh zones used in the cartridge boiling model without CWH ............................................................................................................................... 61 Figure 4.15 Cartridge model with CWH, 800 mL/min, Novec 649, 15 ?C subcool ............................... 63 Figure 4.16 Cartridge model with and without CWH, 800 mL/min, Novec 649, 15 ?C subcool ............................................................................................................................. 64 Figure 4.17 Streaklines for cartridge model without CWH, 800 mL/min, Novec 649, 15 ?C subcool ............................................................................................................................ 65 Figure 4.18 Streaklines for cartridge model without CWH, 800 mL/min, Novec 649, 15 ?C subcool ............................................................................................................................. 66 Figure 4.19 Cartridge model without CWH, 800 mL/min and 1200 mL/min, Novec 649, 15 ?C subcool ................................................................................................................... 67 xi Figure 4.20 Vapor volume fraction contour, cartridge model without CWH, 800 mL/min, 6 W/cm2, Novec 649, 15 ?C subcool ................................................................................ 68 Figure 4.21 Vapor volume fraction contour, cartridge model without CWH, 1200 mL/min, 6 W/cm2, Novec 649, 15 ?C subcool.......................................................... 69 Figure A.1 "Solve? dialogue box in Icepak .......................................................................................... 78 Figure A.2 Fluent launcher .................................................................................................................. 79 Figure A.3 "Create/Edit Materials" dialogue box ................................................................................. 81 Figure A.4 "Fluent Database Materials" dialogue box.......................................................................... 82 Figure A.5 Example of a BC dialogue box, show for a "Pressure Outlet" condition ............................. 83 Figure A.6 ?Create/Edit Mesh Interfaces? dialogue box....................................................................... 83 Figure A.7 "Multiphase Model" dialogue box...................................................................................... 85 Figure A.8 ?Secondary Phase? dialogue box ....................................................................................... 87 Figure A.9 ?Phase Interaction? dialogue box ....................................................................................... 88 Figure A.10 Setting mass transfer mechanisms using the "Phase Interaction? dialogue box.................. 89 Figure A.11 Example of variable monitor ............................................................................................ 91 xii NOMENCLATURE Subscripts l Liquid Phase v Vapor Phase Abbreviations BC Boundary Condition BCGSTAB Bi-Conjugate Gradient Stabilizer BWR Boiling Water Reactor BFBT BWR Full-size Fine-Mesh Bundle Test CCM Compact Conduction Model CEA French Alternative Energies and Atomic Energy Commission CFD Computational Fluid Dynamics CHF Critical Heat Flux CWH Chilled Water Header DNB Departure Nucleate Boiling FW First Wall IGBT Insulated-gate Bipolar Transistor ITER International Thermonuclear Experimental Reactor NC Non-Conformal NUPEC Nuclear Power Engineering Corporation xiii PCB Printed Circuit Board PWR Pressurized Water Reactor PSBT PWR Subchannel and Bundle Tests RNG Re-normalization Group SSE Standard State Enthalpy UDF User Defined Function VOF Volume of Fluid Definitions A Area, m2 Ab Fractional Area of Influence Cp Specific heat, J/kg-K DW Bubble Departure Diameter, m Db Bubble diameter, m f Frequency of Bubble Departure, s-1 g Gravitational Acceleration, m/s2 h Heat Transfer Coefficient, W/m2-K hc Single Phase Heat Transfer Coefficient, W/m2-K hl Liquid Enthalpy, kJ/kg hlv Latent Heat of Vaporization, kJ/kg hv Vapor Enthalpy, kJ/kg Ja Jacob Number k Thermal Conductivity, W/m-K L Characteristic Length, mm Nu Nusselt Number NW Nucleation Site Density, m-2 xiv P Pressure, Pa Pr Prandtl Number ?? Heat Flux, W/cm2 Ra Rayleigh Number Ra* Modified Rayleigh Number Re Reynolds Number S Characteristic Spacing, mm T Temperature, K t Periodic Time, s VW Volume of Bubble, m3 vf Vapor Volume Fraction ?Tsub Tsat ? Tl, Subcool Temperature Difference, K ? Thermal Diffusivity, m2/s ? Thermal Expansion Coefficient, K-1 ? Density, kg/m3 ? Surface Tension, N/m ? Contact Angle, degrees u Velocity, m/s ?? Dynamic Viscosity, kg/m-s ? Kinematic Viscosity, m2/s 1 CHAPTER 1: INTRODUCTION With the increased prevalence of data centers, energy consumption in the information and communication technology (ICT) sector is at an all-time high. From 2005 to 2010, worldwide energy consumption by data centers is estimated to have increased 56% [1]. In 2010, the global power consumption of data centers was estimated to be between 250 and 350 TWh annually [2]. With the progression of microprocessor design allowing for more transistors per chip in an ever decreasing footprint, liquid immersion cooling has seen a resurgence in popularity. Intel recently successfully completed a year long test of Green Revolution Cooling?s (GRC) CarnotJetTM liquid immersion system for server cooling [3]. GRC claims their system, which directly immerses servers in a mineral oil bath, can reduce data center cooling power up to 90-95% [4]. Recently, 3M has introduced a passive two-phase immersion cooling technique that offers a 90% decrease in cooling costs over traditional air cooling [5]. The purpose of this study is to evaluate the ability of the commercially available computational fluid dynamics (CFD) code, ANSYS Fluent 14.0, to model one and two-phase liquid immersion cooling of high-power electronics. 1.1 Numerical modeling of boiling The boiling model used in this study was the two-fluid Eulerian boiling model, known as the RPI Boiling Model, available in ANSYS Fluent. This model uses four main boiling parameters to capture the underlying physics of a boiling system: the nucleation site density, bubble departure diameter, bubble departure frequency, and fractional area of influence. These boiling parameters were originally developed for simulations with water, limiting their effectiveness when other fluids are involved. Variations of the RPI Boiling Model have been used successfully to model boiling in water for nuclear 2 applications [6, 7]. With slight modifications to the boiling parameters this model was also used to model boiling in water for a unique geometry [8]. Narumachi et al. was able to successfully model liquid jet impingement cooling of high power electronics with subcooled R-134 by modifying the boiling parameters [9]. Other researchers have had success using other fluids such as HFE-301 [10], R-113 [11, 12], and R-12 [11, 13], by using similar modifications to the boiling model. One of the goals of this study is to evaluate the applicability of existing boiling models when simulating boiling of a dielectric fluid used for electronics cooling. Using the electronics packaging software ANSYS Icepak, 3D models of electronics were created. Once these models were validated against experimental results, they were imported into ANSYS Fluent, where the boiling models were implemented. Along with the default functions for the boiling parameters in the RPI Boiling model, Fluent also has alternate boiling functions for two of these parameters: bubble departure diameter and nucleation site density. First, combinations of the built in functions for the boiling model parameters were simulated using a 3D model of a multi-chip module under subcooled pool boiling conditions. The results from these simulations were compared to available experimental data. Then, a small form factor, modular liquid immersion cooled server was modeled. This server model, which operates under flow boiling conditions, was simulated in Fluent and evaluated against experimental data using the methods developed with the multi-chip module. 3 CHAPTER 2: LITERATURE REVIEW This section reviews the Eulerian boiling model used in ANSYS Fluent and research that was performed using this model. First, the RPI Boiling Model by Kurul and Podowski will be introduced. Next, research simulations of boiling in water that use this model will be discussed, followed by a review of research simulations of boiling in other fluids. Finally, research simulations that use other methods and models will be discussed. 2.1 Boiling model of Kurul and Podowski (RPI Boiling Model) In this section, the wall boiling model proposed by Kurul and Podowski, known as the RPI Boiling Model, will be introduced as it is presented in the ANSYS Fluent 14.0 Theory Guide [14]. This model divides the total wall heat flux into three components: Convective heat flux, quenching heat flux, and evaporative heat flux. The wall heat flux is written in the following form: ??? = ??? + ??? +??? (2.1) where ??? is the convective heat flux, ??? is the quenching heat flux, and ??? is the evaporative heat flux. The surface of the heated is divided into two subsections, one being the fractional area ??, which is covered by nucleating bubbles, and the other being (1? ??), which is covered by fluid. The convective heat flux is expressed in the following form: ??? = ??(?? ? ??)(1 ???) (2.2) where ?? is the single phase heat transfer coefficient as calculated by Fluent, ?? is the wall surface temperature, and ?? is the liquid temperature. The quenching heat flux, which models the cyclic averaged 4 transient energy transfer resulting from liquid filling the region near the wall immediately after bubble detachment, is expressed in the following form: ??? = 2????? ?? (?? ? ??) (2.3) where ?? is the thermal conductivity of the liquid, t is the periodic time of bubble generation and departure, ?? = ??? ???? is the thermal diffusivity, and again ?? and ?? are the wall and liquid temperatures respectively. The evaporative heat flux is expressed in the following form: ??? = ?????????? (2.4) where ?? is the volume of bubble based on the bubble departure diameter, ?? is the active nucleate site density, ?? is the vapor density, ??? is the latent heat of evaporation, and ? is the bubble departure frequency. In order to close the above equations, the following parameters need to be defined: Fractional area of influence, frequency of bubble departure, nucleation site density, and bubble departure frequency. The fractional area of influence as proposed by Del Valle and Kenning [15] is expressed in the following form: Ab = KNw?Dw24 (2.5) where Dwis the bubble departure diameter, and K is an empirical constant usually set to 4, however it has been found that this value is not universal and may vary between 1.8 and 5. The following relation for this constant has been implemented based on [15]: ? = 4.8?(? ????? 80 ) (2.6) and ????? is the subcooled Jacob number defined as 5 ????? = ??????????? ???? (2.7) where ????? = ???? ? ??. In order to avoid numerical instabilities due to unbound empirical correlation for nucleation site density, the area of influence was restricted in Fluent in the following way: ?? = min (1,??????24 ) (2.8) The frequency of bubble departure as proposed by R. Cole [16] is expressed in the following form: ? = 1? = ?4?(?????)3? ??? (2.9) Nucleation site density is modeled in Fluent using an empirical correlation based on wall superheat that was proposed by Lemmert and Chawla [17], which has the following form: ?? = ??(?? ?????)? (2.10) where n = 1.805 and C = 185; when temperatures are in ?C or K and NW is in m-2. The bubble departure diameter is modeled using empirical correlations proposed by Kurul and Podowski, and is described by the following expression: ?? = min (0.0014,0.0006?? ??? 45.0) (2.11) where 0.0014 and 0.0006 have the units of meters, and 45.0 is a reference temperature difference in ?C or K. 2.2 Numerical simulations using the RPI Boiling Model In this section, different numerical studies in the literature that used the wall boiling model proposed by Kurul and Podowski will be reviewed. To begin with, numerical studies performed using water as a working fluid will be discussed. The original RPI Boiling Model was developed using water as a working fluid, therefore it is important to understand how it performs in such scenarios. Given the goals of this study, it is also of interest to understand how the RPI Boiling Model and variations upon it 6 perform using a dielectric fluid as the working fluid. So lastly, numerical studies in the literature that used other working fluids are discussed. 2.2.1 Numerical simulations of boiling in water In et al.[6] used the commercially available Computational Fluid Dynamics (CFD) codes ANSYS ?CFX 12.1 and ANSYS ?CFX 10 to simulate subcooled boiling flow in a Pressurized Water Reactor (PWR) and a Boiling Water Reactor (BWR), respectively. The CFD analysis of the PWR used four different channel subsets to represent the different areas of flow within the PWR. For each channel subset, the power and inlet temperature was varied and in two subsets, the pressure and inlet temperature also varied. The CFD analysis on the BWR used symmetry to model half of an 8x8 fuel rod assembly. The results of the PWR simulations were compared to experimental benchmark cases based on Nuclear Power Engineering Corporation (NUPEC) PWR Subchannel and Bundle Tests (PSBT). BWR simulations were compared to the BWR Full-size Fine-Mesh Bundle Test (BFBT) benchmark. The CFD simulations for the PWR subchannel predicted a void fraction and fluid density within 10% of experimental benchmarks at low inlet subcooling. Void fraction was over estimated as the inlet subcooling increased. The CFD simulations of the BWR predicted the cross-sectional averaged void fraction of the array within 1% to 10% of the benchmark for the low and high exit quality conditions, respectively. It was noted that the non-drag forces showed a strong influence on the radial distribution of void fraction. Wang et al. [7] used the CFD code ANSYS CFX-12.0 to simulate boiling heat transfer in the secondary side of a steam generator in a nuclear reactor core. The wall boiling model proposed by Kurul and Podowski was applied to model the fluid flow characteristics and heat transfer inside the steam generator tube. Default boiling parameter models were used in these simulations and the numerical predictions were validated against experimental results. Grid independence was found in all cases. Good 7 agreement was found between the predicted values of vapor volume fraction and wall temperatures and experimental data. Once the model was validated, simulations were performed to determine the effect of the tube support plate on the system. Temperature, vapor volume fraction, and liquid pathlines were inspected and found to compare reasonably well with predictions in the literature. It was found that inlet subcooling had a great effect on the thermo-fluid characteristics of the steam generator. An increase in liquid subcooling led to an increase in liquid and vapor velocity, significant increase in void fraction, and temperature variations in the steam generator primary and secondary side H-T. Kim et al. [8] used RPI boiling model by Kurul and Podowski as implemented in the CFD code ANSYS CFX-12.1 to model flow boiling in an inclined channel with a downward-facing heated upper surface. They compared their numerical results to data from experiments they also performed. It was found that while the vapor volume fraction along the channel was well represented by the RPI Boiling Model with default boiling parameters, the wall superheat was not. Therefore, the boiling model parameters were evaluated, and new parameters were introduced that accurately models the underlying physics. The following boiling model parameters in the RPI Boiling Model were modified by H.-T. Kim et al. to more accurately predict the wall superheat: Bubble departure diameter, bubble departure frequency, and nucleation site density. During boiling in an inclined channel with a downward-facing heated surface, the generated bubbles coalesce to create large slugs and slide along the heated surface. This phenomenon was observed in the experiments by H.T. Kim et al. and was also reported by other researchers. To account for this effect in the boiling model, the bubble departure diameter was set to a constant value of 9.5mm, a value determined through visual inspection. 8 Three models for nucleation site density were evaluated. By using high-speed videos recorded from the experiment, and counting the number of active nucleation sites, it was found that the best approximation to experimental data came from the following modification of the Lemmert and Chawla model: Original Lemmert and Chawla model: [17] ?? = [185(?? ?????)]1.805 (2.12) Modified Lemmert and Chawla model: [8] ?? = 0.03[185(?? ? ????)]1.805 (2.13) where NW is in m-2 when temperature are in ?C or K and. 0.03 is a unitless correction factor, developed to model boiling from an inclined, downward facing heated surface and should not be used for an upward facing or vertical heated surface. Finally, the bubble departure frequency was modified to account for the evaporation that occurs under a large slug bubble. H.-T. Kim et al. compared the predicted wall superheat from these new models to the default RPI Boiling model and experimental results, which can be seen in Figure 2.1. Figure 2.1: Wall superheat along the channel length as predicted by: Experiment data, RPI Boiling Model, and the RPI model as modified by H.-T. Kim et al. (Kim Model) [8] 9 After being modified to accurately represent the underlying physics, the wall boiling model did a good job at predicting the wall superheat. The default RPI Boiling model greatly over predicted the wall superheat. 2.2.2 Numerical simulations of boiling in other fluids Narumanchi et al. [9] numerically modeled turbulent jet impingement cooling which involved boiling. The Eulerian multiphase model in ANSYS Fluent was used for all simulations. The multiphase model is described, and a new model for bubble departure diameter is introduced. This model is compared to experimental and analytical results available in the literature. Once the new bubble departure diameter model is implemented into the multiphase model in Fluent, two numerical simulations were done to replicate experiments. Finally, a simulation of an insulated-gate bipolar transistor (IGBT) was performed using the models previously detailed. This simulation, which was of jet impingement cooling with boiling in R-134, was not validated experimentally. Using a balance between buoyancy forces and gravity on a departing bubble and neglecting liquid shear, Narumanchi developed a new model bubble departure diameter. This new model had the following form: ?? = 3[ ??2(? ????)? (?????(?????)? 1/2 ??? ) 4 ] 1/3 (2.14) When comparing this model to a model developed by Kolev, Narumanchi found that the effect of pressure and liquid velocity needed to be accounted for. Empirical correction functions for pressure and fluid velocity were introduced, giving the model the final form of: ?? = 3?1?2[ ??2(? ????)? (?????(?????)? 1/2 ??? ) 4 ] 1/3 (2.15) ?1 = 1+ 0.1503(? ? 10?5 ?1.01352) (2.16) 10 ?2 = ??6.3512|?????| (2.17) where ?? is in m/s, P is in Pascals, ?1 and ?2 are dimensionless. Default boiling parameter models are used for nucleation site density, frequency of bubble departure, and fractional area of influence. Next, this new model was implemented into the multiphase solver in Fluent, and simulations were performed to validate it against experimental results. In these experiments, a submerged water jet with a 3?C subcool impinges on a heated disk. A constant heat flux condition is imposed on the disk, and a boiling cure is generated. In Figure 2.2, the performance of the multiphase solver with new bubble departure diameter model implemented can be seen. Figure 2.2: Boiling curve from experimental data compared to boiling curve from Narumanchi et al. [9] The multiphase model was able to predict temperature rise within 20%. Once the model was validated for water as a working fluid, it was compared to experimental results which used R-113 as a working fluid. In these experiments, R-113 with a 18.5?C subcool is directly impinged onto a circular heated disk. The boiling cure can be seen below for two different flow rates. 11 Figure 2.3: Comparison of Narumanchi model [9] to experimental results for: (a) jet velocity of 0.41 m/s and (b) jet velocity of 11.36 m/s As can be seen in Figure 2.3, the new multiphase model did a god job predicting temperature rise in these cases. These simulations matched experimental results within 10%. A new simulation was performed that modeled an IGBT cooled by submerged jet impingement of R-134a with a 3?C subcool. The purpose of these simulations was to show that boiling jets are beneficial in cooling IGBTs. Because these simulations have not been experimentally verified, the results will not be discussed here. 12 Vyskocil et al. [12] used two commercially available CFD codes to simulate vertical up flow boiling of R12 in a heated pipe. Data sets from previous experiments performed at the French Alternative Energies and Atomic Energy Commission (CEA) in Grenoble, France with the ?DEBORA? test facility were used. In the experiments, the pipe was divided into three sections: the adiabatic entry section 1m in length, the constant heat flux section 3.5m in length, and the adiabatic exit section 0.5m in length. The void fraction, vapor velocity, the mean bubble diameter, and the interfacial area profiles were measured at the end of the heated section. The boiling model developed by Kurul and Podowski at was implemented in NEPTUNE_CFD and Fluent 6. The results from these two codes were comparable and provided reasonable agreement with experimental results. Pipe radial liquid temperature was found to be comparable with NEPTUNE_CFD and Fluent in all cases. The maximum void fraction predicted with the codes was well captured, while the void fraction profile was not. The mean bubble was slightly underestimated in NEPTUNE_CFD. If the same model of interfacial area transport is implemented in Fluent, it over predicts the bubble diameter in the pipe center. This can be corrected by modifying the value of turbulent kinetic energy the model uses. With both codes, grid independence was found. E. Chen et al. [13] used the CFD code ANSYS CFX-10 to simulate subcooled flow boiling experiments of R-113 in a concentric vertical annulus. The boiling model developed by Kurul and Podowski with default boiling model parameters was used to model boiling at the wall. While improvements have been suggested for the Kurul and Podowski model based modification of the boiling model parameters, few studies have looked at the effect of temperature dependent properties. Some studies have been performed on the effect saturation temperature variation had on subcooled flow boiling of water at low pressures, but no further analysis was performed on the influence of variable liquid properties on subcooled flow boiling. 13 Experimental results were used to validate the subcooled flow boiling model. The following liquid properties were varied in this study: Density, thermal conductivity, viscosity, and specific heat. In Figure 2.4, the effect of variable liquid properties can be seen. Figure 2.4: Effect of variable liquid properties of R-113 on (a) predicted void fraction and (b) predicted liquid temperature at the M.P. (Measurement plane, located 0.5m down the heated section of the concentric channel) [13] In most cases, with variable liquid properties, the wall boiling model does a better job at predicting experimental data than it does with constant liquid properties. However, as seen in Figure 2.4, this effect is relatively small. 14 Prabhudharwadkar et al. [11] validated the Eulerian two-phase model by Kurul and Podowski as implemented in ANSYS CFX-11 against experimental results. First, the models for bubble departure frequency and nucleation site density are reviewed, and the models appropriate for the following simulations are selected. Then, the boiling model was validated against experiments results of subcooled boiling of R-12 in a vertical pipe and the most accurate of the remaining boiling parameter models were selected. Finally, using the boiling parameter models as previously recommended, the boiling model was again validated against another set of experimental results which looked at subcooled boiling of R-113 in a vertical pipe. In these simulations, the original model that Kurul and Podowski used for bubble departure diameter was used. Functions for the bubble departure frequency were compared, and it was found that Cole?s model [16], which is the default in CFX, was reasonable for these simulations. Using the experimental data available in the literature, two built in functions for bubble nucleation site density were evaluated: The Lemmert and Chawla model [17] and the Kocamustaffaogullari and Ishii model [18]. The Lemmert and Chawla correlation, which is default in CFX, under-estimated the void fraction for all cases. While the Kocamustaffaogullari and Ishii model predicted the data well in all cases, it was found that its prediction deteriorated at low pressures. Overall the Kocamustaffaogullari and Ishii model produced better results that the Lemmert and Chawla model. It was noted that the best results came from the combination of the Tolubinski and Kostanchuk bubble departure diameter function with Cole?s model for bubble departure frequency and the Kocamustaffaogullari and Ishii model for nucleation site density. After finding the most appropriate combination of boiling parameter models, simulations were performed to determine the accuracy of the built in single phase turbulent wall functions and a new two- phase turbulent wall function. The results from these two wall functions were compared to measured values of turbulent kinetic energy and Reynolds stress. The two-phase turbulent wall functions showed 15 significant improvement over the single phase wall functions in predicting both the turbulent kinetic energy and Reynolds stress in the near wall region. B. Koncar et al. [10] used the wall boiling model in ANSYS CFX-12.1 to predict flow boiling of HFE-301 in a rectangular vertical channel. The wall boiling model by Kurul and Podowski was used with default boiling parameters in CFX. Predictions by these simulations were validated against experimental data. These simulations successfully predicted the main experimental data associated with Reynolds number variation and heat flux. However, an over-prediction of liquid velocity was observed in the boiling region accompanied by an under-prediction of turbulence. Agreement was improved in these areas through the use of a simple wall roughness model. 2.3 Numerical simulations using other models and methods Domalapally et al. [19] used two commercially available CFD codes to model flow boiling inside the first wall (FW) of the International Thermonuclear Experimental Reactor (ITER). In the first approach, the Rohsenow model for nucleate boiling is used along with switching seamlessly to a Volume of Fluid (VOF) approach in the CFD code STAR-CCM+. The second approach used the Bergles- Rohsenow model implemented in ANSYS Fluent 13 through the use of user defined functions (UDFs). To begin with, fully developed flow through a flat channel was modeled, followed by a detailed model of the FW hypervapotron geometry. A good balance between grid independence and computational cost was found, with the maximum cell size comparable to sizes published in works on the same subject. Numerical predictions are compared to experimental results and the performance of both models is compared in terms of accuracy and computational cost. For the flat plate model, the accuracy of both codes was comparable in predicting temperature rise for low to moderate heat flux values, while the Rohsenow model performed better at heat fluxes that were approaching the critical heat flux (CHF). The average error for Rohsenow model in the first case (V= 1 m/s) was ~5% of bulk liquid subcooling for low to moderate heat fluxes, and ~4% at higher heat 16 fluxes. While the average error for the Bergles-Rohsenow was ~8% far from CHF, and ~14% at heat fluxes approaching CHF. This increase in accuracy for the Rohsenow model is attributed to the inclusion of more detailed physics and a large number of free parameters that allow for tuning against experimental data. For the hypervapotron model, only the Rohsenow model was tested against experimental case, resulting in an average error of ~5 ?C or ~5% of the bulk fluid subcooling. 2.4 Conclusions from literature The two fluid, Eulerian boiling model available in ANSYS Fluent, known as the RPI Boiling Model, that was used in this study was introduced. It was found that this model was successfully implemented to simulate boiling in both water and dielectrics. In the cases that modeled dielectrics, modification to either or both the bubble departure diameter and nucleation site density was required. Some researchers modified these parameters by using other functions that are available in Fluent by default, while other implemented their own models developed on physical assumptions or experimental observations. Based on these findings, it was decided that this study would focus on evaluating the ability of ANSYS Fluent, and the RPI Boiling Model, to model two-phase fluid flow in dielectrics using any combination of built in boiling model parameters. 17 CHAPTER 3: SINGLE PHASE SIMULATIONS IN ANSYS ICEPAK In this chapter, single phase simulations performed in ANSYS Icepak which have been previously published [20] will be discussed. In these simulations, a printed circuit board with multi-chip modules is suspended vertically in a quiescent pool of dielectric liquid. Using existing experimental results, the accuracy of the numerical model was verified. Once a robust model was obtained, a parametric study was performed to determine the effect of die spacing and size on temperature rise. Simulations were performed using FC-72 as well as NOVEC 649 and HFE 7100. First, a detailed description of Icepak models used in the single phase simulations is given. Next, initial results are compared to experimental data to validate the model. The effect of die spacing on the five die model is then evaluated, followed by a study focusing on the normalization of the results for predictive purposes. Finally, the effect of Printed Circuit Board (PCB) copper composition will be discussed. 3.1 Single phase model In this section, the development of the five and four die single phase models will be discussed. First, the five die model will be described, followed by a review of the solver settings in Icepak. Then, the convergence criteria and mesh used in the five die model will be discussed. Next, the compact conduction models that were used to model the die will be described. After the five die model has been thoroughly explained the four die model will be introduced. The four die model does not differ drastically from the five die model. 18 3.1.1 Five die model In this section, all aspects of the five die multichip module model will be discussed. First the five die multichip model will be described, followed by a review of the factors affecting the convergence to an accurate solution. Finally, the compact conduction models that were used to model the flip-chip die will be discussed. 3.1.1.1 Five die model description Using ANSYS Icepak, a model was created to represent the test board used in experiments as reported by Sridhar et al. [21]. In these experiments, the thermal performance of five die mounted on a PCB suspended in a quiescent, subcooled liquid pool of FC-72 was tested. Each die was heated uniformly. Due to the fact that the heat flux in all cases was relatively low, a laminar solver was used. To reduce computational time, parallel processing utilizing 10 processing cores was used. Figure 3.1: Five Die Model Configuration 19 The PCB in these simulations measured 150 mm by 150 mm by 1.4 mm thick and was made out of the material FR-4. It contained 1.2 ounce copper coverage on 16% and 14% of the high and low sides, respectively. Each die is modeled using compact conduction models in Icepak, and measures 18 mm by 18 mm (L = 18 mm) and 0.6 mm in height. They were arranged in a cross pattern, centered on the PCB, and spaced 25 mm apart (S = 25 mm). The configuration of the board can be seen in Figure 3.1, which includes a picture of the experimental five die setup that is being modeled. 3.1.1.2 Icepak solver settings for the five die model The following solver settings were used in ANSYS Icepak. First order discretization schemes were used for both momentum and temperature, while the standard scheme was used for pressure. The following values for the under-relaxation factors were used: 1.0 for temperature, body forces, and joule heating potential, 0.7 for momentum, and 0.3 for pressure. The bi-conjugate gradient stabilizer (BCGSTAB) was used for the pressure, temperature, and joule heating potential linear solvers. Double precision was used. 3.1.1.3 Convergence criteria and mesh independence In order to ensure the model was correctly predicting flow patterns and temperature rises properly, appropriate convergence criteria were found through a parametric study. Flow terms converged to 1?10?6, while energy and joule heating converged to 1?10?7. The effect the size of the computational cabinet had on the accuracy of the model was also tested. To compromise between computational time and accuracy, a cabinet measuring 0.25 m x 0.25 m x 0.1 m was used. Mesh independence for this model was also found through a parametric study. Non-conformal (NC) meshes were used in all cases, which contained between 1.5 and 3 million nodes. Icepak?s default mesh type ?Mesher-HD? was used for all NC zones in all cases. A cut plane of one quarter of the non-conformal mesh that was used with the five die model can be seen in Figure 3.2. 20 Figure 3.2: Cut plane of one quarter of the non-conformal mesh used with the five die model 3.1.1.4 Compact conduction models To accurately model the complexity of the flip-chips without drastically increasing computational cost, compact conduction models (CCM) were used in Icepak. Each flip-chip package measured 6.0 mm x 6.0 mm and 0.6 mm tall. Nine of these packages were arranged in a square array to create each 18 mm x 18 mm x 0.6 mm tall die. Power generated was specified on a per-package basis, and combined to equal the desired per-die power generation. The individual packages remained the same in all different size die models. If a larger or smaller die was needed, die would be scaled by 6 mm, or one package. 21 In Icepak, the flip-chip package substrate material was specified as ?Si-Typical?, from the default material library, and measured 0.3 mm thick. No vias were specified in the package. Traces were specified in the package, made out of ?Cu-Pure?, and measured 0.035 mm thick. Top and bottom trace coverage percentage was set at 25, while internal layer coverage was specified to be 100 percent. Fifteen rows of solder balls in both X and Y directions were specified using a peripheral array. Thirteen rows of solder balls were suppressed in both X and Y directions without central thermal balls. Each solder ball was 0.125mm in diameter and was made of eutectic tin-lead solder, ?Solder-pb39.2_sn60.8?, from the material library. The shape of the balls was specified as ?block? in Icepak. ?Si-Typical? in the Icepak material library was used for both the die and die underfill material. The die underfill was 0.1 mm thick. 3.1.2 Four die model In this section, the single phase four die model will be discussed. This model is not drastically different from the five die model. First, the four die model will be described, and its differences and similarities to the five die model will be highlighted. Then, the solver settings, convergence criteria, and mesh settings will be reviewed. 3.1.2.1 Four die model description The only major difference that the four die model has from the five die model is its use of four equally heated die. These die are mounted in a square array centered about the center of the PCB. All four die are equally spaced, have a characteristic length, L, and a characteristic spacing, S. Compact conduction models are the same in both five and four die models, meaning, that each die is made of the same flip-chip packages as defined in earlier. The PCB used in the four die model was of the same size, thickness, and copper composition as that of the five die model. The computational cabinet was also the same size in both four and five die models. In Figure 3.2, the general layout of the four die model can be seen. 22 Figure 3.3: General layout of the four die model 3.1.2.2 Icepak solver settings for the four die model In most areas, the solver settings for the four die model were the same as the five die model. First order discretization schemes were still used for both momentum and temperature. The four die model also still used the standard scheme for pressure. The under-relaxation factors for temperature, body forces, and joule heating potential remained set at 1.0. In the four die however, the under-relaxation factors for momentum and pressure were changed, from 0.7 to 0.55 for momentum, and from 0.3 to 0.2 23 for pressure. Double precision was used in the four die model. The bi-conjugate gradient stabilizer (BCGSTAB) was used for the pressure, temperature, and joule heating potential linear solvers. 3.1.2.3 Convergence criteria and mesh independence Convergence settings for the four die model are the same as those in the five die model. Flow terms converged to 1?10?6, while energy and joule heating converged to 1?10?7. The computational cabinet measured 0.25 m x 0.25 m x 0.1 m, and NC meshes were used in all cases, which contained between 1.5 and 3 million nodes. Icepak?s default mesh type ?Mesher-HD? was used for all NC zones in all cases. A cut plane of one quarter of the mesh used with the four die model can be seen in Figure 3.4. Figure 3.4: Cut plane of one quarter of the mesh used with the four die model 24 3.2 Single phase model validation To validate the numerical model of the five die test board used in experiments by Sridhar et al. [21], the maximum temperature rise from each was compared. A comparison of these data sets can be seen below in Figure 3.3. The heat flux (??) reported here is the total heat dissipated from the die (q), divided by the surface area of that die (L?). This area does not account for the heat spreading in the board or convection from the back of the board. ?T reported here is the maximum observed temperature in the die minus the ambient temperature of the pool. Figure 3.5: Thermal performance of five die board, experiment and simulation Excellent agreement was found between numerical and experimental results. As can be seen in the image taken from numerical simulations in Figure 3.3, the temperature rise in each die was almost equal. 25 3.3 Effect of die spacing Having demonstrated that the single phase model is accurately predicting the thermal performance of the multi-chip module, studies were performed to determine the effect of die spacing on temperature rise. Using variations of the five die model, simulations were performed using two die sizes, 18 and 24 mm , and three dimensionless spacings, S/L = 1, 0.5, 0.25, for several heat fluxes. Attempts were made to discern the effect of die spacing through by plotting the data using simple geometric nondimensionalization. While some useful information comes out of this, it becomes clear that the thermal performance is more complicated than just being a function of normalized die spacing, or heat flux. 3.3.1 Normalized die spacing The data from the simulations of two different die sizes, three different dimensionless die spacings, and range of heat fluxes was plotted in Figure 3.4. In this figure, the maximum temperature rise (?T) is plotted against heat flux (??). In all cases, the relationship between ?? and ?T had the same trend, meaning, all the lines had the same slope. As die spacing decreases, the temperature of the die at a given heat flux increases. Viewing the data in this way does not lead to any other insights as the lines are crisscrossed and intertwined. To understand more about how the model performs thermally, the results were normalized and plotted in Figure 3.5. If a large enough board was used and the die were moved far enough apart, each die would act independently. Meaning, each die would act as if it was the only die mounted on the board. The vertical axis in Figure 3.5 is the temperature rise divided by the temperature rise of a single die mounted on the board. A value of one for this normalized temperature rise would indicate that the thermal performance of the multichip module is the same as a single die mounted on the same board. Any value above one indicates the fractional increase of temperature rise above this case. The horizontal axis in Figure 3.5 is the dimensionless grid spacing, S/L. 26 Figure 3.6: Heat flux as a function of temperature rise, two die sizes, three die spacings, five die model, FC-72 Figure 3.7: The effect of die spacing on temperature rise, five die model, FC-72 27 It can be seen in Figure 3.5 that as die spacing decreases, the temperature rise increases. The simulated thermal image in Figure 3.5 shows that the maximum observed temperature occurs on the center die. This is due to the interaction of the vertical plumes coming from each die. When spacing is sufficiently decreased, there is a point at which the plumes coalesce, resulting in an increase in temperature of the center die. One observation of note from Figure 3.5 is that the die spacing can be significantly closer than 25 mm without significantly increasing die temperature. On this five die board, with 18 mm die heated at 2 W/cm?, S/L can be reduced to a value of about 0.3, or 6 mm die spacing, before the resultant temperature rise exceeds 10% of that seen at 25 mm spacing. Despite this insight, it is still obvious that the thermal performance cannot be easily captured by a simple function of dimensionless spacing, or heat flux. In order to come up with a general function that would predict the thermal performance of this system, other non-dimensional groupings were considered. 3.4 Normalization of results In this section, the results from the five and four die models will be generalized using dimensional similitude. First, the thermal performance of the five die model will be investigated using the Nusselt and Rayleigh numbers. Next, in order to account for the use of different fluids, a Prandtl number effect was incorporated. It was found that using these non-dimensional groups and accounting for the Prandtl number of different fluids, that the thermal performance of a given geometry could be predicted. These correlations were then applied to the four die model. Finally, the effect the copper content of the PCB was determined. 3.4.1 Nusselt-Rayleigh relation In natural convection heat transfer, dimensional similitude will exist between appropriately scaled cases for a given geometry and flow by non-dimensionalizing results. In convection, the Nusselt number is used to define the heat transfer coefficient of a particular system. 28 Nusselt number,Nu ? hLk or Nu???? ? h?Lk (3.1) Where h is the heat transfer coefficient, or ?? is the average heat transfer coefficient, k is the thermal conductivity of the fluid and L is the characteristic length of the geometry. For natural convection in particular, the Nusselt number is usually a function of Rayleigh number and in some cases depends also on the Prandtl number [22]. 3g T LR a y le ig h n u m b e r , R a ??? ?? (3.2) Where g is gravitational acceleration, ? is the coefficient of thermal expansion of the fluid, L is a characteristic length of the geometry, ? is the kinematic viscosity of the fluid and ? is the thermal diffusivity of the fluid. In the following work, L is the size of the die as show in Figure 3.1. ?T is a temperature difference defined as being the maximum temperature in the system minus the temperature of the fluid far from the board. In some cases, Ra*, an alternate for of Ra is used for cases involving constant heat flux [23, 24]. Ra* is defined as: Flux Rayleigh Number,Ra? ? ?????4??? (3.3) In this study, it was found that Ra produced more consistent results that Ra*. The heat transfer as defined by Newton?s Law of cooling: 2qqq h A T h A T L T? ? ? ? ??? (3.4) Where q is the heat dissipated, uniformly, in a single die and ?T is the previously defined temperature rise. Combining equations (3.1) and (3.4) yields the following: 29 qNu TLk?? (3.5) Both the Nusselt number and Rayleigh Number contain the temperature rise, ?T. Using this new Nu and Ra relation, the results from Figure 3.4 are recast and shown in Figure 3.6. It can be seen that for both die sizes and the same S/L ratio, the Nusselt number falls on the same line. Figure 3.8: Nusselt dependence on Rayleigh number, five die model, FC-72 30 3.4.2 Inclusion of Prandtl number Once die spacing had been investigated, it was desired to determine the effect of changing the fluid. FC-72, Novec 649, and HFE 7100 were simulated using two die sizes, 18 and 24 mm, and one normalized die spacing, S/L=1. These results, which can be seen in Figure 3.7, suggest that the Rayleigh number dependence is the same, but there is a Prandtl number effect. FC-72 (blue) with a Prandtl number of 12.35 for this temperature range and Novec 649 (red) with a Prandtl number of 11.96, agree very well. Although, HFE 7100 (green), which has a Prandtl number of 9.63, had the same trend, there was a significant offset. Figure 3.9: Nusselt dependence of Rayleigh number, five die model, S/L=1, three fluids In some cases, the Nusselt number as a function of the Rayleigh number alone will not capture the effect of the Prandtl number adequately. One form of the Nusselt number, as described in [22], captures the effect of both the Rayleigh and Prandtl numbers by writing the function yielding the Nusselt number in the following form: 31 nmNu CRa Pr? (3.6) where, C, n and m are constants for a particular geometry and flow regime. Using the Solver in Microsoft Excel, the data from Figure 3.7 was analyzed to find the values of n, m, and C that resulted in the best agreement for the Nu and Ra relations. The value of m that best fit was found to be 0.319. This value can be taken to be 0.3 without significant negative impact on accuracy. Rearranging equation (3.6), we get: nmNu CRaPr ? (3.7) which should agree well for the dielectric fluids, geometries, and flow regimes of this model. If Nu/Prm, with m=0.3, is plotted and the result is linear, then this will be proven to be true. Nu/Pr0.3 is plotted for three spacing ratios, S/L = 1, 0.5, and 0.25, two die sizes, L = 18 and 24 mm, and three fluids, FC-72, Novec 649, and HFE 7100 in Figure 3.8. 32 Figure 3.10: Non-dimensionalized thermal performance for the five die model, two die sizes, three die spacings, and three fluids For a given S/L ratio, all of the data correlates regardless of fluid or die size. All data representing S/L=1 is plotted using different colors and markers to show this. Data representing S/L=0.5 and 0.25 are plotted using the same marker for simplicity. In Table 3.1, the constants developed for the five die model can be seen. Table 3.1. Constants for use with the five die array n 0.3Nu CRa Pr? S/L C n 1 4.6492 0.1686 0.5 3.3873 0.1853 0.25 2.3754 0.2010 33 3.4.3 Normalization applied to the four die model Using the four die model described in Section 3.1.2.1, simulations were performed for L = 18 and 24 mm, and S/L = 1, 0.5, and 0.25. It was found that these results correlated well using Nusselt and Rayleigh correlation that was used for the five die model. These results can be seen in Figure 3.9. Figure 3.11: Non-dimensionalized thermal performance of four die model, two die sizes, three die spacings, and two fluids In Table 3.2, the constants developed for the four die model can be seen. Table 3.2. Constants for use with the four die array n 0.3Nu CRa Pr? S/L C n 1 4.0517 0.1776 0.5 3.327 0.1878 0.25 2.629 0.1977 34 3.4.4 Issue of model scaling In order for dimensional similitude to hold true, when one dimension in a system is changed, the entire system must be scaled by the same amount. When the die size in either model was scaled up from 18 to 24 mm, the rest of the system was also scaled by 4/3, this included die height, board thickness, and board size. A set of simulations were performed in order to quantify the effect scaling. Using the five die model, simulations were performed using a 24 mm die, but with all other dimensions the same as the model with 18 mm die. The effect this had on the Nusselt number can be seen in Figure 3.10. Using the correlations for the scaled system, it was found that the Nusselt number, therefore the heat transfer coefficient, was less than three percent larger than the non-scaled system. Figure 3.12: Effect of scaling system dimensions on predicted thermal performance, five die model, S/L = 1, FC-72 35 3.5 Effect of board copper The original PCB in the experiments by Sridhar et al. [21] and in the simulation data that has been reported up to this point contained 1.2 ounce copper covering 16% of the high side and 14% of the low side. As the conductivity of copper is much greater that FR-4 used in the board, any change to the amount of copper could have a huge effect on the thermal spreading, and therefore the thermal performance of the system. To quantify this effect, a simulation was performed using the five die model, 18 mm die, S/L = 1, in FC-72 with a PCB that contained twice the copper as the original. This changed the in plane conductivity from 3.84 W/m-K to 7.23 W/m-K. The results of this simulation can be seen in Figure 3.11. Figure 3.13: Effect of twice the copper in the PCB on thermal performance, five die model, 18 mm die, S/L = 1, FC-72 The Nusselt number, therefore the heat transfer coefficient, increases by 30% for a given Rayleigh number. Changing the copper in the PCB had a huge effect on the thermal performance of this 36 system. Therefore, the correlations described here should only be used for systems with PCB containing similar amounts of copper. 37 CHAPTER 4: TWO-PHASE SIMULATIONS In this chapter, two-phase simulations using the RPI Boiling Model will be presented and the results will be discussed. First, the RPI Boiling Model will be reintroduced as it is presented in the Fluent Theory Guide along with the additional built in functions for the boiling model parameters. Next, a preliminary simulation of a uniformly heated vertical flat plate submerged in a quiescent pool of subcooled water will be presented and compared to a well know empirical correlation. Then, a model of a multi-chip module submerged in a quiescent liquid pool of subcooled Novec 649 will be presented and simulated using the RPI Boiling Model. Additional built in functions for boiling model parameters will be evaluated using the multi-chip module and compared to experimental results. Finally, a model of a small form factor modular liquid immersion cooled server will be presented. Boiling will be simulated using the default RPI Boiling Model and compared to experimental results. 4.1 Overview of the two-phase boiling model In this section, the two-phase Eulerian boiling model in ANSYS Fluent that was used in this study will be reintroduced. This model uses four parameters to model the physics of boiling. Built in functions in Fluent for nucleation site density, bubble departure diameter, bubble departure frequency, and fractional area of influence will be introduced. 4.1.1 The RPI Boiling Model as implemented in Fluent In ANSYS Fluent, the Eulerian multi-phase model was used to simulate pool boiling and flow boiling inside a server cartridge. The wall heat flux was modeled using the RPI Boiling Model; it will be reintroduced here as it is presented in the ANSYS Fluent 14.0 Theory guide [14]. The RPI Boiling Model 38 divides the wall heat flux into three parts: The convective heat flux, quenching heat flux, and evaporative heat flux. The wall heat flux is written in the following form: ??? = ??? + ??? +??? (4.1) where ??? is the convective heat flux, ??? is the quenching heat flux, and ??? is the evaporative heat flux. The surface of the heated wall is divided into two subsections, one being the fractional area ??, which is covered by nucleating bubbles, and its converse, (1? ??), which is covered by fluid. The convective heat flux is expressed in the following form: ??? = ??(?? ? ??)(1 ???) (4.2) where ?? is the single phase heat transfer coefficient as calculated by Fluent, ?? is the wall surface temperature, and ?? is the liquid temperature. The quenching heat flux, which models the cyclic averaged transient energy transfer resulting from liquid filling the region near the wall immediately after bubble detachment, is expressed in the following form: ??? = 2????? ?? (?? ? ??) (4.3) where ?? is the thermal conductivity of the liquid, t is the periodic time of bubble growth and departure, ?? = ??? ???? is the thermal diffusivity, and again ?? and ?? are the wall and liquid temperatures respectively. The evaporative heat flux is expressed in the following form: ??? = ?????????? (4.4) where ?? is the volume of bubble based on the bubble departure diameter, ?? is the active nucleate site density, ?? is the vapor density, ??? is the latent heat of evaporation, and ? is the bubble departure frequency. There are four main parameters used in this model to capture the underlying physics of a given system: Bubble departure diameter, bubble departure frequency, nucleation site density, and 39 fractional area of influence. Fluent allows the user to implement built in User Defined Functions or UDFs for all four of these parameters. 4.1.2 Boiling model parameters This section is an overview of the boiling model parameters that are available to a Fluent user to specify within the RPI Boiling Model. First, the bubble departure diameter functions will be introduced followed by the built in function for bubble departure frequency and fractional area influenced by the nucleating bubble. Finally, the two built in functions for nucleation site density will be presented. 4.1.2.1 Bubble departure diameter Within the RPI Boiling Model, there are three built in models that can be used to describe the bubble diameter at departure: Tolubinski-Kostanchuk, Unal, and Kocamustafaogullari-Ishii (K-I). The Tolubinski-Kostanchuk model is used by default in Fluent and has following form: ?? = min (0.0014,0.0006?? ??? 45.0) (4.5) where 0.0014 and 0.0006 are in meters, and 45.0 is in K or ?C. This is the model originally used by Kurul and Podowski in the RPI Boiling Model. The second available model for bubble departure diameter is based on the Unal relationship [25]. Attempts to implement this model were unsuccessful, therefore it will not be discussed in detail. The final model available in Fluent is the K-I model [18]; It has the following form: ?? = 0.0012(??)0.90.0208 ?? ??(? ????) (4.6) where ?? = (?? ???)/?? (4.7) 40 and ? is the contact angle in degrees, and DW is in meters. Equation (4.6) was originally developed for boiling in water, this is reflected in the constants 0.0012 and 0.0208. Attempts were made to implement these models in different boiling scenarios; these results can be found in the following sections. 4.1.2.2 Bubble departure frequency and fractional area of influence In Fluent, only one model is available by default for the frequency of bubble departure and area of influence. The frequency of bubble departure is model by a correlation proposed by R. Cole [16] and is given by equation (2.9), while the fractional area of influence proposed by Del Valle and Kenning [15] is implemented in Fluent by equation (2.8). These models were used for all two-phase simulations performed in this study. While it is also possible to implement UDFs for these parameters, there were no investigations into the effect of using alternate models. 4.1.2.3 Nucleation site density Two built in models are available to predict the density of active nucleation sites in Fluent: the Lemmert-Chawla model [17], and the K-I model [18]. The Lemmert-Chawla model is used by default in Fluent and has the following form: ?? = ??(?? ?????)? (4.8) where n = 1.805 and C = 185 when temperatures are in ?C or K and NW is in m-2. The second model available for nucleation site density is the K-I model, which has the following form: ?? = ?? ? ??2? (4.9) where, Nw? = f(??)rc? ?4.4 (4.10) rc? = 2rc/Dw (4.11) 41 rc = 2?Tsat? vhlv?Tw (4.12) f(??) = 2.157x10?7???3.2(1 +0.0049??)4.13 (4.13) Dw is the bubble departure diameter, and ?? is described by equation (4.7). 4.2 Pool boiling simulations In this section, the simulations performed in a pool boiling environment will be discussed. First, the flat plate model will be introduced and compared to well-known empirical correlation for pool boiling. Then, the multi-chip module model will be introduced, the boiling parameter functions used will be discussed, and the results will be compared to experimental data. 4.2.1 Flat plate model In this section, the first attempt to create a model in Icepak, import it into Fluent, and simulate pool boiling will be discussed. Meant primarily as a test of proper methodology, a vertically oriented flat nickel plate suspended in a subcooled, quiescent pool of water was modeled. The default RPI Boiling Model was used, and compared to a well know analytical boiling model. This model, known as the Rohsenow model, is a semi-empirical model that uses experimentally validated constants to model boiling in different scenarios. The parameters used for the Rohsenow correlation will be detailed in section 4.2.1.4. Good agreement was found between the RPI Boiling Model in Fluent and the Rohsenow correlation. 4.2.1.1 Flat plate model description In this section, the flat plate model used for preliminary boiling simulations will be described. A single vertically oriented, flat nickel plate was suspended in a quiescent pool of water, which had a subcool of 10? C. A symmetry plane was used to cut the computational domain in half. The half nickel plate measured 0.15 m tall, 0.75 m wide and 0.001 m thick. A constant internal generation rate was set 42 that resulted in a surface heat flux of 96,000 W/m2. Due to the turbulent nature of boiling, the re- normalization group (RNG) k-? method was used with standard wall functions. The Eulerian RPI boiling model was used with an implicit scheme for volume fraction parameters. To specify the interaction between phases, models need to be specified for the coefficients for lift, drag and heat transfer. The ?boiling-ishii? model was used for the drag coefficient. The ?boiling-moraga? model was used for the lift coefficient. Lastly, the ?ranz-marshal? model was used for the interphase heat transfer coefficient. One mass transfer mechanism was specified between the liquid and vapor phase to allow for the generation and condensation of vapor. Figure 4.1: Two-phase flat plate model, flat plate (grey), symetry plane (green) 43 In Figure 4.1, an isometric view of the flat plate model can be seen showing the half plate in grey and the symmetry plane in green. To reduce computational time, parallel processing utilizing 10 computer cores was used. 4.2.1.2 Fluent solver settings for the two-phase flat plate model In this section, the solver settings used for the flat plate boiling model are described. Pressure and velocity were coupled with volume fractions. A least squares cell based discretization scheme was used for the gradient, while first order upwind discretization schemes were used in all other cases. The flow Courant number was set to 10, while the explicit relaxation factors for momentum and pressure were set to 0.55 and 0.25 respectively. The under-relaxation factors for density, body forces, vaporization mass, turbulent viscosity, and energy were set to 0.4 while the under-relaxation factors for turbulent kinetic energy and turbulent dissipation rate were set to 0.3. A bi-conjugate gradient stabilizer (BCGSTAB) was used for all applicable variables. Double precision was used. 4.2.1.3 Convergence criteria and mesh independence This model consisted of approximately one million elements, and no non-conformal meshing zones were used. The default Icepak mesher ?Mesher-HD? was used. Tests were performed for this model to ensure that mesh independence was found. Computational cabinet independence was also found, the cabinet that was used measured 0.2 m tall, 0.05 m thick, and 0.1 m wide to the half plane. The simulation was run until the residuals oscillated about a certain point, and an absolute convergence criterion of 1x10-5 was met for all terms. This was done to mimic convergence criteria used in a boiling tutorial provided by ANSYS Fluent. Because this model was used to verify that Fluent could model boiling, these convergence criteria were sufficient for this case. 44 4.2.1.4 Boiling parameter functions used For the flat plate model, the RPI Boiling Model was used with default boiling model parameters. The Tolubinski-Konstanchuk model was used for bubble departure diameter, and is given by equation (4.5). The frequency of bubble departure was modeled using the built in function by R. Cole [16], equation (2.9). Fractional area of influence was modeled by the Del Valle and Kenning [15] function, equation (2.8). The Lemmert-Chawla model used for nucleation site density [17], and is given by equation (4.8). No other combinations of boiling model parameters were used for this simulation. 4.2.1.5 Comparison to Rohsenow?s correlation In this section, the results from the two-phase flat plate simulation will be compared to the temperature rise calculated from the well know Rohsenow correlation. The Rohsenow correlation is a semi-empirical correlation that uses experimentally determined constants along with liquid and vapor properties to predict boiling heat transfer. The Rohsenow correlation has the following form: ?? ????? [ ? ?(?????)] 1/2 = ( 1 ???) 1/? ?????/? ????[???????]? ?? ? 1/? (4.8) where ?" is the wall heat flux, ?? is the dynamic viscosity of the liquid, ??? is the latent heat of vaporization, ? is the surface tension, g is the gravitational acceleration, ?? and ?? are the liquid and vapor densities, respectively, ??? is the specific heat of the liquid, ?? is the wall temperature, ???? is the liquid saturation temperature, ??? is the liquid Prandtl number, ??? is a surface-liquid interaction coefficient, r and s are experimentally determined exponents, that for water equal 0.33 and 1.7, respectively. Although different values for ??? have been determined experimentally, none are available for this exact problem. ??? was taken to be 0.006, which was developed for boiling in water on a vertical nickel tube [26]. Rearranging equation (4.8) to solve for ?? ? ???? yields the following equation: ?????????? = ?? ? ???? = ???? ?? { ??? ???? [ ??(? ????) ]1/2(???)1/? 1?? ? ??/?}? (4.9) 45 For this problem, the Rohsenow correlation yields ?????????? = 3.7 ?. The two-phase flat plate model simulation in Fluent predicted a wall temperature of 379.8 K, which results in ????? = ?? ? ???? = 6.7 ?. The flat plate was about the same temperature on its entire surface, as can be seen in Figure 4.2. Figure 4.2: Two-phase flat plate model temperature contours, 96,000 W/m2, water, 10 ?C subcool 4.2.2 Multi-chip module model In this section, the pool boiling simulations using the two-phase multi-chip module will be presented. First the multi-chip module model will be presented along with the settings used in the multi- phase model. Next, the results from using alternate boiling model parameters will be introduced, followed by the comparison of these results to experimental results. 46 4.2.2.1 Multi-chip module model description In this section, the multi-chip module model that was for pool boiling simulations will be presented. The model consisted of four die mounted on a PCB, suspended vertically in a quiescent, liquid pool of Novec 649 with a 15 ?C subcooling. This model was created to simulate the test apparatus described in [27]. All simulations were performed at heat fluxes above the transition boiling regime, ~3-4 W/cm2. The PCB measured 150 mm by 150 mm by 1.4 mm thick and was made out of the material FR- 4. It contained 1.2 ounce copper coverage on 16% and 14% of the high and low sides, respectively. The modeling of the die will be discussed in a later section. The four die were equally spaced 0.025 m from each other and arranged in a square layout, which was centered about the board. The board layout can be seen in Figure 4.3, which shows the PCB (green) and four die (grey). Figure 4.3: The two-phase multi-chip module In order to see how the model predicted local heat transfer on the die surface, only a single die was heated. This was done to ensure that the heated die would not interact with each other directly, via conduction in the board, or indirectly, by the lower die inducing fluid flow across the upper die. Only the top left die was heated. It was determined that the blocks representing the unheated die had a small 47 enough profile that they did not affect the fluid moving to the heated die in any way, therefore they were removed from the model. The updated multi-chip module with a single heated die can be seen in Figure 4.4. Figure 4.4: Updated multi-chip module model with single heated die The RNG k-? method was used with standard wall functions. When the heat flux becomes sufficiently high, the standard RPI Boiling Model in conjunction with the standard turbulent wall functions became unstable. Because it was desired to simulate heat fluxes higher than the point at which this occurred, to achieve stability, the non-equilibrium boiling model was used in conjunction with the non-equilibrium turbulent wall functions for all simulations above this point. The heat flux at which the model with standard wall functions became unstable at differed for each combination of boiling model parameters tested. A heat flux was simulated just below this point for each combination of boiling parameters, and it was found that the standard and non-equilibrium models produced the same results. To specify the interaction between phases, models need to be specified for the coefficients for lift, drag and heat transfer. The ?boiling-ishii? model was used for the drag coefficient. The ?boiling-moraga? model was used for the lift coefficient. Lastly, the ?ranz-marshal? model was used for the interphase heat 48 transfer coefficient. One mass transfer mechanism was specified between the liquid and vapor phase. To reduce computational time, parallel processing utilizing 10 computer cores was used. 4.2.2.2 Fluent solver settings for the two-phase multi-chip module model In this section, the solver settings used for the multi-chip boiling model are described. Pressure and velocity were coupled with volume fractions. A least squares cell based discretization scheme was used for the gradient, first order upwind discretization schemes were used in all other cases. Under the ?solution controls? tab in Fluent, the flow courant number was set to 10, while the explicit relaxation factors for momentum and pressure were set to 0.5 and 0.2 respectively. The under-relaxation factors for density, body forces, vaporization mass, turbulent viscosity, and energy were set to 0.6 while the under- relaxation factors for turbulent kinetic energy and turbulent dissipation rate were set to 0.4. BCGSTAB was used for all applicable variables. Double precision was used. 4.2.2.3 Convergence criteria and mesh independence This model consisted of approximately 1.4 million elements, and no non-conformal meshing zones were used. The mesh was generated using the default Icepak mesher ?Mesher-HD?. A cut plane of the mesh used can be seen in Figure 4.5. 49 Figure 4.5: Cut plane of mesh used with multi-chip module model A coarser mesh consisting of ~1.2 million elements and finer mesh consisting of ~1.7 million elements were tested to ensure mesh independence. A variation of less than 1 ?C in the maximum predicted temperature was seen between the mesh used in this study and the fine mesh consisting of ~1.7 million elements. The maximum observed temperature predicted with the coarse mesh varied ~4 ?C from that which was predicted using the mesh in this study. It was determined that this model used in this study was sufficiently mesh-independent, while remaining relatively computationally inexpensive. Computational cabinet independence was also found, the cabinet that was used measured 0.15 m tall, 0.025 m thick, and 0.15 m wide. For these simulations, convergence was determined by monitoring the area-weighted average temperature of the die surface. When it was clear that the temperature of the heated die was no longer significantly changing, the simulations were stopped. 50 4.2.2.4 Multi-chip model die modeling In the single phase simulations, Compact Conduction Models (CCMs) were used to model the heated die in both five and four die models. Attempts to use the same compact conduction die models in the two-phase simulations as were used in the single phase simulation were unsuccessful. Poor model stability was observed in all cases using CCM, none resulted in a converged solution. To simplify the die model, a silicon block was specified in Icepak of the appropriate size. Constant internal heat generation was specified in the block die model, this was changed for each case to get the desired heat flux. No apparent stability issues were encountered from using the die-block model. The die block model measured 0.024 m x 0.024 m and was 0.006 m tall. ?Si-Typical? in the Icepak material library was used as the material for the silicon die block. 4.2.2.5 Boiling model parameters used One way to change how the RPI Boiling Model predict heat transfer for fluids other than water is to use alternate correlations for the parameters that govern it. In this section, the effects of changing the combination of boiling parameters used within the boiling model will be discussed. In the literature ([9], [11]) success was found using the RPI Boiling Model to different fluid-surface combinations by changing the correlations for the nucleation site density or bubble departure diameter. Using the built in correlations available in ANSYS Fluent for bubble departure diameter and nucleation site density and the single heated die multi-ship module model, three sets of simulations were performed which will be referred to as CASE A, B, and C. In all simulations in this study, the frequency of bubble departure was modeled using the default built-in function by Cole [16], equation (2.9), and the fractional area of influence was modeled by the Del Valle and Kenning [15] function, equation (2.8), which is also the default. First, a set of simulations, referred to as CASE A, were performed using the default boiling parameters in Fluent for nucleation site density, by Lemmert and Chawla [17], equation (4.8), and bubble departure diameter, by Tolubinski and 51 Kostanchuk, equation (4.5). Next, a set of simulations were run using the default functions for nucleation site density, fractional area of influence, and frequency of bubble departure in conjunction with the K-I bubble departure diameter function [18], equation (4.6), referred to as CASE B. Finally, the default functions for bubble departure diameter, fractional area of influence, and frequency of bubble departure were used in combination with the K-I nucleation site density function [18], equation (4.9), to perform a set of simulations. This set of simulations is referred to as CASE C. The result from these simulations will presented in the following sections. 4.2.2.6 Comparison to experimental results In this section, the results from the three sets of simulations described in section 4.2.2.5 will be compared to experimental data received via personal communication from B. Ramakrishnan et al. [28] Results from CASE A will be presented. The results from this set of runs can be seen in Figure 4.5. 52 Figure 4.6: Comparison of experimental data to numerical simulations of the single heated die multi-chip module model, CASE A, Novec 649, 15?C subcool The default boiling parameters did a poor job at predicting temperature rise in the die. In the numerical model, as heat flux increased, the die temperature increased rapidly until ~ 9 W/cm2, where the rate at which temperature increased started to decrease. This trend was not seen in the experiments, where the rate at which temperature increased continued to increase with heat flux. Also, within the boiling regime the temperatures that were predicted with the numerical model were nowhere close to that which was found experimentally. The second set of simulations, referred to as CASE B, which used default boiling parameters except for bubble departure diameter, for which the K-I correlation was used, showed good agreement with experimental results within the heat flux range of interest. The results from the CASE B can be seen in Figure 4.6. 0 2 4 6 8 10 12 14 0 5 10 15 20 25 30 35 40 45 He at Flu x, q" (W /cm ?) Wall Superheat, ?T = Twall - Tpool (K) Experimental Numerical, Default Boiling Parameters 53 Figure 4.7: Comparison of experimental data to numerical simulations of the single heated die multi-chip module model, CASE B, Novec 649, 15?C subcool The predicted temperature rise agreed well with experimental results at low to moderate heat fluxes using this combination of boiling parameters. At higher heat fluxes, above ~ 10 W/cm2, the predicted temperature rise deviates from that which was found experimentally. The boiling curve that was generated with these boiling parameters more closely resembled the experimentally generated curve than the default boiling parameters did. As the heat flux increased, the rate at which temperature increased slowly increased. In the CASE C, which used default boiling parameters except for nucleation site density, for which the K-I model was used, also showed good agreement with experimental data. The results from CASE C can be seen in Figure 4.7. 0 2 4 6 8 10 12 14 0 10 20 30 40 50 He at Flu x, q" (W /cm ?) Wall Superheat, ?T = Twall - Tpool (K) Experimental Numerical, K-I Bubble Departure Diameter 54 Figure 4.8: Comparison of experimental data to numerical simulations of the single heated die multi-chip module model, CASE C, Novec 649, 15?C subcool The trend of the boiling curve that was predicted in CASE C showed good agreement with the experimental curve. As can be seen in Figure 4.7, the trend generated from the numerical simulations has a positive offset of ~ 2 K from the experimental curve. Meaning, the temperature predicted at a given heat was ~ 2 K higher than that which was found experimentally. The temperature contours for CASE A, B, and C all showed that the temperature of the die was essentially uniform. A temperature contour of CASE C, heated at 4 W/cm2, showing this uniformity, can be seen in Figure 4.8. 0 2 4 6 8 10 12 14 0 10 20 30 40 50 He at Flu x, q" (W /cm ?) Wall Superheat, ?T = Twall - Tpool (K) Experimental Numerical, K-I Nucleation Site Density 55 Figure 4.9: Temperature contour, CASE C, 4 W/cm2, Novec 649, 15?C subcool The volume fraction of the vapor phase on the surface of the die was represented reasonably well in all cases. CASE A and CASE C both predicted the maximum vapor volume fraction on the surface to be between 1.2-3.5% at 4 W/cm2, while CASE B under predicted it at ~ 0.5 % for the same heat flux. In all cases, the volume fraction on the die surface was non-uniform, with the highest predicted value occurring at the hottest spot on the die. A contour plot of the maximum vapor volume fraction on the die surface for CASE C, heated at 4 W/cm2, can be seen in Figure 4.9. In these simulations, the velocity vector s show trends that were expected under pool boiling conditions. In Figure 4.10 the velocity vectors fro CASE C, heated at 4 W/cm2, are shown on a cut plane located about the center of the heated die. As is shown in this vector plot, the fluid surrounding the heated die and board slowly creeps to the die surface. The fluid is then heated by the die and accelerates away eventually exiting out of the top of the computational cabinet. The boundary layer is being well represented and the no-slip condition is being enforced on all surfaces. 56 Figure 4.10: Contour of vapor volume fraction on the die surface, CASE C, 4 W/cm2, Novec 649, 15?C subcool 57 Figure 4.11: Velocity vector plot near heated die, CASE C, 4 W/cm2, Novec 649, 15?C subcool 4.3 Flow boiling simulations In this section, simulations with flow boiling will be introduced. In these simulations, a sealed cartridge containing filled with flow loop of subcooled liquid Novec 649 housing a multi-chip module was modeled. Boiling from four heated die was simulated using default boiling parameters. Results from these simulations were compared to experimental data via personal communication with J. Gess et al. [29] Good agreement was found between the observed and predicted maximum temperature rise of the four die. 4.3.1 Cartridge model In this section, the cartridge model that was used in the flow boiling simulations will be introduced. Solver settings, convergence criteria, and mesh settings will be discussed. Next, an overview 58 of the boiling model parameters that were used in these simulations will be given. Finally, the results from these simulations will be compared to available experimental data. 4.3.1.1 Cartridge model description A modular, small form factor server cartridge as described in [30] was modeled numerically. The computational cabinet of this model was generated to the dimensions of the interior of the cartridge used in the experimental tests. A flow loop circulated Novec 649 with a 15 ?C subcooled into and out of the cartridge. The PCB was a unique shape representing that which is used in [30], it contained 1.2 ounce copper coverage on 16% and 14% of the high and low sides, respectively. The heated die were modeled using silicon blocks in Icepak. The blocks were defined as ?Si-Typical? in the Icepak material library. A unique flow distributor designed by J. Gess was modeled to direct inlet flow over the heated die. Numerical simulations were performed assist in design and optimization of the flow distributor. The final designs for the flow distributor were obtained through personal communication with J. Gess [29]. A set of simulations were performed to determine the effect the chilled water header (CWH) that was present in experiments by J. Gess et al. [30] had on the system. The cartridge boiling model with the CWH can be seen in Figure 4.10. 59 Figure 4.12: Cartridge boiling model, four die (light grey), PCB (green), flow distributor (dark grey), CWH (pink), Novec 649, 15 ?C subcool It was found that the CWH only slightly subcooled the bulk liquid inside the cartridge, as a majority of the heated liquid that came in contact with the CWH immediately exited the model through the outlet. This slight subcool only resulted in a decrease in maximum observed temperature of ~ 0.5 ?C from the temperature predicted with the cartridge model without the CWH. The cartridge boiling model without the CWH can be seen in Figure 4.11. Results from these simulations will be presented in section 4.3.1.5. 60 Figure 4.13: Cartridge boiling model, four die (light grey), PCB (green), flow distributor (dark grey), no CWH, Novec 649, 15 ?C subcool The RNG k-? model was used with non-equilibrium wall functions to model turbulence. The non-equilibrium boiling model was used to model the wall heat flux due to boiling. 4.3.1.2 Fluent solver settings for the two-phase cartridge model In this section, the solver settings used for the cartridge boiling model are described. Pressure and velocity were coupled with volume fractions. A least squares cell based discretization scheme was used for the gradient, first order upwind discretization schemes were used in all other cases. The flow courant number was set to 5, while the explicit relaxation factors for momentum and pressure were set to 0.55 and 0.25 respectively. The under-relaxation factors for density, body forces, vaporization mass, turbulent viscosity, and energy were set to 0.6 while the under-relaxation factors for turbulent kinetic energy and turbulent dissipation rate were set to 0.3. BCGSTAB was used for all applicable variables. Double precision was used. 61 4.3.1.3 Convergence criteria and mesh independence This model consisted of approximately 1 million elements. Non-conformal meshing zones were used around the flow distributor and inlet, the outlet, and each of the heated die. The mesh was generated using the default Icepak mesher ?Mesher-HD?. Once a working mesh was established, a parametric study was performed on the mesh to ensure mesh independence by varying the maximum mesh cell size by ?25% for all non-conformal zones, the refined mesh consisted of ~1.4 million nodes, while the coarse mesh consisted of ~500,000. No significant improvement was found though refining the by reducing the maximum mesh cell size. It was found that the use of fine non-conformal mesh zones around the heated surfaces allowed for a sufficient accuracy without drastically increasing computational cost. Figure 4.14: Cut plane of non-conformal mesh zones used in the cartridge boiling model without CWH In Figure 4.12, a cut plane of the mesh that was used can be seen, showing each of the refined non- conformal meshing zones. For these simulations, convergence was determined by monitoring the area- weighted average temperature of the die surface. When it was clear that the temperature of the heated die was no longer significantly changing, the simulations were stopped. 62 4.3.1.4 Boiling model parameters used For these flow boiling simulations, a converged solution was found for only one set boiling model parameters. This was done using the default boiling parameters: Lemmert-Chawla nucleation site density function [17], equation (4.8), Tolubinski-Kostanchuk bubble departure diameter function, equation (4.5), Cole [16] bubble departure frequency function, equation (2.9), and the Del Valle and Kenning function for fractional area of influence [15], equation (2.8). Good agreement was found with these boiling model parameters, therefore the convergence issues with the other sets of parameters were not deeply investigated. Results from the simulations using default boiling model parameters will be presented in the next section. 4.3.1.5 Comparison to experimental results In this section, the results from the cartridge boiling model will be compared to experimental data received via personal communication with J. Gess et al. [29] First, the effect of the CWH was evaluated. A pair of simulations was performed with the cartridge model that included the CWH. In these simulations, the four die were heated to 4, 6, 8, and 10 W/cm 2, while flow was introduced to the cartridge at 800 mL/min. These simulations showed good agreement with experimental results. Results from these simulations are compared to experimental data in Figure 4.12. 63 Figure 4.15: Cartridge model with CWH, 800 mL/min, Novec 649, 15 ?C subcool While these simulations did a good job at predicting maximum temperature rise of the die, they were computationally expensive, requiring over 48 hours to achieve convergence in each case. It was desired to reduce computational cost of the model while still providing good agreement to experimental results; this was accomplished by removing the CWH at top of the cartridge. This new model, which is shown in Figure 4.11, was simulated for the same flow rate, 800 mL/min, and a range of heat fluxes. Results from these simulations are compared to experimental data and predictions from the model with the CWH present can be seen in Figure 4.13. 0 2 4 6 8 10 12 14 16 18 20 0 10 20 30 40 50 60 70 He at Flu x, q" ( W/c m? ) Temperature Rise, ?T = Twall - Tinlet ( K ) Experimental, 360 mL/min Experimental, 660 mL/min Experimental, 960 mL/min Experimental, 1250 mL/min Numerical, with CWH, 800 mL/min 64 Figure 4.16: Cartridge model with and without CWH, 800 mL/min, Novec 649, 15 ?C subcool As can be seen in Figure 4.13, removing the CWH from the numerical model only resulted in a slight increase in maximum observed die temperature. In the cartridge model with the CWH, after the liquid shoots up over the heated die it reaches the CWH and a majority of this fluid then exits the cartridge through the outlet. A small amount of the liquid that is affected by the CWH is recirculated through the cartridge, only slightly subcooling the bulk fluid inside. Because of this, the cartridge model with the CWH predicted temperature rises ~0.5 K below that which was predicted with the cartridge model without the CWH. Streaklines from the cartridge model with CWH can be seen in Figure 4.14 showing this effect. Likewise, streaklines from the cartridge model without the CWH can be seen in Figure 4.15, again showing a majority of the fluid that reaches the top of the cartridge does not recirculate but exits through the outlet. 0 2 4 6 8 10 12 14 16 18 20 0 10 20 30 40 50 60 70 He at Flu x, q" ( W/c m? ) Temperature Rise, ?T = Twall - Tinlet ( K ) Experimental, 360 mL/min Experimental, 660 mL/min Experimental, 960 mL/min Experimental, 1250 mL/min Numerical, with CWH, 800 mL/min Numerical, without CWH, 800 mL/min 65 Figure 4.17: Streaklines for cartridge model with CWH, 800 mL/min, Novec 649, 15 ?C subcool 66 Figure4.18: Streaklines for cartridge model without CWH, 800 mL/min, Novec 649, 15 ?C subcool Removing the CWH from the model resulted in a six fold reduction in computational time, without drastically reducing the accuracy. Good agreement was found between experimental data and predictions using the cartridge model without the CWH. Next, the effect of inlet liquid flow rate was investigated. Using the cartridge model previously developed without the CWH, simulations were performed at 1200 mL/min for a range of heat fluxes. Results from these simulations can be seen in Figure 4.14. 67 Figure 4.19: Cartridge model without CWH, 800 mL/min and 1200 mL/min, Novec 649, 15 ?C subcool Increasing the flow rate at the inlet to 1200 mL/min resulted in a reduction of maximum observed die temperature of ~1.5 ?C from the temperature predicted at 800 mL/min for a heat fluxes. The trends shown in Figure 4.14 do not match those seen in the experimental data. In the experiments, at high heat fluxes, little improvement was found through increasing the flow rate. In the simulations, increasing the flow rate resulted in a distinct shift in temperature for all heat fluxes. To shed light on this discrepancy, the contours of vapor volume fraction were investigated. Contour plots of vapor volume fraction for simulations with flow of 800 mL/min and 1200 mL/min can be seen in Figure 4.15 and Figure 4.16, respectively. 0 2 4 6 8 10 12 14 16 18 20 0 10 20 30 40 50 60 70 He at Flu x, q" ( W/c m? ) Temperature Rise, ?T = Twall - Tinlet ( K ) Experimental, 360 mL/min Experimental, 660 mL/min Experimental, 960 mL/min Experimental, 1250 mL/min Numerical, without CWH, 800 mL/min Numerical, without CWH, 1200 mL/min 68 Figure 4.20: Vapor volume fraction contour, cartridge model without CWH, 800 mL/min, 6 W/cm2, Novec 649, 15 ?C subcool 69 Figure 4.21: Vapor volume fraction contour, cartridge model without CWH, 1200 mL/min, 6 W/cm2, Novec 649, 15 ?C subcool As can be seen in Figure 4.20 and Figure 4.21, the boiling model is predicting very low vapor volume fraction. The low predictions of vapor volume fraction on the die surface indicate that the system is still dominated by convective heat loss. Equation (2.1) shows that the wall heat flux is made up of 2 boiling terms and one convective term. By increasing the flow rate, the liquid velocities throughout the system increased. At a flow rate of 1200 mL/min, the maximum velocity was 0.48 m/s, while at 800 mL/min, this value was only 0.32 m/s. This increase in velocity led to a significant increase in the convective heat loss with respect to the heat lost through the mechanism of boiling. Increasing the flow rate results in a shift in the boiling curve for all heat fluxes rather than the data collating to a single line as the heat flux increases for all flow rates. 70 CHAPTER 5: CONCLUSIONS 5.1. Single phase simulations Numerical simulations of the single phase multi-chip module model described in Chapter 3 showed excellent agreement with experimental results. Using a common form of the Nusselt number, the results were normalized for the geometry and dielectric fluid used. Suggestions for the constants used in these correlations were made for the four and five die model. This common form of the Nusselt number can also be manipulated into a form, such that given a fluid, geometry, and an operating heat flux, the temperature rise in the die can be found. These relationships can be used to find trends and rough approximations of the temperature rise ? heat flux relations for a range of geometries. 5.2. Two-phase pool boiling simulations A numerical model of a dielectric liquid immersion cooled multi-chip module was created in ANSYS Icepak, imported into ANSYS Fluent, and simulated under subcooled pool boiling conditions. The Eulerian RPI Boiling Model was used in Fluent to model the wall heat flux. Multiple built-in boiling functions are available in Fluent for different boiling parameters of the RPI Boiling Model; different combinations of these functions were tested. Three sets of simulations, referred to as CASE A, B, and C, were performed with different function combinations and the results were compared to experimental data. In CASE A, using the default boiling parameters, the thermal behavior of the die was not properly modeled. The temperature rise of the die was not well predicted and the general trend of the boiling curve did not resemble that which was found experimentally. In this case, the vapor volume fraction on the die surface was reasonably represented. In CASE B, the trend of the boiling curve roughly matched 71 experimental results, while the maximum temperature rise of the heated die was well predicted for low to moderate heat fluxes, 3-9 W/cm2. The vapor volume fraction on the surface of the heated die was not well represented in this case. In CASE C, the trends of the boiling cure showed good agreement with that which was found experimentally. However, the predicted temperature rise of the heated die was offset from the experimental values by ~ 2 ?C for all simulated heat fluxes. The vapor volume fraction on the surface of the heated die was reasonably well represented in this case. 5.3 Two-phase flow boiling simulations In the two-phase flow boiling simulations, the RPI Boiling Model was implemented with default boiling parameters. Two sets of simulations were performed with this model, each at a different inlet liquid flow rate. The temperature rise predicted with these models showed good agreement with experimental results. The trend of the boiling cure matched fairly well with experimental data. In the experiments, as the heat flux increased, the temperatures eventually coalesced to a single line regardless of inlet liquid flow rate. This behavior was not well captured by the cartridge boiling model. The vapor volume fraction on the surface of the heated die that was seen in these simulations was not representative of that which was seen in the experiments. 5.4. Recommendations for future work The focus of this study has been to evaluate the ability of the built in functions for the boiling parameters used in the RPI Boiling Model. The default and alternate functions that were used to model nucleation site density and bubble departure diameter are both semi empirical functions developed for certain fluids under certain conditions. Additional functions for these parameters that have been developed for dielectric fluids in these conditions should be sought out and implemented into the Fluent solver. In the pool boiling CASE C, the boiling curve was well predicted with a slight offset in die temperature rise. Building off this model, a more accurate representation of the bubble departure 72 diameter could be implemented to try to shift the curve to match the experimental data. One possible way to accomplish this is by using measured values of bubble departure diameter and implementing them into the RPI Boiling Model. While the flow boiling simulations were fairly successful at predicting the temperature rise seen in the experiments, the vapor volume fraction was not well represented and the trend of the simulated boiling curve did not match the experimental data. The convergence issues that prevented the use of the CASE C boiling parameter set with the cartridge boiling model should be investigated, as under pool boiling conditions this parameter set predicted vapor volume fractions on the surface of the die reasonably well. 73 REFERENCES [1] Jonathan Koomey, (August 2011). Growth in data center electricity use 2005 to 2010, http://www.analyticspress.com/datacenters.html [2] Mills, M.P., (2013, August). The Cloud Begins with Coal August [PDF]. Available: http://www.tech- pundit.com/wp-content/uploads/2013/07/Cloud_Begins_With_Coal.pdf [3] Judge, Peter. "Liquid Cooling For Data Centre Servers Passes Intel Tests."TechWeekEurope UK. N.p., 5 Sept. 2012. Web. 14 Mar. 2014. [4] Green Revolution Cooling. ?Cost savings & Comparison?, March 2014, http://www.grcooling.com/data-center-cost-savings [5] 3M. (2012). Playing It Cool. (May, 2012) [PDF]. Available: http://multimedia.3m.com/mws/mediawebserver?mwsId=SSSSSuH8gc7nZxtU5xtvMx_eevUqe17zH vTSevTSeSSSSSS--&fn=Novec%20Immersion%20Cooling.pdf [6] In, W. K., C. H. Shin, and C. Y. Lee. "CFD Simulation of Subcooled Boiling Flow in Nuclear Fuel Bundle." Paper presented at 7th International Conference on Computational Fluid Dynamics (ICCFD7), Big Island, Hawaii, July 9-13, 2012 [7] Wang, Chenglong, et al. "Numerical prediction of subcooled wall boiling in the secondary side of SG tubes coupled with primary coolant." Annals of Nuclear Energy 63 (2014): 633-645. [8] Kim, Hyoung-Tak, et al. "Flow Boiling in an Inclined Channel With Downward-Facing Heated Upper Wall." Heat Transfer Engineering 35.5 (2014): 492-500. 74 [9] Narumanchi, Sreekant, et al. "Numerical simulations of nucleate boiling in impinging jets: Applications in power electronics cooling." International Journal of Heat and Mass Transfer 51.1 (2008): 1-12. [10] Kon?ar, Bo?tjan, and Marko Matkovi?. "Simulation of turbulent boiling flow in a vertical rectangular channel with one heated wall." Nuclear Engineering and Design 245 (2012): 131-139. [11] Prabhudharwadkar, D., Bertodano, M. L., Buchanan Jr., J., ?Assessment of the heat transfer model and turbulent wall functions for CFD simulations of subcooled and saturated boiling?, 7th International Conference on Multiphase Flow (ICMF 2010), Tampa, FL, May 30-June 4, 2010. [12] Vyskocil, L., and J. Macek. "Boiling flow simulation in Neptune CFD and Fluent codes." XCFD4NRS, Grenoble, France (2008). [13] Chen, Erfeng, Yanzhong Li, and Xianghua Cheng. "CFD simulation of upward subcooled boiling flow of refrigerant-113 using the two-fluid model." Applied Thermal Engineering 29.11 (2009): 2508-2517. [14] ANSYS, ANSYS Fluent. "14.0 Theory Guide." ANSYS Inc. (2011). [15] Del Valle, V. H., Kenning, D. B. R., "Subcooled flow boiling at high heat flux". International Journal of Heat and Mass Transfer. 28(10). 1907?1920. 1985. [16] Cole, R., "A Photographic Study of Pool Boiling in the Region of the Critical Heat Flux". AIChE J. 6. 533?542. 1960. [17] Lemmert, M., Chawla, L. M., ?Influence of flow velocity on surface boiling heat transfer coefficient in Heat Transfer in Boiling?. E. Hahne and U. Grigull, Eds., Academic Press and Hemisphere, New York, NY, USA. 1977. [18] Kocamustafaogullari, G., Ishii, M., "Foundation of the Interfacial Area Transport Equation and its Closure Relations". International Journal of Heat and Mass Transfer. 38. 481?493. 1995. [19] Domalapally, Phani, et al. "CFD analysis of flow boiling in the ITER first wall." Fusion Engineering and Design 87.5 (2012): 556-560. 75 [20] Knight, R. W., Fincher, S., Bhavnani, S. H., Harris, D. K., Johnson, R. W., ?A Numerical Study of Single Phase Dielectric Fluid Immersion Cooling of Multichip Modules?, IMAPS 2012, San Diego, CA, Sept. 11-13, 2012. [21] Sridhar, A., Styslinger, S., Duron, C., Bhavnani, S.H., Knight, R.W., Harris, D., Johnson, R.W., ?Cooling of High-Performance Server modules Using Direct Immersion?, Proceedings of the ASME 2012 Summer Heat Transfer Conference, Rio Grande, Puerto Rico, July 8-12, HT2012-58433, 2012. [22] Gebhart, B., Jaluria, Y, Mahajan, R.L. and Sammakia, B., ?Buoyancy-Induced Flows and Transport?, Hemisphere Publishing Corporation, New York, 1988. [23] Misale, M. and Bergles, A.E., ?The Influence of Channel Width on Natural Convection and Boiling Heat Transfer with Simulated Microelectronic Components?, Experimental Thermal and Fluid Science, vol. 14, pp. 187-193, 1997. [24] Fujii, T. and Fujii, M., ?The Dependence of Local Nusselt Number on Prandtl Number in the Case of Free Convection Along a Vertical Surface with Uniform Heat Flux?, Int. J. Heat Mass Transfer, vol. 19, pp. 121-122, 1976 [25] Unal, H. C. "Maximum Bubble Diameter, Maximum Bubble Growth Time and Bubble Growth Rate During Subcooled Nucleate Flow Boiling of Water Up To 17.7 MN?m2 ". Int. J. Heat Mass Transfer. 19. 643?649. 1976. [26] Carey, Van P. "Liquid-vapor phase-change phenomena ? 2nd edition". New York: Taylor and Francis Group, LLC, (2008). [27] Ramakrishnan, B., Bhavnani, S. H., Gess, J., Knight, R. W., Harris, D. K., and Johnson, R. W., 2014, ?Effect of System and Operational Parameters on the Performance of an Immersion-Cooled Multichip Module for High Performance Computing,? 30th Annual Thermal Measurement, Modeling, and Management Symposium, SemiTherm, San Jose, CA, March 9-13, 2014 [28] Ramakrishnan, B., (personal communication, June 11, 2013) [29] Gess, J., (personal communication, October 10, 2013) 76 [30] Gess, J., Bhavnani, S. H., Ramakrishnan, B., Johnson, R. W., Harris, D. K., Knight, R. W., Hamilton, M., and Ellis, C., 2014, ?Investigation and Characterization of a High Performance, Small Form Factor, Modular Liquid Immersion Cooled Server Model,? 30th Annual Thermal Measurement, Modeling, and Management Symposium, SemiTherm, San Jose, CA, March 9-13, 2014 [31] Fluent, ANSYS. "14.0 User's Manual." ANSYS Inc., Canonsburg, PA (2011). [32] Moraga, F. J., Bonetto, R. T., Lahey, R. T., "Lateral forces on spheres in turbulent uniform shear flow". International Journal of Multiphase Flow. 25. 1321?1372. 1999. [33] Ranz, W. E., Marshall, Jr., W. R., "Evaporation from Drops, Part I and Part II". Chem. Eng. Prog.. 48(4). 173?180. April 1952. 77 APPENDIX A ? BOILING MODEL CREATION IN ANSYS FLUENT: A GUIDE Introduction: This document is mean to be a review of the methodologies used to simulate boiling heat transfer in ANSYS Fluent using three dimensional models created in ANSYS Icepak. The modeling and meshing tools in ANSYS Icepak are used to create the model, then the model is opened in Fluent and the multiphase solver is enabled. This is possible because Icepak and Fluent use the same file type for their models; the case (.cas) file. Due to the time consuming nature of preprocessing in Fluent with the multi-phase model enabled, it is important to follow a certain procedure to reduce this time. Therefore, the order in which steps are presented here has been optimized to reduce time spent preprocessing. Icepak Model Creation: Models should be created in Icepak and simulated in single phase before switching over to Fluent and enabling the multiphase solver. This should be done to ensure that the model is well built and does not demonstrate any stability issues. Once the single phase model has been created and proven to be working as desired, the case file needs to be created. The solve dialogue box should be used to create the case file for Fluent. The box labeled ?Don?t start solver? should be checked when creating a case file to import into Fluent. Also, when e using non-conformal meshes make sure to uncheck ?Merge zones when possible?. If merged zones are used in conjunction with certain non-conformal meshes and the multiphase solver in Fluent, it is possible that improper relationships will be created. For example, if non-conformal mesh zones are being 78 used around multiple cell zones with differing constant heat generation conditions in place, the heat generation will be set to the same value. This was not observed to be a problem if the model was simulated in Fluent using the single phase solver. An example of the ?solve? dialog box in Icepak can be seen in Figure A.1. Figure A.15: "Solve? dialogue box in Icepak 79 Starting Fluent: When starting Fluent, the appropriate settings need to be set before the programs launches. From the Fluent launcher, the solver precision, serial or parallel processing (and how many cores to use), the display options, and the dimension of the model being used can be set. The Fluent Launcher is shown in Figure A.2. Figure A.2: Fluent Launcher 80 If parallel processing is being used, the number of cores (processes) that are set in the Processing Options section should be set to the same number of cores that were used in solving the problem in Icepak. Select 3D in the Dimension section, Double Precision was used for all multiphase simulations. Model Manipulation Before Enabling Multiphase Solver: The following is a summary of the steps that need to be taken in Fluent before the multiphase solver is enabled. 1) Open Case File File > Read > Case? Open the appropriate case file in Fluent that was created in Icepak. 2) Specify Fluid Material Materials > liquid_name > Create/Edit? Modify the fluid properties to the appropriate values for the primary fluid (liquid). Drop down boxes can be used to set fluid properties to piecewise, linear, or user defined functions. Along with these options, Boussinesq relationship can also be set for liquid density. The Create/Edit Materials dialogue box can be seen in Figure A.3. For the secondary phase (vapor) it is necessary to create a new material. This is done by making a copy of water-vapor from the Fluent Database into the case file and modifying it to represent the fluid needed. Click Fluent Database? to bring up the Fluent Database Materials box, select water-vapor form the list and click copy. Once the fluid is copied into the case file, change the name appropriately and enter the properties of the desired vapor. The Fluent Database Materials dialogue box can be seen in Figure A.4. Note: 81 When changes to the water-vapor material are first made, there is a prompt to overwrite the water-vapor material. Click Ok. This is only a local change and will not alter the material saved in the Fluent Database. Figure A.3: "Create/Edit Materials" dialogue box 82 Figure A.4: "Fluent Database Materials" dialogue box 3) Boundary Conditions Boundary Conditions > boundary_name The type of boundary condition (BC) that is set in Icepak cannot be changed in Fluent. However, the values governing these BCs, i.e. constant temperature, heat flux, turbulence, fluid temperature, etc., can be changed. Most of the values are set correctly by Icepak; those that are not need to be specified by the user. This can be done by selecting the BC in question and clicking Edit? An example of a BC dialogue box can be seen in Figure A.5. 83 Figure A.5: Example of a BC dialogue box, show for a "Pressure Outlet" condition Note: When modeling natural convection and using a ?pressure-outlet? boundary condition, the default turbulent intensity may be too high for some cases. Changes to these values might be needed to ensure stable convergence. 4) Mesh Interfaces Mesh Interfaces > Create Mesh Interfaces If non-conformal meshes are used, the interfaces of the non-conformal zones need to be merged before a solution can be calculated. The draw mesh function should be used to display the particular mesh zone interfaces before merging to ensure that the correct zones are merged. In general, the adjacent interfaces that need to be merged will be next to each other in the list. The ?Create/Edit Mesh Interfaces? 84 dialogue box can be seen in Figure A.6. Figure A.6: ?Create/Edit Mesh Interfaces? dialogue box Enabling Multiphase Solver: Models > Multiphase > Edit? For subcooled boiling simulations of models created in Icepak and imported into Fluent, the Eulerian Boiling model was used. The ?Multiphase Model? dialogue box can be seen in Figure A.7. 85 Figure A.7: "Multiphase Model" dialogue box The following is an overview of the Eulerian boiling models available in ANSYS Fluent, as presented in the ANSYS Fluent 14.0 Theory Guide [14]. a) RPI Boiling Model: Best suited for sub-cooled nucleate boiling simulations. The solver does not keep track of the temperature of the secondary phase (vapor), rather it is set at the saturation temperature. b) Non-Equilibrium Boiling: Best suited for departure nucleate boiling (DNB). It uses the same fundamental equations for boiling heat flux that is used in the RPI Boiling Model. However, the non-equilibrium model solves for the temperature of the secondary phase (vapor). c) Critical Heat Flux: Used to model critical heat flux (CHF). CHF is characterized by a sharp reduction of local heat transfer coefficients and the excursion of wall surface temperatures. 86 Fluent extends the RPI boiling model from the nucleate boiling regime to CHF and post dry-out conditions. Model Manipulation After Enabling Multiphase Solver: 1) Material Properties Once the multiphase solver has been turned on, the standard state enthalpy needs to be set for both the liquid and the vapor phases. For boiling simulations it is important to set the difference between standard state enthalpy (SSE) for liquid and vapor phase to the latent heat of vaporization. This was accomplished by setting the SSE of the liquid phase to 0, and the SSE of the vapor phase to the latent heat of vaporization of the fluid. For models created in Icepak, the solid material properties will already be written into the case file. However, you can make changes to the solid properties in the same manner as the vapor phase or by using the materials already in the Fluent Database. If any different materials are used, the properties of the individual cell zones must be changed in Cell Zone Conditions box. 2) Phase Interaction Once the material has been set, the primary and secondary phases need to be assigned. By default, the default fluid that was used in Icepak will be set to phase-1 in Fluent. The newly created vapor phase should be set to phase-2. In general, the Sauter-Mean bubble diameter should be used for bubble diameter in multiphase simulations. If needed a constant value or user defined function can be used for this parameter. The ?Secondary Phase? dialogue box can be seen in Figure A.8. 87 Figure A.8: ?Secondary Phase? dialogue box Once the phase material has been set for the appropriate phases, the equations governing phase interaction need to be defined. This is done through the ?Phase Interaction? dialogue box, which can be seen in Figure A.9. 88 Phase > Interaction Figure A.9: ?Phase Interaction? dialogue box Using the ?Phase Interaction? dialogue box, relationships for drag, lift, heat, mass, interfacial area, and surface tension between phases need to be set. The following relationships should be set to their corresponding parameters: ? Drag ? ?boiling-ishii? ? Lift ? ?boiling-moraga? ? Heat ? ?ranz-marshall? ? Surface Tension ? appropriate value ? Mass ? 1 Mass Transfer Mechanism ? Interfacial Area ? ?ia-symmetric? For the mass transfer mechanism, ?From Phase? should be set to phase-1 (liquid) and ?To Phase? should be set to phase-2 (vapor). An example of this can be seen in Figure A.10. 89 Figure A.10: Setting mass transfer mechanisms using the "Phase Interaction? dialogue box Click Edit under the Mechanism heading to access the correlations used for boiling model parameters such as nucleation site density, bubble departure diameter, saturation temperature, frequency of bubble departure, area of influence, and saturation temperature. To model boiling in fluids other than water, the saturation temperature needs to be set to that of the desired fluid. It might be appropriate to use alternate boiling model parameters in certain scenarios. 3) Boundary Conditions When using a ?pressure outlet? BC, attention should be paid to the backflow temperatures for both the primary and secondary phases under the Thermal tab of ?Pressure Outlet? dialogue box. By default, the secondary phase backflow temperature is set to the same temperature of the primary phase. For the RPI boiling model, the secondary phase temperature is kept at saturation. In some cases, a discrepancy in these temperatures could lead to solution stability issues. However, some simulations were successfully run without altering these temperatures. This should be handled on a case by case basis. 90 4) Solution Methods Select Full Multiphase Coupled from the drop down menu and select Coupled with Volume Fractions. To begin with, default spatial discretization schemes should be used for all simulations. 5) Solution Controls For boiling simulations, it is best to keep the Courant number set between 1 and 20. It is recommended to set it at an initial value of 10, and adjust accordingly to improve stability. In the ANSYS Fluent 14.0 User?s Guide [31] it is suggested that the explicit relaxation factors for Pressure and Momentum should both be set to 1. However, with relaxation factors for pressure and momentum set this high, multiphase simulations proved to be unstable. Stability was found when Pressure was set between 0.2 and 0.4 and Momentum was set between 0.5 and 0.8. All other relaxation factors were set between 0.25 and 0.8. In Advanced Options, the biconjugate-gradient stabilizer (BCGSTAB) should be to be turned on to improve stability. Notes: a) It was found that the best methodology for quickly stabilizing simulations is to set the relaxation factor for each flow variable to the lowest end of the suggested range. This should be revisited once a stable solution is found as increasing these parameters can decrease overall computational time. b) In the case that higher order spatial discretization schemes are used, high order term relaxation should be considered to decrease computational time while still benefitting from a higher order schemes. An under-relaxation parameter of 0 corresponds to using only first order schemes. Likewise, setting the under-relaxation parameter to 1 corresponds to using the full higher order discretization schemes. 91 6) Monitors Monitors > Surface Monitors > Create In order to ensure that a solution of a boiling simulation has converged, in most case additional criteria need to be considered. For example, in pool boiling simulations surface temperature can be used in conjunction with traditional criteria for convergence. Additional monitors can be added for many variables that might be of interest through the use of the Monitors tab. An example of a monitor that can be created in Fluent can be seen in Figure A.11. Figure A.11: Example of variable monitor 92 7) Solution Initialization Solution Initialization > Initial Values The initial vapor temperature should be set to the saturation temperature; liquid temperature needs to be set to the bulk fluid temperature that was specified in Icepak. Turbulent kinetic energy should be set to 0.001 and adjusted if the solution fails to converge. Turbulent dissipation rate should be set to 1 to begin with, and adjusted to improve the initial convergence if needed. Notes: a) It was observed that when a ?block? was specified in Icepak, imported into Fluent, and used in a multiphase simulation, that its initial condition would be set to the temperature of the secondary (vapor) phase. This could be problematic if the object was not specified to be generating heat, because it will take a considerable amount of time for the solution to truly converge. b) If a model uses multiple unheated blocks, it is beneficial to use two different case files to generate the initial condition and run the simulation. Using the Cell Zone Conditions tab, the temperature of unheated block should be set to the bulk fluid temperature. These conditions should be implemented before the multi-phase solver is turned on. This model should be simulated for ~500 iterations and the case and data files should be saved. Then using the actual model, which should be a different case file, load the data file and continue the simulation. Doing so will reduce overall computational time, as it eliminates bad initial conditions that Fluent will automatically implement. 93 8) Run Calculation In order to run a simulation, the number of iterations needs to be specified. Enter the desired number of iterations in the Number of Iterations box. Click Calculate. Once the solution has adequately converged, the solution can be stopped and saved. Notes: a) If the solution is ?paused?, Fluent will not automatically save the case and data files. This needs to be done manually. b) It is worth noting that the solution can be stopped at any time and resumed at the same point. Variables in the model can be modified if desired while the simulation is ?paused?. c) In some cases, it can be beneficial to adjust convergence criteria during the simulation to speed up or improve convergence. This can be done by ?pausing? the simulation, altering the values as needed, and resuming the simulation. 94 APPENDIX B ? SINGLE PHASE RESULTS In this appendix, results from the single phase simulations will be presented. Table B.1contains the data used for model validation in Figure 3.5. Tables B.2 and B.3 contain temperature data from Figure 3.7 for three different heat fluxes. Tables B.4-B.8 contain temperature, Nusselt, and Rayleigh data for multiple heat fluxes that was plotted in Figures 3.8-3.10. Table B.1: Single heated die baseline for five die model L = 24 mm L = 18 mm W/cm2 ?Tsingle q ?Tsingle 2 39.53 2 35.15 1 22.06 1 19.53 0.25 6.82 0.25 5.98 Table B.2: The effect of die spacing on temperature rise, five die model, L = 18 mm, FC-72 ?? = 2 W/cm2 ?? = 1 W/cm2 ?? = 0.25 W/cm2 S S/L ?T ?T/?Tsingle ?T ?T/?Tsingle ?T ?T/?Tsingle 25 1.39 35.47 1.01 19.94 1.02 6.28 1.05 24 1.33 35.55 1.01 19.98 1.02 6.29 1.05 23 1.28 35.74 1.02 22 1.22 35.83 1.02 20.15 1.03 6.35 1.06 21 1.17 35.87 1.02 20.18 1.03 6.36 1.06 20 1.11 35.96 1.02 20.24 1.04 6.38 1.07 19 1.06 36.04 1.03 20.29 1.04 6.40 1.07 18 1.00 36.12 1.03 17 0.94 36.21 1.03 20.40 1.04 6.45 1.08 16 0.89 36.32 1.03 15 0.83 36.42 1.04 20.53 1.05 6.50 1.09 95 14 0.78 36.55 1.04 20.60 1.05 6.53 1.09 13 0.72 36.76 1.05 12 0.67 36.86 1.05 20.79 1.06 6.60 1.10 11 0.61 37.03 1.05 20.90 1.07 6.64 1.11 10 0.56 37.18 1.06 21.00 1.08 6.68 1.12 9 0.50 37.39 1.06 21.13 1.08 6.73 1.13 8 0.44 37.63 1.07 21.28 1.09 6.79 1.13 7 0.39 37.93 1.08 21.47 1.10 6.86 1.15 6 0.33 38.27 1.09 21.68 1.11 6.94 1.16 5 0.28 38.67 1.10 21.93 1.12 7.04 1.18 4 0.22 39.31 1.12 22.31 1.14 7.17 1.20 3 0.17 40.22 1.14 22.83 1.17 7.35 1.23 2 0.11 41.83 1.19 23.77 1.22 7.64 1.28 1 0.06 44.77 1.27 25.42 1.30 8.16 1.36 Table B.3: The effect of die spacing on temperature rise, five die model, L = 24 mm, FC-72 ?? = 2 W/cm2 ?? = 1 W/cm2 ?? = 0.25 W/cm2 S S/L ?T ?T/?Tsingle ?T ?T/?Tsingle ?T ?T/?Tsingle 25 1.04 40.61 1.07 22.94 1.04 7.28 1.03 24 1.00 40.67 1.07 22.98 1.04 7.29 1.03 23 0.96 40.73 1.07 23.01 1.04 7.31 1.03 22 0.92 40.79 1.08 23.06 1.05 7.33 1.03 21 0.88 40.84 1.08 23.09 1.05 7.34 1.03 20 0.83 40.95 1.08 23.15 1.05 7.37 1.04 19 0.79 41.02 1.08 23.20 1.05 7.39 1.04 18 0.75 41.12 1.09 23.26 1.05 7.41 1.04 17 0.71 41.20 1.09 23.31 1.06 7.43 1.04 16 0.67 41.34 1.09 23.40 1.06 7.46 1.05 15 0.63 41.45 1.10 23.46 1.06 7.49 1.05 14 0.58 41.58 1.10 23.54 1.07 7.52 1.05 13 0.54 41.70 1.11 23.63 1.07 7.55 1.06 12 0.50 41.83 1.11 23.70 1.07 7.58 1.06 11 0.46 42.02 1.12 23.82 1.08 7.63 1.06 10 0.42 42.20 1.12 23.93 1.09 7.67 1.07 9 0.38 42.44 1.13 24.08 1.09 7.72 1.07 8 0.33 42.65 1.14 24.21 1.10 7.78 1.08 7 0.29 42.94 1.15 24.40 1.11 7.85 1.09 6 0.25 43.25 1.16 24.60 1.12 7.92 1.09 5 0.21 43.69 1.18 24.87 1.13 8.03 1.11 96 4 0.17 44.36 1.20 25.27 1.15 8.17 1.12 3 0.13 45.22 1.23 25.79 1.17 8.35 1.14 2 0.08 46.82 1.27 26.70 1.21 8.64 1.18 1 0.04 50.12 1.35 28.52 1.29 9.18 1.27 Table B.4: Temperature, Nusselt, and Rayleigh data, five die model, FC-72 L = 24 mm ?? (w/cm2) 2 1 0.25 0.125 S/L = 1 ?T (K) 40.21 22.21 6.80 3.80 Nu 209.45 189.60 154.83 138.43 Ra 7.24E+07 4.00E+07 1.22E+07 6.85E+06 S/L = 0.5 ?T (K) 40.80 22.63 7.05 3.94 Nu 206.39 186.04 149.28 133.51 Ra 7.35E+07 4.07E+07 1.27E+07 7.10E+06 S/L = 0.25 ?T (K) 43.32 24.17 7.62 4.29 Nu 194.40 174.19 138.05 122.76 Ra 7.80E+07 4.35E+07 1.37E+07 7.72E+06 L = 18 mm ?? (w/cm2) 2 1 0.25 0.125 S/L = 1 ?T (K) 35.74 19.68 5.96 3.30 Nu 176.74 160.47 132.38 119.46 Ra 2.71E+07 1.49E+07 4.53E+06 2.51E+06 S/L = 0.5 ?T (K) 36.79 20.46 6.34 3.52 Nu 171.68 154.32 124.61 112.13 Ra 2.79E+07 1.55E+07 4.81E+06 2.67E+06 S/L = 0.25 ?T (K) 39.72 22.26 6.99 3.91 Nu 159.01 141.87 112.90 101.04 Ra 3.02E+07 1.69E+07 5.31E+06 2.97E+06 97 Table B.5: Temperature, Nusselt, and Rayleigh data, five die model, Novec 649 L = 24 mm ?? (W/cm2) 1 0.25 0.125 S/L = 1 ?T (K) 21.71 6.65 3.72 Nu 187.35 152.89 136.64 Ra 3.96E+07 1.21E+07 6.79E+06 S/L = 0.5 ?T (K) 22.11 6.89 3.85 Nu 184.00 147.53 131.92 Ra 4.04E+07 1.26E+07 7.04E+06 S/L = 0.25 ?T (K) 23.59 7.45 4.19 Nu 172.43 136.58 121.42 Ra 4.31E+07 1.36E+07 7.64E+06 L = 18 mm ?? (W/cm2) 1 0.25 0.125 S/L = 1 ?T (K) 19.26 5.84 3.24 Nu 158.44 130.69 117.84 Ra 1.48E+07 4.49E+06 2.49E+06 S/L = 0.5 ?T (K) 19.99 6.19 3.44 Nu 152.65 123.19 110.83 Ra 1.54E+07 4.77E+06 2.65E+06 S/L = 0.25 ?T (K) 21.72 6.83 3.82 Nu 140.45 111.72 99.93 Ra 1.67E+07 5.26E+06 2.94E+06 98 Table B.6: Temperature, Nusselt, and Rayleigh data, five die model, HFE 7100 L = 18 mm ?? (W/cm2) 2 1 0.25 0.125 S/L = 1 ?T (K) 33.03 18.22 5.55 3.09 Nu 157.98 143.19 117.57 105.44 Ra 2.21E+07 1.22E+07 3.71E+06 2.07E+06 Table B.7: Temperature, Nusselt, and Rayleigh data, four die model, FC-72 L = 24 mm ?? (W/cm2) 2 1 0.25 0.125 S/L = 1 ?T (K) 39.46 21.78 6.74 3.77 Nu 213.40 193.29 156.14 139.65 Ra 7.10E+07 3.92E+07 1.21E+07 6.79E+06 S/L = 0.5 ?T (K) 39.91 22.14 6.91 3.87 Nu 211.02 190.18 152.35 135.94 Ra 7.18E+07 3.99E+07 1.24E+07 6.97E+06 S/L = 0.25 ?T (K) 41.79 23.29 7.33 4.13 Nu 201.52 180.81 143.57 127.59 Ra 7.52E+07 4.19E+07 1.32E+07 7.43E+06 L = 18 mm ?? (W/cm2) 2 1 0.25 0.125 S/L = 1 ?T (K) 34.71 19.31 6.00 3.35 Nu 181.96 163.55 131.62 117.88 Ra 2.64E+07 1.47E+07 4.56E+06 2.54E+06 S/L = 0.5 ?T (K) 35.53 19.82 6.21 3.48 Nu 177.74 159.36 127.23 113.28 Ra 2.70E+07 1.51E+07 4.71E+06 2.65E+06 S/L = 0.25 ?T (K) 37.43 20.95 6.61 3.73 Nu 168.75 150.70 119.42 105.89 Ra 2.84E+07 1.59E+07 5.02E+06 2.83E+06 99 Table B.8: Temperature, Nusselt, and Rayleigh data, four die model, Novec 649 L = 24 mm ?? (W/cm2) 1 0.25 0.125 S/L = 1 ?T (K) 21.61 6.63 3.70 Nu 188.25 153.33 137.32 Ra 3.94E+07 1.21E+07 6.76E+06 S/L = 0.5 ?T (K) 21.85 6.82 3.81 Nu 186.13 149.16 133.38 Ra 3.99E+07 1.24E+07 6.96E+06 S/L = 0.25 ?T (K) 23.13 7.28 4.07 Nu 175.90 139.68 124.79 Ra 4.22E+07 1.33E+07 7.44E+06 L = 18 mm ?? (W/cm2) 1 0.25 0.125 S/L = 1 ?T (K) 19.10 5.84 3.24 Nu 159.73 130.66 117.56 Ra 1.47E+07 4.50E+06 2.50E+06 S/L = 0.5 ?T (K) 19.61 6.06 3.37 Nu 155.56 125.89 113.00 Ra 1.51E+07 4.67E+06 2.60E+06 S/L = 0.25 ?T (K) 21.01 6.55 3.66 Nu 145.19 116.43 104.15 Ra 1.62E+07 5.04E+06 2.82E+06 100 APPENDIX C ? TWO-PHASE RESULTS In this appendix, the results from the two-phase simulations are presented. Table C.1 contains the heat flux and temperature data that was plotted in Figures 4.6-4.8. Table C.2 contains the temperature and heat flux data that was reported in Figures 4.15, 4.16, and 4.19. Table C.1: Pool boiling simulation results from the multi-chip module model CASE A CASE B CASE C ?? (W/cm2) ?T (K) ?? (W/cm2) ?T (K) ?? (W/cm2) ?T (K) 3 22 3 31.6 4 38.15 5 24.1 4 33.95 6 41.2 6 24.9 6 38.2 8 42.65 7 25.7 8 41.15 10 44.05 7.5 26.1 9 42.7 12 44.8 8.5 26.8 12 46.4 9.5 28.5 10.5 30.8 12 34.9 Table C.2: Flow boiling simulation results from the cartridge model Numerical, without CWH Numerical, with CWH 800 mL/min 1200 mL/min 800 mL/min ?? (W/cm2) ?T (K) ?? (W/cm2) ?T (K) ?? (W/cm2) ?T (K) 12 53.25 12 51.85 6 43.55 10 51.55 10 50.45 4 46.85 8 49.85 8 48.55 6 47.45 6 46.15 4 44.15 4 42.6 101 APPENDIX D ? FLUID PROPERTIES In this appendix, the fluid properties used in all simulations is presented. Liquid properties are presented for all fluid used for single phase simulations, vapor properties are included for fluids used in two-phase simulations. Table D.1: Fluid properties used for simulations Liquid Properties ?? (kg/m3) ? (1/K) ??? (J/kg-K) ?? (W/m-K) ?? ( kg/m-s) ?? (m2/s) Pr hl (kJ/kg) FC-72 1680.0 1.56E-03 1100.00 0.06 6.40E-04 3.81E-07 12.35 NA Novec 649 1575.1 1.80E-03 1108.10 0.06 6.40E-04 4.06E-07 11.96 0 HFE 7100 1480.0 1.80E-03 1166.00 0.07 5.70E-03 3.85E-06 9.63 NA Water 998.2 4.60E-05 4182.00 0.60 1.00E-03 1.00E-06 2.2 0 Vapor Properties ?? (kg/m3) ??? (J/kg-K) ?? (W/m-K) ?? ( kg/m-s) ? (N/m) hv (kJ/kg) Novec 649 5.56 906.40 0.00582 1.84E-07 0.009222 88.0 Water 0.55 1563.08 0.0261 1.34E-05 0.0589 2259.2 102 APPENDIX E ? INTERFACIAL MOMENTUM TRANSFER MODELS FOR BOILING The following is a summary of the available interfacial momentum transfer models used with the RPI Boiling Model as described in the ANSYS Fluent 14.0 Theory Guide [14]. The coefficient for interfacial lift force (boiling-moraga) is given by the following correlation [32]: CLift = { 0.00767 ? > 6000 ?(0.12 ?0.2e? ? 3600)e ? 3x10 ?7 6000 ? ? ? 1.9x105 ?0.002 ? < 1.9x105 (F.1) where ? = RebRev, the lift coefficient combines the opposing action of two lift forces: the classical aerodynamics lift force resulting from the interaction between bubble and liquid shear, and the lateral force resulting from interaction between bubble vortexes shed by bubble wakes. Here Reb = ?lDb|ul?uv|? l is the bubble Reynolds Number, and Rev = Db 2?l(?ul) ?l is the bubble shear Reynolds Number. Here Db is the bubble diameter, which is a function of local subcooling, ?T = Tsat ?Tl. Db = { 0.00015 ?Tsub > 13.5K 0.00015 ? 0.00001?Tsub 0 ? ?Tsub ? 13.5K 0.0015 ?Tsub < 0.0 (F.2) The heat transfer from the bubble interface to the liquid phase (ranz-marshall model) is defined by the following: q?lt = Aihsl(Tsat ?Tl) (F.3) Ai is the interfacial area given by: Ai = 6(1?vf)min (vf,vf,crit)D b(1?min (vf,vf,crit)) (F.4) 103 where vf is the volume fraction of vapor, vf,crit is the critical volume fraction of vapor, vf,crit = 0.8. hsl is the heat transfer coefficient based on the Ranz-Marshall correlation [33], and is given by the following: hsl = klD b (2 +0.6Re0.5Pr0.33) (F.5)