COMPUTATIONS OF DROP DYNAMICS WITH HEAT TRANSFER By MARIANNE FRANOIS A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2002
Copyright 2002 by Marianne Franois
To my parents and my sister.
iv ACKNOWLEDGMENTS My appreciation and respect go to my advisor, Profe ssor Wei Shyy, for introducing me to the field of computational fluid dynamics with moving boundaries and for giving me the opportunity to study bubble and drop dynamics for my Ph.D . dissertation. His continuous support and patience he lped me to achieve this wo rk. I would like to express my special thanks to Professor Jacob Chung fo r his help and encouragement. I sincerely thank Professor Ulrich Kurzweg, Professor Re nwei Mei and Professor Corin Segal, for serving as members of my supervisory committee. I would like to thank Professor H. S. Udaykumar for helpful research discussion during my first year at the University of Florida. I acknowledge th e collaboration with Daniel Lrstad, from Lund University, Swed en, on the work presented in chapter 5, section 5.1.2. I would like to thank all my co lleagues of the comput ational thermo-fluid group, in particular, Narcisse Nâ€™Dri and Inan c Senocak. I also extend my thanks to the individuals in the Department who have helped me one way or another during my graduate studies. I would like to acknowledge the financial suppor t of the Air Force Office of Scientific Research (AFOSR) and National Aeronautics and Space Administration (NASA) for this research through grants awarded to my advisor. I also than k Zonta International for its financial help through the Amelia Earhart fe llowship I was awarde d for the academic year 2000-2001.
v Finally, I would like to thank my parents and my sister for their continuous support and encouragement over the years of my studies . To them, I dedicate this dissertation.
vi TABLE OF CONTENTS page ACKNOWLEDGMENTS.................................................................................................iv LIST OF TABLES...........................................................................................................viii LIST OF FIGURES.............................................................................................................x NOMENCLATURE..........................................................................................................xv ABSTRACT...................................................................................................................xxii i CHAPTERS 1 INTRODUCTION...........................................................................................................1 1.1 Motivation................................................................................................................. 1 1.2 Overview of Numerical Methods for Multiphase Flow............................................3 1.3 Objectives................................................................................................................. 5 1.4 Structure of the Dissertation.....................................................................................7 2 MICRO-SCALE THERMAL MANAGEMENT DEVICE...........................................10 2.1 Thermodynamic Analysis.......................................................................................10 2.2 Phase Change Analysis...........................................................................................18 2.3 Analysis of the Ejection Mechanism......................................................................22 3 LITERATURE REVIEW..............................................................................................25 3.1 Droplet Ejection......................................................................................................25 3.1.1 Droplet Ejector Devices a nd Physical Observations...................................25 3.1.2 Simulation of Droplet Ejection....................................................................29 3.2 Droplet Impingement..............................................................................................38 3.2.1 Experimental Observations of Droplet Impingement..................................38 3.2.2 Simulations of Droplet Impingement...........................................................44 3.3 Condensation and Evaporation of a Single Drop on a Flat Surface.......................49 3.3.1 Dropwise Condensation...............................................................................51 3.3.2 Dropwise Evaporation.................................................................................59
vii 4 THEORETICAL FORMULATION A ND COMPUTATIONAL METHODS............67 4.1 Theoretical Formulation.........................................................................................67 4.1.1 Validity of Continuum Model......................................................................67 4.1.2 Governing Equations...................................................................................68 4.1.3 Interfacial Boundary Conditions..................................................................69 4.1.4 Formulation with an Imme rsed Boundary Interface....................................72 4.1.5 Nondimensionalization of the Governing Equations...................................73 4.1.6 Dimensionless Parameters within the Micro-Scale Cooling Device...........75 4.2 Computational Methods..........................................................................................77 4.2.1 Fluid Flow Field Equation Solver................................................................78 4.2.2 Immersed Boundary Technique...................................................................95 4.3 Summary...............................................................................................................107 5 DYNAMICS OF RISING DROP AND BUBBLE......................................................109 5.1 Rising Bubble and Drop by Buoyancy.................................................................109 5.1.1 Rising Bubble in a Viscous Liquid............................................................110 5.1.2 Rising Oil Drop in Water...........................................................................120 5.2 Drop Motion within the Micro-Scale Cooling Device.........................................126 5.3 Summary...............................................................................................................130 6 DROP IMPACT DYNAMICS WITH HEAT TRANSFER........................................131 6.1 Initial and Boundary Conditions...........................................................................131 6.2 Results with a Static Contact Angle Model..........................................................134 6.2.1 Impact of Water and Ink Dr op on Polycarbonate Surface.........................134 6.2.2 Parametric Variation Effects on the Drop Dynamics................................136 6.3 Results with a Dynamic Contact Angle Model....................................................149 6.3.1 Static versus Dynamic Contact Angle.......................................................149 6.3.2 Comparison with Experiment....................................................................155 6.4 Heat Transfer Analysis.........................................................................................161 6.4.1 Heat Transfer for Re=300, We=60 and =60 Case..................................161 6.4.2 Parametric Variation Effects on the Heat Transfer....................................165 7 SUMMARY AND FUTURE WORK.........................................................................174 7.1 Summary and Conclusions...................................................................................174 7.2 Future Work..........................................................................................................177 LIST OF REFERENCES.................................................................................................178 BIOGRAPHICAL SKETCH...........................................................................................190
viii LIST OF TABLES Table page 2.1 Ideal cycle characteristics for T3=50 C, T1=20 C..................................................15 2.2 Cycle with a compressor efficiency of 0.8..............................................................16 2.3 Effect of irreversib ility on the coefficient of performance for R134a....................16 3.1 Summary of the numerical si mulations of droplet ejection....................................37 3.2 Relevant numerical studies of a wate r liquid droplet impac ting on a flat solid surface.....................................................................................................................49 3.3 Summary of the previous analytical studies in single dr oplet condensation..........58 3.4 Classification and description of the dynamics and heat tran sfer of a droplet evaporation..............................................................................................................60 3.5 Previous studies on low s uperheat dropwise evaporation.......................................66 4.1 Properties of saturated liquid-va por at 1 atm from Carey (1992)...........................76 4.2 Liquid to vapor ratio of the fluid prop erties for water and R-12 at 1 atm and at saturation temperature.............................................................................................76 4.3 Range of dimensionless parameters in droplet ejection/impingement problem.....76 4.4 Multigrid performance for a cavity lid driven flow for Re=100.............................88 4.5 Natural convection results for Ra=104, 105, 106 compared to benchmark results..92 4.6 Effect of the density ratios for a spherical drop in st atic equilibrium.....................105 5.1 Bubble aspect ratio at t=1.5 as a function of Re and We........................................119 5.2 Comparison of drag coefficients as a function of Re and We................................120 5.3 Physical properties of S.A.E. 10W oil drops and water at T=23oC from Klee and Treybal (1956)........................................................................................................120 5.4 Computational cases for the rising oil drop in water..............................................121
ix 6.1 Effect of the grid size and backgr ound fluid, i.e. gas properties on the drop spreading coefficient at t= 1.0 and 5.0 for Re=200, We=20 and =60 ..................137 6.2 Computational fluid properties...............................................................................150 6.3 Contact angle values of water on a wax coated surface from Fukai et al. (1995)..150 6.4 Cases considered for the comparison between the present numerical method and experimental results................................................................................................155 6.5 Computational thermal properties...........................................................................162
x LIST OF FIGURES Figure page 2.1 Schematic of the reversed Rankine cycle or vapor-compression refrigeration cycle........................................................................................................................12 2.2 Temperature-entropy (T-s) diagram and pressure-enthalpy (P-h) diagram of the ideal vapor compression refrigeration cycle...........................................................12 2.3 Schematics of the conceptual micro-scale cooling device......................................14 2.4 Coefficient of performance (COP) and compression ratio (P2/P1) versus T=(T3-T1) with T3 =constant=50 C for ideal cycle..............................................17 2.5 3-D schematic perspective of a liquid droplet on a solid surface...........................18 2.6 Evaporation/condensation time scale versus droplet radius...................................21 2.7 Schematic of the three st ages in the droplet motion...............................................23 2.8 Total energy requirement versus drop let size for R-12 with U=1m/s and H=1cm.23 2.9 Effect of gravity on en ergy required for ejection vers us droplet radius for R-12 with U=1m/s and H=1cm........................................................................................24 3.1 Illustration of different droplet ejector devices.......................................................26 3.2 Ejection of ethylene glycol droplets produced by a piezoelectric inkjet type device from Shield, Bogy and Talke (1987)...........................................................27 3.3 Rapid atomization of a 0.1 ml water droplet forced by st epped actuation at 1080Hz, from Vukasinovic, Glezer and Smith (2000)...........................................29 3.4 Instantaneous shapes, str eamlines and pressure contour s of droplet ejection from the numerical simulation re sults of Fromm (1984)................................................31 3.5 Drop shape, length (L), and rod vertical di splacement (zâ€™) versus time t for Re20 We , from Wilkes and Basaran (2001)...................................................35 3.6 Pendant drop undergoing oscillations for Re20 We , from Wilkes and Basaran (2001)........................................................................................................36
xi 3.7 Schematic of a sessile droplet on a flat surface......................................................39 3.8 The impact of a n-heptane droplet on a heated stainless stee l surface (view at an angle of 30 ) for Re=2300 and We=43 from Chandra and Avedisian (1991)........43 3.9 Numerical simulation results from Zhao et al. (1996a) of the fluid dynamics and heat transfer phenomenon during the im pingement of a water droplet for Re=500, We=10 and Fr=630..................................................................................48 3.10 Cycle of dropwise condensation (adapted from Tanasawa, 1991). ........................52 3.11 Schematic of the droplet evapor ation lifetime and boiling curves. ........................59 4.1 Schematic of a computational cell arrangement.....................................................79 4.2 Illustration of the discretization notations...............................................................81 4.3 Streamlines plot at steady state of the cavity lid driven flow for Re=100 on a 102 x 102 grid.........................................................................................................87 4.4 Convergence path of the pressure equation for a cavity lid driven flow for Re=100 at the first time step...................................................................................87 4.5 Illustration of the non-uniform grid 114 x 114.......................................................91 4.6 Natural convection patterns for Pr=0.71 obtained with a 114 x 114 grid...............91 4.7 Velocity profile for Poisseuille flow, comparison of exact solution with numerical solution...................................................................................................94 4.8 Error 2-norm of the radial velocity versus grid spacing for Poisseuille flow.........95 4.9 Illustration of a comput ational domain composed of two immiscible fluids..........96 4.10 Discrete Heaviside step function in 1D for three diffe rent grid spacing with the interface located at x=0.5........................................................................................98 4.11 Discrete delta function in 1D for thr ee different grid spaci ng with the interface located at x=0.5.......................................................................................................100 4.12 Marker points considered for the estimation of the force at point P.......................101 4.13 Fluid points considered for the evaluati on of the velocity at marker point X........102 4.14 Effect of grid refineme nt on the error for the validation test case of a droplet in static equilibrium....................................................................................................106 4.15 Pressure contours and velocity field of a drop in static equilibrium. drop/fluid=10, t=10-4, R/x=10.............................................................................106
xii 4.16 Illustration of one leve l local grid refinement.........................................................108 5.1 Schematic of the computational se tup for a single bubble rising by buoyancy......112 5.2 Multigrid performance on the computations of a bubble rising by buoyancy........113 5.3 Rising speed and aspect ratio versus time for Re =100, We=4, Fr=1, density ratio=100, viscosity ratio=1, on 202 x 42 (d/h=20) and 402 x 72 grids (d/h= 40).113 5.4 Effect of the density ratio (d.r.) on the convergence path of the pressure equation for Re=100, We=4, Fr=1 at the first time iteration on 402 x 72..............114 5.5 Instantaneous bubble shapes at time t=0, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5 for Re=100, We=4 and Fr=1 for density ra tios of 10, 100 and 1000, and viscosity ratio of 1..................................................................................................................115 5.6 Instantaneous bubble shapes at time t=0, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5 for Re=100, We=4 and Fr=1 for viscosity ra tios of 1, 10 and 100, and density ratio of 10........................................................................................................................11 6 5.7 Effect of the density ratio on the risi ng speed and aspect ratio versus time for Re=100, We=4, Fr=1, viscosity ratio=1.................................................................116 5.8 Effect of the viscosity ratio on the rising speed and as pect ratio versus time for Re=100, We=4, Fr=1, density ratio=10..................................................................117 5.9 Streamlines in a reference frame a ttached to the bubble. Re=100, We=4 and Fr=1 for the three density ratios cases at t=1.5.......................................................118 5.10 Computed rising bubble shap es at time t=1.5 as a function of Re and We with a density ratio of 100 and viscosity ratio of 1............................................................119 5.11 Effect of wall proximity on instan taneous rising drop sh apes at the same nondimensional time t=0, 2, 4, 6, 8, 10 of a 5 mm oil drop in water......................122 5.12 Effect of wall proximity on rising speed in percentage..........................................122 5.13 Drop rising speed and asp ect ratio versus time for d=5 mm..................................123 5.14 Effect of grid refine ment on rising speed and as pect ratio at nondimensional time t=10.................................................................................................................123 5.15 Effect of the drop diam eter on the instantaneous shap es of a rising oil drop in water at t=0, 2, 4, 6, 8, 10, 12.................................................................................125 5.16 Rising speed and aspect ratio versus drop diameter...............................................126 5.17 Schematic of the computationa l setup for the droplet dynamics............................127
xiii 5.18 Effect of droplet diam eter on droplet dynamics.....................................................128 5.19 Effect of initial velocity on droplet dynamics.........................................................128 5.20 Instantaneous shapes for a R-12 droplet.................................................................129 5.21 Instantaneous shapes for water droplets.................................................................129 6.1 Initial and boundary conditions for th e droplet impingement problem with heat transfer....................................................................................................................132 6.2 Water droplet impinging on polycarb onate surface for Re=3200 and We=30.......135 6.3 Ink droplet impingi ng on polycarbonate surface for Re=2300 and We=190.........136 6.4 Droplet shapes at different time inst ants for varying static contact angle..............141 6.5 Droplet shapes at different time in stants for varying impact velocity....................142 6.6 Droplet shapes at different time in stants for varying droplet viscosity..................143 6.7 Droplet shapes at different time instants for varying surface tension.....................144 6.8 Parametric variation effects on th e spreading coefficient versus time...................145 6.9 Comparison between the values of th e maximum spreading coefficient obtained from the present simulation and from th e energy correlation of Eq. (3.13)...........146 6.10 Pressure contours and streamlines at di fferent time instants for Re=100, We=20, =60 .......................................................................................................................148 6.11 Illustration of the dyna mic contact angle model.....................................................150 6.12 Effect of grid refinement on spre ading coefficient and maximum height for Re=300, We=60 with the static contact angle model.............................................152 6.13 Effect of grid refinement on spre ading coefficient and maximum height for Re=300, We=60 with the dynami c contact angle model........................................152 6.14 Instantaneous drop shapes at same inst ant for the static contact angle model (left column) and dynamic contact angle mode l (right column) for Re=300, We=60 and a grid with d/h=40............................................................................................153 6.15 Spreading coefficient and normalized maximum height versus dimensionless time t* for Re=300, We=60 with a stat ic contact angle model and a dynamic contact angle model and a grid with d/h=40...........................................................154 6.16 Contact angle versus dimensionless tim e t* for Re=300, We=60 and a grid with d/h=40.....................................................................................................................154
xiv 6.17 Comparison of the instantaneous drop sh apes of the present numerical results with the numerical and experimental resu lts of Pasandideh-Fard et al. (2001) of a water drop impinging on stainless steel...............................................................157 6.18 Comparison of the spreading coefficien t of present numerical results with the experimental results of Pasa ndideh-Fard et al. (2001)............................................158 6.19 Contact angle versus dimensionless c ontact line velocity a nd versus time during the course of the present computation for Re=2900, We=47 case of PasandidehFard et al. (2001).....................................................................................................158 6.20 Comparison of the present numerical results with the experimental and numerical results of Fuka i et al. (1995)..................................................................160 6.21 Contact angle versus dimensionless c ontact line velocity a nd versus time during the course of the present computation fo r Re=6200, We=128 case of Fukai et al. (1995)......................................................................................................................161 6.22 Instantaneous droplet shapes and isothermals every T=0.1 during droplet impact on a constant temperature flat surface for Re=300, We=60, =60 ...........164 6.23 Wall heat flux (q) versus radial position (r) at different time instants for Re=300, We=60, =60 ..........................................................................................165 6.24 Instantaneous droplet shapes and isothermals every T=0.1 during droplet impact on a constant temperature flat surface for Re=100, We=20, =60 ...........166 6.25 Drop shapes and corresponding wall heat flux (q) versus radial position (r) at different time instants for Re=100, We=20, =60 ................................................167 6.26 Effect of the contac t angle on wall heat flux..........................................................169 6.27 Effect of the impact velocity on wall heat flux.......................................................170 6.28 Effect of the droplet viscosity on wall heat flux.....................................................171 6.29 Effect of the surface tension on wall heat flux.......................................................172 6.30 Parametric effects on the Nusselt number..............................................................173
xv NOMENCLATURE A area C curve representing the interface CD drag coefficient Cp specific heat COP coefficient of performance D diameter of computational cylinder in chapter 5 and 6 D diffusion terms in chapter 4 d drop diameter d distance equal to twice the grid spacing d=2h dbase basal diameter of drop in contact with surface after impact dmol. diameter of a molecule Ein energy input 2dE dissipation energy 1k E kinetic energy before impact 1 s E, 2 s E surface energy before and after impact, respectively F convective flux in chapter 4 Fd drag force fN natural frequency g gravitational acceleration
xvi H height between condenser surface and evaporator surface in chapter 2 H height of comput ational cylinder in chapter 5 and 6 h enthalpy in chapter 2 h grid spacing in chapter 4, 5, and 6 h heat transfer coefficient in chapter 3 I indicator function k thermal conductivity L characteristic length scale M molecular weight M mass of object (drop or bubble) m mass m mass flow rate m.f.p. mean free path N Avogadroâ€™s constant Nx, Ny number of grid points in the axia l and radial direction, respectively Nin, Nout number of computational grid cells in and out of drop, respectively Nu Nusselt number n unit normal vector P pressure p pressure q heat flux q heat flux Q total heat
xvii Q heat transfer rate H Q heat rejected L Q heat absorbed or cooling load r, x axisymmetric coordinate R radius R1, R2 principal radii of curvature Rh radius of hemispherical drop Rd radius of spherical drop 32dhRR S source term s entropy in chapter 2 s arclength T temperature t time t* dimensionless time, t*=tU/d t unit tangential vector tc condensation time scale te evaporation time scale tsp. drop spreading time scale U drop vertical velocity U drop impact speed U bubble rising speed in an infinite domain u bubble rising speed u velocity vector
xviii V interface velocity vector V drop volume Wd drop weight inW work input x, y Cartesian coordinate xi vapor mole fraction at the liquid-air interface xa air-vapor mole fraction Greek Symbols thermal diffusivity spreading coefficient coefficient of volume expansion ratio of specific heat any material property incremental difference thickness of disk shape drop after impact Dirac delta function liquid to solid conductivity ratio diffusion coeffiecient dimensionless maximum drop height contact angle interface curvature
xix latent heat dynamic viscosity kinematic viscosity number pi density surface tension viscous stress any field variable vis dissipation function streamfunction volume interpolation weight Superscripts 0 initial g grid level n time indice Subscripts 0 initial 1 fluid 1 2 fluid 2 1,2,3,4 vapor-compression cycle stages
xx E,W,N,S neighboring cell centered points e,w,n,s face cell points i,j cyclical indices A advancing C constant G finest grid level g grid level g gas H high k interfacial marker indice L low l liquid M indice for the four surrounding points max maximum min minimum n cell indice P cell centered point R receding r radial direction S static or at equilibrium s solid s first derivative w ith respect to arclength ss second derivative with respect to arclength
xxi sat saturation v vapor w wall x axial direction Nondimensional numbers Bi Biot number Bi=hr/kl Bo Bond number Bo=We/Fr Ca Capillary number Ca=We/Re Eo Etvs number Eo=g d2/ Fr Froude number Fr=U2/(gd) Ja Jakob number Ja=Cp T/ Kn Knudsen number Kn=m.f.p./L Mo Morton number Mo=g 4 /( 23) Oh Ohnesorge number Oh=We/Re Pe Peclet number Pe=Re.Pr Pr Prandtl number Pr= Cp/k Ra Rayleigh number Ra= g (Thot-Tcold)L3/( ) Re Reynolds number Re= Ud/ R/W R/W=1/Oh We Weber number We= U2d/ Gibbs-Thomson number = Tsat/( T l d)
xxii Operators gradient . divergence
xxiii Abstract of Dissertation Pres ented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy COMPUTATIONS OF DROP DYNAMICS WITH HEAT TRANSFER By Marianne Franois August 2002 Chair: Dr. Wei Shyy Major Department: Aerospace Engineeri ng, Mechanics and Engineering Science In computer, electronics and avionics areas , with continuous re duction in physical dimensions and increasing power density of the devices, therma l management is a critical issue. To help address this need, a microscale cooling device concept based on dropwise condensation and spray cooling is proposed. Compared to alternative singleand multiphase heat transfer techniques, dropwise phase change offers hi gher heat transfer rates. In the present research, thermodyna mic cycle and component analys es are first presented to demonstrate the performance of the proposed device. Then, a comprehensive literature review is made on each main physical mech anism involved in the micro-scale device, namely, dropwise condensation and evaporati on, drop ejection, a nd drop impingement. The drop dynamics related to impingement and heat transfer is investigated based on detailed multiphase fluid flow computations. To facilitate the study, the immersed boundary method is adopted to handle th e moving boundary dynamics, where the interface is represented by marker points and advected in a Lagrangian framework. The
xxiv mass, momentum, and energy c onservation equations are solv ed on a fixed (Eulerian) Cartesian grid using a s econd-order projection method along with the multigrid technique. Numerous singleand multi-phase fluid flow problems are computed and assessed against the published e xperimental, analytical and computational information. In particular, the impacts of la rge property ratios between pha ses and the grid resolution on the computational efficiency and accuracy are evaluated quantitatively. In simulating the drop impingement, both static and dynamic cont act angle models are adopted. It is shown that high heat transfer rates occur during the early stage of im pingement, and maximum wall heat fluxes are located near the contac t line. The effects of inertia, wettability, viscosity and surface tension on drop spreading, recoiling, as well as heat transfer characteristics are also highlighted.
1 CHAPTER 1 1 INTRODUCTION 1.1 Motivation There is a growing interest in micro-sc ale fluid dynamics and heat transfer for applications including flow control, flight vehicle prot ection, micro-propulsion, and thermal management. This interest is mo tivated by recent advanc es in micro-device fabrication, computational t echniques, and diagnostic tools, which have enabled us to explore the micro-scale multiphase physics an d engineering practic es. New designs in sensors, actuators, and micr o devices such as valves, pumps, turbines, and heatexchangers are being adopted in propulsion, vehicles, and energy sectors (Ho and Tai, 1998; Ldfdahl and Gad-El-Hak, 1999; Bis hop et al. 2001; Sentur ia, 2001; Tien, 1997). In computer, electronics and avionics areas, with con tinuous reducti on in physical dimensions and increasing power density of the devices, therma l management is a critical issue (Semiconductor Industry Association [SIA], 1999; Sarno and Moulin, 2001). Cooling systems with a higher rate of heat tr ansfer than traditional techniques such as fans and fin arrays are essential to support fu rther development in these critical technical areas. Since the two-phase flow with phase-chang e offers an efficient way of transferring more energy across a smaller temperature grad ient, micro-scale syst em designs utilizing condensation and evaporation, like micro heat pipes (Pet erson et al., 1997), vibrationinduced droplet atomization heat transfer cell (Heffington et al ., 2002) and miniature â€œrefrigeratorsâ€ (Shannon et al., 1999; Ashraf et al., 1999), are being actively pursued.
2 In multiphase heat transfer, both filmwi se and dropwise mechanisms have been explored (Collier and Thome, 1994). Despite the fact that many investigators have evaluated general drop condensation and grow th characteristics (e.g., Peterson et al., 1992; Majumdar and Mezic, 1999), its mechanis m is not so well understood as to support robust engineering design under varying opera ting conditions. Therefore, even though dropwise condensation offers h eat transfer rates at least an order of magnitude higher than filmwise condensation, th e latter concept has been more frequently used in industrial applications (Collier and Thome, 1994). Ho wever, even if one accepts a lower heat transfer rate, certain operating conditions are not suitable for filmwi se condensation. For example, for spacecraft under much reduced gravity, filmwise condensation becomes unstable, and cannot function as a continuous process due to a lack of condensate removal mechanism. Hence, in a microgravity environment, it seems more advantageous to promote dropwise condensation than to pr omote filmwise condensation. Furthermore, as will be demonstrated later, condensation and ejection of drops with diameters up to few hundreds micron-meters are not sensitive to the gravity level, making both design and testing of microgravity a pplications easier to conduct in an earthbound environment. Besides microgravity applicati ons, a highly effective mode of heat transfer is spray cooling (Liu et al., 1999; Ya ng et al., 1996). Ghodbane and Holman (1991) and Holman and Kendall (1993) have reported heat fluxes in excess of 6 x 105 W/m2 for horizontal spray of Freon-113. It seems clear that a device combining dropwise condensation and spray cooling can be very attractive for ther mal management applications for both earthbound and microgravity applicatio ns. In this work, a novel conceptual device based on dropwise condensation and evaporation, two hi ghly efficient heat transfer modes, is
3 proposed. With the growing development of mu ltiphase micro-scale devices, there is an imperative need to develop analytical and num erical models not onl y to help design but also to understand the associated physics. Fo r liquid flow in microdevice for drop having a diameter greater than 20 m, a continuum approximation is valid. At the continuum level, the governing equations of a thermo-f luid problem are the conservation equations of mass, momentum and ener gy. These governing equations are derived starting with Newtonâ€™s law and the thermodynamic laws. The thermo-fluid dynamic theoretical formulation of two-phase flow is discussed in Delhaye (1974) and Ishii (1975). Several numerical methods have been developed ov er the past 20 years to solve the mass, momentum and energy conservations equa tions involving an interface be tween two fluids. These methods are presented next. 1.2 Overview of Numerical Methods for Multiphase Flow Several numerical techniques for simulating multiphase physi cs have been proposed in the literature, each with its own advantages and disadvantag es (Crank, 1984; Floryan and Rasmussen, 1989; Shyy et al., 1996; Shyy et al., 2001). These techniques have been categorized as surface tracking or Lagrangi an methods, volume tracking or Eulerian methods, and combined Eulerian-Lagrangian methods (Shyy et al., 1996; Shyy et al., 2001). In the Lagrangian methods, the grid is co ntinuously configured to conform to the shape of the interface. The tw o main advantages of Lagrangian methods are (i) the interface is explic itly tracked, and (ii) the interfacial boundary conditions apply at the exact location of the interface. The main disadvantages with Lagrangian methods are the issues related to mesh distortions, geomet ric conservation law, and topological changes
4 (Shyy et al., 1996). As discussed later in chap ter 3, the Lagrangian techniques have been employed in drop dynamics by Wilkes et al. ( 1999), Wilkes and Basaran (2001), Fukai et al. (1993, 1995, 1997) and Zhao et al. (1996a). The Eulerian methods are base d on a fixed grid formulati on, and hence the interface between the two phases is not explicitly tracked but is reconstructed from the properties of appropriate field variables, such as fluid fractions. Approp riate procedures need to be developed to deduce the interface location based on the volume fraction information. Also, the interfacial boundary conditions need to be manipulated to cast in the governing transport equations. For example, in the continuum surface force (CSF) method of Brackbill et al. (1992), the surface tension forces are treated as a body force into the momentum equations. A main advantage of th e Eulerian approach is the capability of computing the field variables without having to deal with the grid generation issue. The principal disadvantages of Eu lerian methods are (i) the difficulty in calculating the position, and associated geometric quantitie s such as curvatures, of the interface accurately, and (ii) the possible smearin g of boundary information due to the manipulation of the interfacial boundary conditions (Shyy et al., 199 6). The Eulerian methods include some well adopted technique s such as the volum e of fluid method (VOF) (Hirt and Nichols, 1981; Kothe and Mj olsness, 1992; Francois , 1998; Scardovelli and Zaleski, 1999; Gueyffier et al., 1999) and the level-set method (Osher and Sethian, 1988; Sethian, 1996; Sussman et al., 1994; Sussman and Smereka, 1997). Methods combining the above two methods , referred to as the mixed EulerianLagrangian methods, track the interface explic itly on a fixed grid by means of marker particles. In this category, there are th e immersed boundary me thod (Peskin, 1977;
5 Udaykumar et al., 1997) or the front track ing method (Unverdi and Tryggvason, 1992; Juric and Tryggvason, 1998; Tryggvason et al., 2001), and the cut-cell method (Udaykumar et al., 1999; Ye et al., 1999; Ye et al., 2001). In the immersed boundary or front tracking methods a single set of govern ing equations for the two phases, e.g., the liquid and gas phases, is solved, and the surf ace tension force is cast into the momentum equation as a body force (Peskin, 1977; Brackb ill et al., 1992). The interface is of finite thickness within which the material propertie s vary smoothly. In the cut-cell method, individual governing equations are solved in each phase, and the interfacial boundary conditions are applied explicitly. There are no universally superior methods in efficiency, robustness and accuracy. All of the above techniques have been successf ully applied to diverse multiphase flow problems, as will be reviewed in ch apter 3 for the study of drop dynamics. 1.3 Objectives The main objectives of the present research are to propose a c onceptual design of a micro-scale cooling device util izing dropwise modes of heat transfer and to develop a numerical method able to simulate accurately drop dynamics with heat transfer relevant to the present micro-scale cooling device. Specifically, the objectives are to present a micro-scale device and perform a simplified analysis to examine the feasibility of the device develop a stable, fast and accurate numerical method based on the immersed boundary method to simulate unsteady two-dimens ional two-phase ax isymmetric flow validate the numerical met hod by studying benchmark cases for both the single and two phase flow solver and by comparing th e numerical solutions to experimental results investigate the accuracy and convergence charac teristics for two-phase flow with large property ratios across the phase boundary
6 study the dynamics of isothermal rising bubble and drop in an enclosed domain investigate the effect of in ertia, wettability, surface tensi on, and viscosity on the study of the fluid dynamics and heat transfer of an impinging drop on a flat surface held at constant temperature relevant to the micro-scale device The numerical methodology adopted herein is the immersed boundary technique, originally developed by Peskin (1977) and later extended by Unverdi and Tryggvason (1992), Udaykumar et al. (1997) and Kan et al. (1998). The immersed boundary technique is a mixed Eulerian-Lagrangian met hod, in the sense that the flow is solved using an Eulerian approach on a fixed Cartes ian mesh, and the interface is represented by a discrete set of points and adve cted in a Lagrangian way. In this method, a single set of conservation equations valid for both phase s is solved. The inte rface conditions are incorporated into the gover ning field equation (momentum and energy conservation laws) as source terms. As a result, the interface is considered to be of small non-zero thickness within which the values of the properties change smoothl y. This method is well suitable for solving bubble and drop dynamics that often encounter larg e deformation as has been proved in the existing litera ture. Among the published stud ies, Unverdi and Tryggvason (1992) have computationally investigated the rising of bubble by buoyancy, Kan et al. (1998) the bio-cell dynamics of leukocyt e, and Juric and Tryggvason (1998) boiling flows. This method is chosen for the following reasons: it can track the interface explicitly with marker particles on a fixed mesh and can handle large interfacial deformation. it is simple to implement and less dema nding in computational operations than a Lagrangian method (Dandy and Leal, 1989) or a fixed-grid sharp interface approach such as the cut-cell method (Ye et al., 1999, 2001). To our knowledge, the immersed boundary met hod is not known to have been adopted to investigate the problem of drop impingement with heat transfer. Previous studies in
7 drop impingement have employed a Lagrangian approach (Fukai et al., 1993, 1995, 1997; Zhao et al., 1996a) and Eulerian method such as VOF (Pasandideh-Fard et al. 1996, 2001; Bussmann et al., 1999, 2000). Furthermore in these studies onl y the liquid drop dynamics is modeled. In this work, the interf ace is treated as immersed, the liquid and gas interactions are considered, and the large pr operty jumps typical of water-steam fluid are accounted for. Special treatment is necessary at the contact line, i. e., where the interface is in contact with the solid surface. This re search will contribute to the understanding of the drop dynamics with and without heat tr ansfer, as well as the development of the immersed boundary method for drop impingement . The present simula tion will directly help refine the design of the micro-scale cooling device that uses an actuated membrane condenser. 1.4 Structure of the Dissertation The novel conceptual micro-scale cooling devi ce is presented and analyzed in chapter 2. The micro-scale device utilizes both dr opwise condensation and evaporation regimes, two highly efficient modes of heat tran sfer. Background and performance evaluation based on a simplified analysis of this mi cro-scale cooling device is offered. A thermodynamic analysis and a component an alysis are performed and serve as a feasibility study. For the component analysis a single droplet system where an artificial nucleation site is produced at the center of the membrane is considered. The physical mechanisms involved in this micro-scale devi ce are (1) the condens ation and evaporation of a droplet on a surface, (2) the ejection of a droplet by an actuated membrane, (3) the interaction of the droplet with the surr ounding vapor after detachment, and (4) the impingement of a droplet on a flat surface.
8 In chapter 3, a literature review is made on each main physical mechanism involved in the micro-scale device, namely, dropwise c ondensation and evaporation, drop ejection and drop impingement. The mechanisms of dropwise condensation and evaporation are described. Existing analytical models on dr opwise phase change are reviewed. Existing technologies of drop ejector are summarized. Experimental observations and previous computations of drop ejection and drop impingement are presented. The theoretical formulation and computational methods ar e discussed in chapter 4. The assumptions, governing equations and boundary conditions at an interface are presented. The theoretical formulation involvi ng an immersed interface is given. The fluid flow solver, based on the projecti on method (Chorin, 1968; Ye et al. 1999), is presented along with the discretization tec hniques. A multigrid technique (Shyy, 1994; Udaykumar, 2001) is adopted to accelerate th e rate of convergence of the pressure Poisson equation. The buoyancy driven fl ow inside a square cavity caused by differentiated heating is studied to validate th e single-phase thermo-fluid flow solver. The Poisseuille flow is investigated to validate the axisymmetric formulation. The immersed boundary technique (Peskin, 1977; Udaykumar et al., 1997) is then presented. The test case of a spherical drop in static equilibrium is performed to assess the performance and accuracy of the multiphase flow solver. In chapter 5, the dynamics of drop a nd bubble rising by buoyancy is studied for varying flow regimes and compared to previo us numerical and experimental studies in order to further validate the tw o-phase flow solver. The dyna mics of a rising drop with characteristics relevant to drops in the micro-device is investigated.
9 The fluid dynamics and heat transfer of an impinging drop on a flat solid surface held at constant temperature are st udied in chapter 6. Both st atic and dynamic contact angle models are implemented. Computations co rresponding to published experimental cases are performed. The effects of inertia, surface tension, liqui d wettability a nd viscosity on the drop dynamics with heat transfer are examined. This dissertation ends with chapter 7, whic h contains a summary and conclusion of the contributions of the present work, and reco mmendations for future research efforts.
10 CHAPTER 2 2 MICRO-SCALE THERMAL MANAGEMENT DEVICE In this chapter, an original micro-scale device based on the refrigeration cycle and utilizing dropwise phase change is proposed . Recently, design concep ts of miniaturized cooling systems have been proposed based on the refrigeration vapor-compression cycle (Shannon et al., 1999; Ashraf et al., 1999). In particular, Sha nnon et al. (1999) have used an electrostatic diaphragm with valves to perform compression, whereas Ashraf et al. (1999) have used a centrifugal compressor. In both studies, conventio nal cross-flow heat exchangers are employed as condenser and evaporator. In the proposed micro-scale cooling device, an actuated membrane is adopted as the condensing surface, i.e., condenser, as well as the ejec ting device. Therefore the drop lets ejected serve the dual purpose of maintaining dropwise condensation and creating a spray for highly efficient cooling. The key elements of the system are the use of membranes for compression, condensation, and droplet ejection. The feas ibility study of this micro-scale cooling device consists of (a) a thermodynamic analysis (b) a phase change analysis, (c) a component analysis. For the component anal ysis a single droplet system where an artificial nucleation site is produced at the center of the membrane is considered. 2.1 Thermodynamic Analysis The thermodynamic analysis is performed based on a refrigeration system utilizing the reversed Rankine vapor-compression cycle. Th e vapor-compression re frigeration cycle is shown in Figure 2.1 and their associated temp erature-entropy (T-s) diagram and pressureenthalpy (P-h) diagram are shown in Figure 2.2. The ideal vapor-compression cooling
11 system consists of four processes (Cenge l and Boles, 1998): 1-2 isentropic compression in a compressor, 2-3 constant pressure heat rejection in a conde nser, 3-4 isenthalpic throttling in an expansion devi ce, 4-1 constant pressure heat absorption in an evaporator. The coefficient of performance denoted COP is defined as the ratio of the cooling load ( L Q) to the work input (inW): 14 21 L inmhh Q COP Wmhh (2.1) where m is the mass flow rate and h the enthalpy. For reference purposes, the maximum possible COP is that of the Carnot refrigeration cycle, for which H H L LTQ TQ (2.2) thus, 1 1ideal H LCOP T T (2.3) where TH is the temperature at the high temperature sink ( T3) and TL is the temperature at the low temperature source ( T1).
12 Figure 2.1: Schematic of the reversed Ranki ne cycle or vapor-compression refrigeration cycle. Figure 2.2: Temperature-entropy (T-s) diagra m and pressure-enthal py (P-h) diagram of the ideal vapor compression refrigeration cycl e. (a) Temperature-entropy (T-s) diagram, (b) pressure-enthalpy (P-h) diagram. The cycle of the proposed micro-scale cooli ng device is shown in Figure 2.3. It is composed of four processes and of four distinct components: compressor, condenser, evaporator and expansion orifi ce. In this system a vapor-c ompression refrigeration cycle is followed. The saturated vapor is first compressed via membrane actuation. Following Condenser Evaporator Compressor Expansion Valve 1 2 3 4 HQ inW LQ(a) (b) T-s diagramP-h diagram Saturated liquid line Saturated vapor line P=const Compressed liquid region T s 1 3 2 4 Superheated vapor region Saturated liquid+vapor region Critical point Critical point Saturated liquid line Saturated vapor line Superheated vapor region P h 1 3 2 4 Compressed liquid region Saturated liquid+vapor region T=const
13 the compression process, the superheated vapor enters the condenser. The condenser surface is also a membrane but made of hi ghly conductive and elastic materials to reduce heat transfer resistance and to facilitate droplet ejection. Furthermore, the membrane could be coated with a thin layer of mate rial to promote dropwise condensation. The membrane is cooled to provide necessary subcooling. Once th e droplets grow to desired sizes, the membrane is actuated to eject the dr oplets. A synchronized valve can be used to expand the working fluid. The droplet would then impinge on the surface of the evaporator (such as an electronic device) to provide cooling by evaporation. The vapor from droplet evaporation would then be colle cted and sent back to the compressor to complete the cycle. The ejection of droplets, which is virtually gravity independent, serves the dual purposes of droplet remova l for maintaining dropw ise condensation and creating a spray for highly efficient cooling proc ess. There are several advantages of such a system: (1) higher reliability because there are no mechanical moving parts, (2) closed re-circulation which is self-contained with no need for extra maintenance or user intervention, (3) higher scalab ility because a system of va rying sizes and performance level can be made in a modular manner, (4 ) lower noise and vibration than standard condenser, (5) higher rate of energy transfer per unit weight and size (higher efficiency), (6) reduced sensitivity to the gravitational effect.
14 HQ Heat rejection Heat absorption LQ Valves Condenser (2-3) Evaporator (4-1) Expansion (3-4) Membrane-EjectorinW Droplets Cold side H ot side Compressor (1-2) (a) (b) Figure 2.3: Schematics of the conceptual micro-scale cooling device. (a) Simple schematic and (b) cross s ectional view schematic. For the sample calculation the cold side is at T1=20C and the hot side is at T3=50C, i.e., T=(T3-T1)=30C. We consider four different refrigerants: water, R-12, R134a and FC-72. R-12 is the dich lorodifluoromethane CCl2F2 of the CFC family commonly named Freon. It is fully halogenated and causes dama ge to the ozone layer. R-12 is replaced by R-134a, a chlorine free refrigerant, the tetrafluoroethane CF3CH2F. FC-72 is a synthetic fluorinert. It is ozone safe, non-toxic and non-flammable. The most useful property for
15 FC-72 is that it is highly dielectric and therefore an excellent electrically insulating fluid for interfacing with electroni c devices. Table 2.1 summarizes the important variables of the cycle based on T=(T3-T1)=30C for the four different fluids . Note that for FC-72, at the exit of the compressor the working fluid is at the state of saturated two-phase with a temperature at T2=T3=50C because of the specia l properties of FC-72. Calculations have been performed with three different T=(T3-T1)=20C, 30C and 40C where T3 is kept at 50oC. The effect of T=(T3-T1) on the COP and the compression ratio for the four refrigerants are shown in Figure 2.4 (a) a nd (b) respectively. The compression ratios for water are larger than those for R-134a because of their different properties, but the COP values are found in the range of 8.0 to 8.5 for all four fluids. Table 2.1: Ideal cycle characteristics for T3=50C, T1=20C. Fluid P1 (kPa) P2 (kPa) P2/P1 COP T2 (C) Water 2.34 12.32 5.3 8.5 165.5 R-12 567 1219 2.2 8.2 53.4 R134a 572.3 1318 2.3 8.0 52.8 FC-72 23.5 79.5 3.4 8.1 50.0 The actual compression process involves irreversibilities such as friction, which increase the entropy production and heat loss. In the computation we have considered a compressor efficiency of 80% and the results are presented in Table 2.2. The COP values are dropped to between 6.5 and 6.8. Other th an the non-isentropic compression, the actual vapor-compression cycle differs from the ideal one mostly due to the irreversibilities occuring in the components. Two common sources of irreversibilitie s are pressure drops and heat transfer across a temperature differen ce. However in our system irreversibilities due to pressure drop are not dominant since we have liquid droplet s. We only consider
16 irreversibilities du e to heat tran sfer through a temperat ure difference of 3C. The COP has been evaluated for ideal and actual vapor-compression cycles. Results for the different cycles are shown in Table 2.3. COP d ecreases as a result of irreversibility due to heat transfer, and due to the compressor effi ciency of 80%, the CO P is about 5 in the worst case for T=(T3-T1)=30C. Table 2.2: Cycle with a compressor efficiency of 0.8. Fluid P1 (kPa) P2 (kPa) P2/P1 COP T2 (C) Water 2.335 12.32 5.3 6.8 201 R12 567 1219 2.2 6.6 57.5 R134a 572.3 1318 2.3 6.4 56.3 FC-72 23.5 79.5 3.4 6.5 50.0 Table 2.3: Effect of irreve rsibility on the coefficient of performance for R134a. Case Fluid R134a COP 0 Ideal Carnot refrigeration cycle 10.8 1 Ideal Rankine refrigeration cycle (isentropic expansion s3=s4) 8.2 2 Ideal vapor-compression refrig eration cycle (Figure 2.2) h3=h4 8.0 3 Case 2 with compressor efficiency =0.8 6.4 4 Case 2 with irreversibility due to heat transfer (through a temperature difference of 3C) 6.4 5 Combination of Case 3 and Case 4 5.1
17 COP=QL/Win4 5 6 7 8 9 10 11 12 13 14 15 101520253035404550(T3-T1) Water R-12 R-134a FC-72 Reference point (a) Coefficient of Performance Compression ratio P2/P11 2 3 4 5 6 7 8 9 10 11 101520253035404550(T3-T1) Water R-12 R-134a FC-72 Reference point (b) Compression ratio Figure 2.4: Coefficient of performance (COP) and compression ratio (P2/P1) versus T=(T3-T1) with T3 =constant=50C for ideal cycle. (a) COP, (b) P2/P1.
18 In the following, an estimation of two key el ements related to the main components of the present device is made. Specifically, th e condensation and evaporation time scales and the power required for ejection are evaluated. 2.2 Phase Change Analysis The condensation and evaporation time scal es are evaluated based on the same basic approach. Therefore as an illustration we pr esent herein only the condensation time scale estimation. The condensation time scale is based on the following assumptions: â€¢ The heat flux is uniform and constant between the droplet and the membrane surface â€¢ The droplet is in the shape of a segmente d sphere during growth and a hemisphere at the moment of ejection (see Figure 2.5) â€¢ The droplet embryo is formed almost in stantaneously after th e departure of its predecessor with the critical radius, and the is othermal growth is much faster than the isobaric growth. Therefore the following growth times are ba sed on the isobaric growth. R h Figure 2.5: 3-D schematic perspective of a liquid droplet on a solid surface. The total heat transfer Q during the lifetime of a sessile droplet is its mass m at ejection times the latent heat , since the sensible heat is assumed negligible: Qm (2.4) with 314 23hmR (2.5)
19 where is the density, and Rh is the radius of the hemispherical droplet at the moment of ejection. Since the area in contact between the droplet and the condenser surface during the growth is changing with time, we use half of the base area of the hemisphere in contact with the membrane at the moment of ejection as the average contact area for heat transfer analysis between the sessile droplet and the membrane. The average heat transfer area A in contact is 22h R A (2.6) The condensation time scale tc is h c R QQ t QAqq (2.7) where Q is the average heat transf er rate between the sessile droplet and the membrane and q is the corresponding average h eat flux. In the calculation and are taken at T=50C for the estimate of the conde nsation time scale and at T=20C for the estimate of the evaporation time scale. Based on the resu lts of Sadhal and Plesset (1979) on dropwise condensation, we consider values of heat flux ranging from 104 to 106 W/m2 as the droplet is growing. The results of the conde nsation and evaporation time scales presented in Figure 2.6 are for different heat fluxes an d droplet radii. The time scales of water droplet growth or evaporati on are ten times larger than those of R-12, R-134a, and FC72, because the latent heat for water is ten times larger than those of R-12, R-134a and FC-72. Again the estimates are produced fr om conservative assumptions and lower bond values such as the droplet si ze after the isothermal growth were neglected. The above time estimates are based on droplet growth fr om zero size in the isobaric period. Most likely, the heat transfer ra tes would be higher than 105 W/m2 on a special membrane. For
20 a heat flux of 105 W/m2 and for droplet radi us greater than 300 m the time scale is about 1s. So with an ejection frequency of 1 Hertz the heat transfer ra te can be at least 0.01 W. If the heat flux is ten times hi gher, the frequency will be ten times higher, as well as the heat transfer rate. For a heat flux of 106 W/m2 and for droplet radius greater than 300 m the time scale is about 0.1 s, i.e., a frequency of 10 Hertz, and the heat transfer rate can be higher than 0.1 W. Depending on a systemâ€™s re quired heat transfer rate, the number of droplets to be condensed simultaneously, a nd therefore the device dimension, can be determined. For example, for a droplet with a radius of 300 m and a heat flux of 106 W/m2, a 3 W system will require 30 droplets, and a circular surface area of radius 1.7 mm. In conclusion, we can say that for la rger droplets, the system requires lower frequencies of cycle operation.
21 Evaporation/condensation time scale versus droplet radius for water 0.01 0.1 1 10 100 1000 0 100 200 300 400 500 droplet radius [ m] q =104 W/m2 q =105 W/m2 q =106 W/m2 Evaporation / Condensation time scale versus droplet radius for R-12 0.001 0.01 0.1 1 10 100 0 100 200 300 400 500 droplet radius[ m] q =104 W/m2q =105 W/m2q =106 W/m2 (a) water (b) R-12 Figure 2.6: Evaporation/condensation time scale versus droplet radius, for (a) water, and (b) R-12.
22 2.3 Analysis of the Ejection Mechanism The power required by the membrane to eject the droplet is now evaluated by considering a simple energy balance. The en ergy requirement to eject the droplet is approximated as the sum of th e kinetic energy increase of the droplet, the potential energy increase of the droplet, the energy di ssipation due to the drag, and the contact surface energy due to the surface tension force. The droplet is assumed to be spherical after it has departed from the membrane. The drag is estimated assuming Stokes flow. 221 6 2indhERUHRmUmgH (2.8) where Ein is the energy input required to eject the se ssile droplet, U the droplet velocity at impact on the evaporator surface, H the height between the condenser surface and the evaporator surface, Rh the radius of the hemisphere droplet at ejection, Rd the radius of the droplet after ejection, 32dhRR , the viscosity of the fluid in the saturated vapor phase. A schematic of the system is given in Figure 2.7 along with th e three stages of the droplet motion. In stage 1, the droplet after it has grown to the desire size is ejected, then in stage 2, the droplet follows an upward motion and finally, in stage 3, the droplet impinges on a flat surface. An estimation of the energy required by the membrane to eject R-12 droplet formed on its surface for a given impact veloc ity U and distance be tween condenser and evaporator H is presented in Figure 2.8. The energy input re quired to eject is dependent essentially on the radius of the droplet and its impact velocity. For a smaller droplet radius the contact energy is dominant whereas for a larger droplet the kinetic energy is dominant. The effect of gravity on the energy required for ejection is small, as shown in Figure 2.9. The system is therefore gravity insensitive. The power requirement for
23 ejection is the frequency times th e energy input required for ejection Ein. Larger droplets require more energy for ejection. Stage 1 Stage 2 Stage 3 H U Fd=6 RU Wd=mg Figure 2.7: Schematic of the three stages in the droplet motion. Energy requirement for ejection versus droplet radius1.0E-13 1.0E-12 1.0E-11 1.0E-10 1.0E-09 1.0E-08 1.0E-07 1.0E-06 0100200300400500 dro p let radius [ m ] Kinetic energy [J] Potential energy [J] Dissipation due to drag [J] Contact energy [J] Total energy for ejection [J] Figure 2.8: Total energy requirement vers us droplet size for R-12 with U=1m/s and H=1cm. For smaller droplet radius the cont act energy is dominant whereas for larger droplet the kinetic energy is dominant.
24 Gravity effect on energy for ejection versus droplet radius1.0E-11 1.0E-10 1.0E-09 1.0E-08 1.0E-07 1.0E-06 0100200300400500droplet radius [ m] against gravity with no gravity along with gravity Figure 2.9: Effect of gravity on energy require d for ejection versus droplet radius for R12 with U=1m/s and H=1cm. The above thermodynamic and simplified com ponent analysis demonstrate that the proposed micro-scale thermal management system is a feasible concept and offers very attractive features. However, there are majo r engineering issues that need to be addressed. In particular, fabr ication of micro-scale valves and actuators suitable for our needs is a major challenge in micro-fluidics. It is expected that with rapid development of the MEMS technologies, the proposed concept can be manufactured and employed in the foreseeable future. In the meantime, developm ent of numerical simulation tool will help improving the design and refining the thermo-fluid analysis of this new thermal management device. Specifically, in this wo rk, a numerical techni que is developed to study drop impingement with heat transfer. Details will be presented in chapter 6. In the next chapter, a literature review on each of the drop mechanisms involved in the microscale device is given.
25 CHAPTER 3 3 LITERATURE REVIEW In this chapter, the fundamental mechanisms and a literature review on some of the existing analytical and numerical studies on droplet ejection and dr oplet impingement on a flat surface and the condensation and evapora tion of a single droplet on a flat surface, are presented. 3.1 Droplet Ejection 3.1.1 Droplet Ejector Devices and Physical Observations In the literature, there have been studies on micro-fluidic dr op ejector devices, especially in areas related to inkjet printing devices and in micro-pipetting of biological liquids. Those devices employ di fferent techniques, including: (i) piezoelectric device connected to a re servoir, ejecting dr ops through a nozzle (Fromm, 1984; Shield et al., 1987; Szita et al., 2001; Percin et al., 1997; Grassia, 1999). In Fromm (1984) and Shield et al. (1987) for inkjet te chnology and Szita et al. (2001) for aspiration and dispensing mi cropipette, the piezoelectric device is mounted around a glass tube that has a small diameter nozzle at one of its extremity. An application of a voltage pulse to the piezoelectric results in a sudden pressure pulse that creates a contraction of the tube wa lls. In Grassia (1999) the piezoelectric device drives a r ubber valve. In Percin et al. (1997), one wall of a fluid reservoir is made of a piezoelectric-actuate d flexible membrane having an orifice at its center. (ii) heater-generated bubble serving as a piston to expel drops (Tseng et al., 1998; Rembe et al., 2000; Chen et al., 1997; Peeters and Verdonckt-Vandebrooeck, 1997). (iii) acoustic waves created by piezoelectri c force (Badie and Frits de Lange, 1997; Badea et al., 1998; Zhu and Kim, 1998; Dijksman, 1999). The piezoelectric force creates disturbance in a liquid cavity that breeds capillary instabilities at the free surface.
26 (a) (b) (c) Glass tube Driving voltage Piezoelectric Fluid reservoir Nozzle Figure 3.1: Illustration of diff erent droplet ejector devices: (a) piezoelectrically actuated squeezing tube from Shield, Bogy and Talke (1987), (b) thermal ink jet from Tseng et al. (1998), and (c) acoustic waves from Zhu a nd Kim (1998). Reprinted with permission. An illustration of the above three types of dr oplet ejector device is given in Figure 3.1. Photographs of the experiments made with the device of Figure 3.1 (a) from Shield et al. (1987) are shown in Figure 3.2 for ethylene gl ycol. Typical pulse duration for this device lasts 10 to 20 s, with pulse amplitude from 50 to 90 V. A quick evaluation with a reference velocity taken to be 1 m/s for the cas es presented in Figure 3.2 indicates that Re is about 17 and We about 0.85. The photographs of Figure 3.2 are take n for four different pulse durations and amplitudes: (a) 10 s and 90 V, (b) 15 s and 70 V, (c) 20 s and 70 V, and (d) 20 s and 90 V. 10 s, the shortest pulse durati on, corresponds to a resonance condition. In case (a), one observes from Fi gure 3.2 (a) the formation and ejection of a single drop. The photographs of Figur e 3.2 (b) are for a pulse of 15 s, which does not
27 correspond to a resonance condition; accordi ngly, two smaller drops are formed. In Figure 3.2 (c), the pulse duration is 20 s, twice the resonance condition and drop ejection occurs in a way similar to that s hown Figure 3.2 (a) for a lower amplitude pulse of 70 V. In Figure 3.2 (d), with the same pulse duration and a higher amplitude compared to Figure 3.2 (c), a satellite drop is formed. (a) (b) (c) (d) Figure 3.2: Ejection of ethylene glycol droplets produced by a piezoelectric inkjet type device from Shield, Bogy and Talke (1987) , with pulse of (a) 90.0 V for 10 s, (b) 70.0 V for 15 s, (c) 70.0 V for 20 s, and (d) 90.0V for 20 s. Reprinted with permission.
28 In the micro-scale device, the use of a membra ne to eject droplets is proposed as seen in chapter 2. Flexible membrane s or diaphragms have also be en utilized in micro-fluidic devices as valves, pumps, and compressors. Membranes can be actuated by electrostatic, piezo-electric or thermal mechanisms to pressu rize a fluid in a cavity. Saif et al. (1999) has proposed an analytical model for electr ostatic actuation membranes. Jerman (1994) has presented a thermal actuation of bimetallic diaphragm valve system, and Guo and Lee (1999) have reported issues on silicon membra ne fabrication and testing. Smith et al. (1998) and Vukasinovic et al. (2000) have investigated a vibration induced droplet atomization system for thermal management applications (Vukasinovic et al., 2001; Heffington et al., 2002). In their system a piezoelectric transdu cer is attached to a thin metallic diaphragm which vibrates. A liquid dr oplet of volume about 0.1 ml is placed on this diaphragm. When the diaphragm vibr ates around 1000 Hz, the vibration induces capillary waves on the free surface of the dropl et which grow in amplitude until the entire droplet bursts into a small cloud of atomized droplets. A volume of 0.1 ml corresponds to a sessile droplet of radius 5 mm. Sequential photographs of their experiment are shown in Figure 3.3.
29 t/T=0 5.4 10.8 17.3 324.0 14.0 11.9 20.5 32.4 Figure 3.3: Rapid atomization of a 0.1 ml water droplet fo rced by stepped actuation at 1080Hz, from Vukasinovic, Glezer and Smith ( 2000). Image courtesy of B. Vukasinovic. Reprinted with permission. 3.1.2 Simulation of Droplet Ejection Numerical simulations of a single drop form ation and ejection have been performed by Fromm (1984), Adams and Roy (1986), Shield et al. (1986, 1987), Badie and Frits de Lange (1998) for the drop-on-demand mechanism in ink-jet printing, and by Wilkes et al. (1999), Wilkes and Basaran (2001) and Gueyffi er et al. (1999) for a pendant drop. Fromm (1984) has developed a combined Eu lerian-Lagrangian method to numerically solve the axisymmetric Navier-Stokes equations for incompressible, free-surface flows. The streamfunction-vorticity formulation is adopted for the computations. The freesurface is represented by a Lagr angian net of points that move with the fluid. Given a pressure and velocity history at the nozzle, the time evolution of a drop as it grows, detaches from the nozzle, and deforms until a spherical shape of constant velocity is reached are calculated. Fromm (1984) empl oys as the nondimensional parameter R/W which is in fact the inverse of the Ohnesorge number:
30 1 2 21 Rr WOh (3.1) Numerical simulations are pe rformed for different values of Oh: 1/2, 1/3,1/4 for constant peak pressure of 80 units, and differe nt peak pressure pulse units: 60, 80, 100 for constant droplet Oh of 1/3. Figure 3.4 (a) shows the results for Oh=1/3, and Figure 3.4 (b) for Oh=1/4. The velocity of the leading edge droplet is investigat ed. It is found that a droplet with a smaller Oh number travels faster than that with a larger Oh (a droplet with larger Oh being more viscous). Detachment or break-off happens earlier for a droplet with a smaller Oh. After detachment a slight increase in the droplet leading edge velocity is observed. The higher velociti es distribution is identified in the droplet tail. When a higher peak pressure pulse is ap plied, a satellite drop emerges at later times. The satellite drop has higher average velocity than the main drop in all cases; furthermore, it is found that the main drop and satell ite drop undergo oscillations.
31 Figure 3.4: Instantaneous shapes, streamlines and pressure contours of droplet ejection from the numerical simulation results of Fromm (1984). (a) Oh=1/3, and (b) Oh=1/4 for a 80 unit peak pressure pulse. Streamlines are quasi parallel to the horizontal axis, and pressure contours are quasi perpendicula r to it. Reprinted with permission. Adams and Roy (1986) have proposed a mode l of drop detachment based on the onedimensional incompressible axisymmetric equa tions written in a Lagrangian framework. The system of first order differential eq uations is solved numerically by using MacCormackâ€™s algorithm. Their results agre e qualitatively with the computations of Fromm (1984). The one-dimensional model predic ts an earlier break-off and faster drop velocities compared to Frommâ€™s model. Shield et al. (1986, 1987) have also solved numerically the one-dimensional incompressible flow equations in an Eulerian formulation with M acCormackâ€™s algorithm.
32 They have considered two schemes. The firs t scheme, scheme A, takes into account the radial inertia, and the other one, scheme B, neglects it. The nozzle exit velocity history derived from the driving pressure history is employed as the driving boundary condition. Simulations are performed for a Reynolds num ber of 2.5, and different Weber numbers of 1, 5 and 20. The characteristic length is taken to be the nozzle exit radius. The calculation starts with a hemis pherical initial droplet shape. It is found that scheme A with inertia effect gives better results at higher We numbers and is able to model oscillations, whereas scheme B at low We gives a better deta chment profile. Their results are in accordance with those of Fromm (1984 ). While their numerical simulation can predict the overall phenomena , quantitative computation re quires a more comprehensive model than a one-dimensional a pproach. In a study of a dr op-on-demand inkjet system by Badie and Frits de Lange (1997) , a numerical simulation of a drop ejection is performed using a finite element method. The axis ymmetric Navier-Stokes equations for incompressible flow are solved. The veloci ty profile obtained from hydroacoustical calculations is prescribed. At the free su rface the dynamical and kinematical boundary conditions excluding surface tension gradient s are applied. The special shape of a pendant drop is observed, followed by neckin g (or constriction). The velocity vectors close to the nozzle are oriented in the negative direction of the axial axis, and the velocity vectors in the drop leading edge are oriented in the positive di rection leading to break-off. Wilkes et al. (1999) have studied comput ationally and experimentally the dynamics of a drop formation from a capillary tube of radius 0.16 cm into an ambient gas, referred to as the pendant drop problem. They have solved the axisymmetric Navier-Stokes equations for an incompressible Newtonian fluid with the dynamical and kinematical
33 boundary conditions at the free-surface. The set of equations is solved using a Galerkin/finite element algorithm on a multi-region mesh that conforms to and evolves with the drop shape changes. The mesh is composed of quadrilateral elements obtained by tessellation of the domain. Their computationa l model is able to capture both the gross features of the phenomenon, such as the limiting drop length at breakup and the volume of the primary drop, and its fine features, such as a micro-thread that develops from the neck region of a viscous drop that approach es breakup. It is observed that the flow reverses in the middle of the neck region, exhibiting a higher pressure there. Their numerical simulations show for the first time that the interface of a viscous drop overturns before it ruptures. They have co mputed the limiting drop length and primary drop volume over a wide range of nondimensiona l parameters. It is found that as the viscous effects increase, the length of the thread increases significantly, and the volume of the drop decreases almost linearly. Wilkes and Basaran (2001) have studied computationally drop ejection from an oscillating rod. They have employed the same numerical approach as in Wilkes et al. (1999). Several numerical simulations are perf ormed to study the effect of the forcing amplitude and frequency on the drop deformation until break-up. It is found that with increasing amplitudes at a constant fre quency, the drop elongates more, causing the formation of a liquid thread, and often a satellit e droplet at a later st age. For most of the frequency, drop ejection occurs during the second oscillation period, but as the frequency continues to increase, drop ejection eventually is delayed to a later period. In their study, the characteristic length scale is the rod radius ( R), and the characteristic time scale is the capillary time 3 ctR . A sequence of drop formation until break-up along with the
34 forcing amplitude and resulting drop length is shown in Figure 3.5 for an initially hemispherical drop, Re20We , dimensionless forcing amplitude A=0.280 and frequency =3.5. Plots of instantaneous mesh de formation, pressure contours and velocity vectors are shown in Figure 3.6 (a), (b) and (c), respectively, for the same condition. The pressure contours shown in Fi gure 3.6 (b) vary primarily in the axial direction when the drop has a maximum elongati on, i.e., at t=1.022 and t=3.036, and vary in the radial direction when the drop is at a minimum elongation, i.e., at t=1.701. The velocity field of Figure 3.6 (c ) is plotted in a fixed refere nce frame. Velocity vectors follow the direction of the rod motion at t=1.02 2, and t=1.701. At t=3.036, as the rod has moved in the upward direction, there is a large pressure gradient along the upward direction, and a small gradient s in the downward dire ction near the narrowest point of the neck. The pressure gradients cause a strong fl ow out of the neck in the upward direction and a weak flow in the downward direction as s een in the velocity fi eld of Figure 3.6 (c). This causes the neck to furthe r shrinks in diameter and leads to the formation of a drop. At t=3.539, while the fluid is entrained by the downward motion of the rod, it can no longer enter the thinning neck . Accordingly, the drop is at a point of break-up.
35 Figure 3.5: Drop shape, length (L), and rod ve rtical displacement (zâ€™) versus time t for Re20We , from Wilkes and Basaran (2001). Initially the drop is a hemisphere. Figure courtesy of E. Wilkes. Reprinted with permission. Gueyffier et al. (1999) have also numerically simulated the pendant drop dynamics using a three-dimensionsal two-phase flow so lver. Their solver computes the full NavierStokes equations for the incompressible fl ow on a fixed grid using the VOF method (volume of fluid). The simulation is perfor med in a rectangular box. The density ratio between the two fluids is 133, Re=4.2, We=9.216 x 10-5, Bo=1. Their numerical simulation is compared with photographs of Peregrine et al. (1990). The agreement is generally good, except that the simulation ove restimates the volume of the main drop. They have also performed axisymmetric si mulations. The axisymmetric simulations are found in good agreement, and do not overestim ate the drop volume. A summary of the numerical simulation stud ies presented above is offered in Table 3.1.
36 (a) adaptive meshes (b) pressure contours (c) velocity field Figure 3.6: Pendant drop undergoing oscillations for Re20We , from Wilkes and Basaran (2001). (a) Grid deformation, (b) pr essure contours, and (c) velocity field. Initially the drop is a hemisphere. Figure courtesy of E. Wilkes. Reprinted with permission.
37 Table 3.1: Summary of the numerical simulations of droplet ejection. Problem References Formulation and Method Major Results Fromm (1984) axisymmetric Navier-Stokes equations incompressible, free surface flow vorticity-streamfunction formulation combined Eulerian-Lagrangian method Oh= 1/2, 1/3, 1/4 higher drop velocity for small Oh after break-off, slight increase in the velocity at the leading edge satellite drops at later time satellite drop velocity higher than main drop Adams and Roy (1986) 1D problem incompressible, free surface flow Lagrangian approach Mac-Cormackâ€™s predictorcorrector scheme Re/We=2.5, 5, predict break-off earlier than Frommâ€™s model for Re/We=5 good agreement until break-off, then beyond break-off 1D simulation indicates higher velocities viscous effects retard the liquid column Shield, Bogy and Talke (1986,1987) 1D problem incompressible, free surface flow Radial inertia terms included in scheme A and omitted in scheme B Eulerian approach Mac-Cormackâ€™s predictorcorrector scheme Re=2.5, We=1, 5, 20 nozzle radius= 50 m reference velocity 1 m/s reference time 50 s comparable results to Frommâ€™s model scheme A give better results at higher We numbers and model oscillations scheme B at low We gives cleaner break up profile Drop on demand Badie and Frits de Lange (1997) axisymmetric Navier-Stokes equations incompressible, free-surface flows Finite element method Re=188, We=36 nozzle diameter= 50 m velocity vectors in opposite directions during necking phenomenon Wilkes, Phillips, and Basaran (1999) axisymmetric Navier Stokes equations incompressible, free surface flow Galerkin/finite element multiregion mesh conforming to drop shape tube radius=0.16 cm 0.037 < Re < 6.6 3.8 x 10-4 < We < 5.2 x 10-4 0.34 < Bo < 0.48 as viscous effects increase, thread length increases, and drop volume decrease Wilkes and Basaran (2001) same as Wilkes et al. (1999) moving reference frame to simulate an oscillating rod larger amplitude causes longer drop elongation, and will cause drop satellite results for Re/20 We shown in Figure 3.5and Figure 3.6 Pendant drop Gueyffier, Li, Nadim, Scardovelli, and Zaleski (1999) 3D and axisymmetric NavierStokes equations Two-phase flow Eulerian fixed grid VOF method Re=4.2, We=9.216 x 10-5, Bo=1 density ratio = 133 axisymmetric model gives better agreement to Peregr ine et al. (1990) experiments than 3D solver
38 3.2 Droplet Impingement 3.2.1 Experimental Observations of Droplet Impingement The system of a drop laying on a flat surface schematically shown in Figure 3.7, is characterized by three interfaces: liquidvapor, liquid-solid, and solid-vapor. The intersection of all three phases is a closed line and is called the contact line. If the drop is axisymmetric, then the contact line is the ba sal circle of the droplet. When seen in two dimensions this contact line is a point, calle d the contact point or tri-junction point. The wettability of a liquid reflects its affinity for solids. The wettability of the liquid is quantified by the contact angle , defined as the angle betwee n the liquid-vapor interface and the solid surface. For value of between 0 and 90, the liquid is said wetting and for value of between 90 and 180, the liquid is said non-wetting. When the system is at equilibrium we have th e Young-Dupr relation : coslvsvsl (3.2) with the contact angle, the surface tension, subscripts s for solid, l for liquid and for v gas-vapor. This relation can be derived either by considering a force balance in the horizontal direction or by a thermodynamic analysis using the concept of minimum energy. The derivation can be found in Carey (1992). The contact angle in the YoungDupr relation is called the static or equilib rium contact angle. The equilibrium contact angle is often experimentally determined by visualization techni ques (Dussan V.,1979). With a moving contact line, the problem beco mes dynamic. The contact line is advancing when its speed U is positive, i.e. U>0, and re ceding when U<0. The extrapolated value of the contact angle in the limit as 0 U with U>0 is called the advancing contact angle A , and with U<0 is called the receding contact angle R (Dussan V.,1979). Experiments
39 have shown that there is a difference betw een the advancing and receding contact angle values. The advancing contact angle A is typically larger than the receding contact angle R. This difference is referred to as the cont act angle hysteresis. Th e contact line does not appear to move when lies between , R A (Dussan V.,1979). The contact angle hysteresis is acknowledged to be a conse quence of three fact ors: (1) surface inhomogeneity, (2) surface roughness, (3) im purities on the surface (Carey, 1992). Furthermore, the contact angle has been obs erved to increase with temperature (Chandra and Avedisian, 1991; Lee et al., 2001; Nagai and Carey, 2001). Since it is often difficult to quantify these factors, an unambiguous expression of the contact angle under a dynamic condition is difficult. vapor solid liquid Figure 3.7: Schematic of a sessile droplet on a flat surface. Many studies on droplets impinging on very hot surface (above 400C), resulting in film boiling have been conducted, among th ose are the studies of Wachters and Westerling (1966), Bolle and M oureau (1976) and Bernardin et al. (1997). These studies are reviewed in Aydin and Yang (2000). When a liquid droplet free of any attach ment impacts onto a solid surface, it can splash, spread or rebound depending on the im pact velocity, droplet diameter, liquid and
40 solid surface properties. The behavior of a droplet impinging on a surface is quite complicate to understand and model because it involves imbalances between the inertia, surface tension and viscous forces. Zapalowicz (1995) in his qualitative analys is of the process of wetting a surface by a water droplet describes the effect of the Weber number as follow: (i) For small Weber number the droplet spread s until the liquid reaches equilibrium. (ii) At intermediate Weber numbers, the droplet after reaching a momentarily equilibrium, shrinks and diminishes its contact area until a next equilibrium is reached. (iii) For high Weber number below the critical Weber number, the droplet retracts and rebounds from the surface. (iv) When the Weber number is larger than th e critical value, the droplet breaks into smaller droplets upon impact, al so referred as splashing. The Weber number represents mainly the e ffect of the droplet size and impact velocity. The critical Weber number or sp lashing threshold not only depends on the droplet properties but also on the surface pr operties. Surface roughness is an important factor, and several splashing threshold correlations are reviewed in Bussmann (2000). The second important factor is the surface temperature as qualitatively discussed in Zapalowicz (1995). Chandra and Avedesian (1991) have studied the effect of surface temperature on the impact dynamics. They have experimentally investigated the defo rmation of n-heptane (C7H16) droplets for We=43 and Re=2300 impinging on a polished stainless steel surface held a constant temperature for temp erature ranging from 24C to 250C. The experiments are carried out at a constant ambient pressure. Th e droplet of initial diameter
41 of 1.50 mm is impacted at 0.93 m/s. The boiling point of n-heptane is 98.4C. Selected photographs of Chandra and Avedisian (1991) are shown in Figure 3.8. Once the droplet is in contact with the surface, it spreads radially and forms a flat disc shape. The spreading rate is driven by the liquid inertia, and is slowed by the viscous effects. The surface tension causes the drop to retract to an equilibrium shape. The drop can also oscillate between equilibr ium shapes. The natural frequency fN of the recoiling process of a drop on a solid surface is (Fukai et al., 1993; Chandra and Avedisian, 1991; Kaviany, 1994): 1/2 2316Nf d (3.3) The kinetic energy of the droplet before im pact is either conve rted into potential energy (due to surface tension) or dissipated. Th e lateral flow reversal at the late stage of the deformation process underlines the signif icance of the surface tension with a low We. A simple energy-based model has been deri ved by Chandra and Avedisian (1991) to obtain an expression for the maximum spread ing coefficient. The kinetic energy of the droplet before impact is given by 121 2k E mU (3.4) where m the mass of the droplet of density and diameter d is 31 6 md (3.5) The droplet surface energy before impact is 12 s E d (3.6)
42 The droplet after impact is assumed to have a disk shape of thickness and diameter dbase such that its surface energy is 221 1cos 4sbaseEd (3.7) The dissipation energy 2d E due to the viscous effect is approximated as: 2dvissp E t V (3.8) where vis is the dissipation function repr esented by the viscous effect: 2 visU (3.9) V is the liquid dropl et volume, and tsp is an estimation of the time the droplet of diameter d takes to spread at speed U: spd t U (3.10) Thus, the dissipation energy is expressed as: 22 36dUdd E U (3.11) The total energy is conserved, therefore 1122kssd E EEE (3.12) which leads after some algebr a to the following expression: 423 1cos40 23 We Ca (3.13) where is the spreading coefficient defined as:
43 based d (3.14) Figure 3.8: The impact of a n-heptane droplet on a heated stainless steel surface (view at an angle 30) for Re=2300 and We=43, from Chandra and Avedisian (1991). (a) Tw=24C, (b) Tw=90C, (c) Tw=104C. Image courtesy of S. Chandra. Reprinted with permission from The Royal Society. Recently, Kim and Chun (2001) have studied experimentally and analytically the recoiling of water, ink and silicone oil dropl et of about 3 mm diameter on polycarbonate
44 and silicon oxide solid surfaces which corresp onds to the following ranges of Weber and Ohnesorge numbers and contact angle values : 30 < We < 582, 0.0017 < Oh < 0.109, and 6.2 < < 87.4. Lee et al. (2001) have characterized the heat tr ansfer of an impinging PF-5060 single droplet using a micro-scale h eater array at low superheat (below nucleation) and intermediate superheat (where nuc leate boiling occurs ). The Reynolds number in their experiment is about 650, and Weber number about 11. They have investigated the evolution of the contact angl e, droplet diameter, a nd heat transfer with time. Their results show two parts in the va porization process: firs t with an effective transient heat transfer coefficient, and second with a constant heat transfer coefficient. They have suggested the initia l transient to be affected by a combination of the heat conduction into the liquid, the oscillatory mo tion of the droplet, and external diffusive vapor region. 3.2.2 Simulations of Droplet Impingement In the literature, the fluid dynamics problem of a droplet impacting on a flat surface has often been treated as an axisymmetric free surface problem. Simulations have been performed for initially spherical drop from th e time the drop is in contact with the surface until the liquid comes to rest. Different num erical approaches have been employed to solve this problem. The Eulerian formulati on with volume tracking method is employed in Harlow and Shannon (1967), Hatta et al . (1995), Bussmann et al. (1999, 2000) and Pasandideh-Fard et al. (2001). The Lagrangian formulation, where the interface is tracked directly by using adaptive grid that follows the interface deformation, has been adopted by Fukai et al. (1993, 1995, 1997) and Zhao et al. (1996a). Other techniques combining
45 Eulerian and Lagrangian formulations have also been applied by Zheng and Zhang (2000). Harlow and Shannon (1967) are the first to have solved numerically the fluid dynamics of a liquid drop splashing on a flat plate. They adopt the axisymmetric NavierStokes equations using the Marker-and-Cell (MAC) algorithm (Harlow and Welch, 1965; Welch et al., 1966). Since they neglect surf ace tension and viscous forces, the drop spreads with no constraint. Hatta et al. ( 1995) have followed a m odified MAC method to obtain solution of the Navier-Stokes equati ons including the viscous and surface tension forces. The results of their numerical si mulation are in good agreement with their experiment results for a water droplet and w ith those of Chandra and Avedisian (1991) for an n-heptane droplet at 24C at an earlier stage. For the water droplet case (We=62.5, Re=1500), after impingement, the droplet spr eads on the solid surfa ce in a thin disc shape, forming a doughnut shape before reco iling. For the simulation with n-heptane droplet case (We=43, Re=2300), the droplet sp reads until reaching a maximum diameter. However, no recoiling process is observed. Bussmann et al. (1999) have developed a three-dimensiona l volume tracking model to study the impact of a droplet onto an inclined surface and onto a sharp edge. Their numerical results are in good agreement with their experiment. Bussmann et al. (2000) have investigated the effect of fingering and splashing of a dr oplet impinging a solid surface by imposing a perturbation function. Pasandideh-Fard et al. (2001) have experimentally investig ated and modeled the fluid dynamics and heat transfer phenomena occurr ing in the substrate and droplet during droplet impingement on stainles s steel at different temperature varying from 50C to
46 120C (below nucleation) for different impact velocity, which corresponds to Re between 1200 to 9000 and We between 7 and 450. They have made use of an Eulerian rectangular grid. Two energy equations, one for the liqui d and the other one for the substrate are solved. Their numerical results agree well w ith their experiment. They have concluded that the effect of an increasing impact velo city makes the droplet spread more, therefore increasing the wetted area across which heat transfer takes place. Fukai et al. (1993) have developed a fi nite element technique to solve the axisymmetric Navier-Stokes equations based on a Lagrangian formulation of a droplet impinging on a flat surface. They accurately track the droplet deformations using deforming triangular-elements. The artificial comp ressibility method is implemented to compute the flow solution. Their model accoun ts for the presence of viscous, gravitation and surface tension forces. The normal a nd tangential stress balance conditions are satisfied on the free surface of the droplet. At the contact point or tri-junction point a net interfacial force is imposed and the axial velocity is se t to zero. Their numerical simulation is able to predict the formation of a propagating ring structure due to the liquid mass accumulation at the periphery of the thin disk, as well as recoiling and subsequent oscillation. The effect of impact velocity , droplet diameter, surface tension, and liquid properties on the fluid dynamics of the splashin g of the droplet is studied. In Fukai et al. (1995), coordinated experimental and numeri cal investigations are performed. The contact angle hysteresis, determined experime ntally, serves as input to the numerical model. Their numerical results are in acco rdance with their experimental results. In Fukai et al. (1995), three numerical si mulations are performed using the same technique, for an n-heptane droplet of diamet er 1.5 mm with impacting velocity 0.93 m/s,
47 corresponding to the experimental cases of Chandra and Avedisian (1991). Three surface temperatures are considered: 24C, 100C, and 150C. The surface temperature effect is accounted for in the simulation by assigning different contact angle values. Zhao et al. (1996a) present a coupled numerical study of the fluid dynamic and heat transfer phenomena during the impingement of a drop let on a flat surface. A model for heat transfer in the liquid and substrate is added to the fluid dynamics model of Fukai et al. (1993, 1995). Two energy equations, one for the li quid and the other for the substrate, are solved by enforcing temperature and heat flux continuity at the interface. Their numerical results of a water droplet impi nging on glass and aluminum are shown in Figure 3.9. It is found that the heat transfer occu rs mainly in radial direction, and that the heat transfer time scale is comparable to the droplet deformation time scale. Their numerical model is validated with experimental re sults in Zhao et al. (1996b). Zheng and Zhang (2000) have developed a new free surfa ce capturing method in three dimensions to simulate the spreading and so lidification of a droplet on a surface. The solidification interface is captured by the ad aptive grid, with the free surface interface tracked by the level-set method. At the contact line both a fixed contact angle condition and a mixed condition (fixed contact area and fixed angle) are implemented. Among their numerical results, a case of a water drople t impinging vertically on a flat surface for Re=100 and Fr=2.5 at an early stage is pres ented. Their results are in agreement with other two-dimensional numerical simulations . For example, the splat thickness is observed thinner at the center, and thicker at the periphery.
48 (a) Glass substrate (b) Aluminum substrate Figure 3.9: Numerical simulati on results from Zhao et al. (1996a) of the fluid dynamics and heat transfer phenomenon during the im pingement of a water droplet for Re=500, We=10 and Fr=630 on (a) a glass substrate and on (b) an alum inum substrate. Reprinted with permission from Elsevier Science. Finally, numerical studies simulating the im pact of a droplet on a thin liquid film recovering a flat surface have also been performed. Weiss and Yarin (1999) have formulated the problem as an axisymmetric potential flow problem since the viscous terms are neglected and solve it using boun dary integral method. Rieber and Frohn
49 (1999) have performed a full numerical si mulation of the Navier-Stokes in threedimensions. Table 3.2: Relevant numerical studies of a water liquid droplet impacting on a flat solid surface. The initial droplet diameter (prior of impact) is denoted by d0 and the impact speed by u0. References Model Simulation Main Results Fukai, Zhao, Poulikakos, Megaridis, Miyatake (1993) axisymmetric conservation equations free surface flow Lagrangian formulation finite-element method adapting triangular-element grid to drop deformation artificial compressibility method Case 1 d0 = 0.69 mm, U = 1.46 m/s Re=500 We= 10, Fr=630 Case 2 d0 = 0.2 mm U = 1 m/s Re=100, We=1.4 Fr=1020 We number is the most important factor influencing the behavior of the splat radius. results in good agreement with experiments at all times Zhao, Poulikakos, Fukai (1996a) same as Fukai et al. (1993) with heat transfer model energy equation in liquid and solid substrate mesh dimension based on thermal penetration depth d0 = 0.69 mm U = 1.46 m/s Re=500, We=10, Fr=630, Pr=7 Aluminum and glass surface sequence shown in Figure 3.9 heat transfer occurs mainly in radial direction simulations show that heat transfer time scale is comparable to droplet deformation time scale PasandidehFard, Aziz, Chandra, Mostaghimi (2001) axisymmetric Navier-Stokes for incompressible flow free surface flow energy equation in liquid and solid VOF method Eulerian fixed grid Experimental dynamic and receding contact angle values are input in the simulation d0=2.0 mm U =1.3m/s Re=2900, We=47 Stainless steel surface results in very good agreement with experiment greater impact velocity, greater droplet spread, greater wetted area where heat transfer takes place 3.3 Condensation and Evaporation of a Single Drop on a Flat Surface The condensation or evaporation of a single dr oplet on a flat surface, is referred in the literature as dropwise condens ation and dropwise evaporation (Sadhal et al., 1996). From a macroscopic point of view, condensation repr esents the change of phase from the vapor state to the liquid state. Evaporation is the reverse process of condensation, i.e. evaporation represents the change of phase from the liquid stat e to the vapor state. During
50 condensation the droplet gains mass, and durin g evaporation it depletes mass. From a microscopic point of view condensation ta kes place when vapor molecules impinge on the liquid surface and remain in the liquid. Evaporation takes place when the molecules escape from the liquid surface into the vapor, and are freed from the intermolecular force field acting in the liquid phase (Tanasawa, 1991). Condensation occurs when the number of molecules crossing a unit area of the vaporliquid interface per unit time into the liquid phase exceeds the number crossing toward the vapor phase. The reverse occurs during evaporation. The system is at equilibrium when the two fluxes are equal (Collier and Thome, 1994). Molecular dynamics studies an d computer simulations of condensation and evaporation have been performed by Fujik awa et al. (1995), Matsumoto et al. (1995), Carey et al. (1997), Matsum oto and Fujikawa (1997). In the heat transfer mechanism of a droplet on a flat surface that is cooled (or heated), the condensation (or evaporation) mechanism takes place at the liquid-vapor interface. The heat transfer in the solid results is non-uniform solid surface temperature and surface heat flux. High heat flux valu es are located next to the co ntact line, because of smaller liquid thickness. Microscopic models accounting for the physics of contact line have been developed and coupled with macroscopic models (Koba yashi et al., 1994; Wayner, 1993, 1994; Lay and Dhir, 1995; Son et al., 1999). The microsco pic model considers the existence in the vicinity of the contact line of very thin nonevaporative liquid film of the order of tens of nanometer in thickness. This film thickness is of the order of molecules size. Wayner (1994) has proposed a physical model for the contact line motion on solid surface due to condensation and evaporation of the liquid in the microscopic region based on kinetics
51 and interfacial concepts. Equati ons for the velocity, heat flux and superheat in the contact line region as a function of the change of the contact angle are derived. In the process of nucleate boiling of a single bubble, Son et al. (1999) have perf ormed a complete numerical simulation by coupling the solutio n of the microscopic and macroscopic region. A level-set method with phase change is employed in the macroscopic region. 3.3.1 Dropwise Condensation Heterogeneous nucleation occurs on a wall when the later is cooled as in the case of a surface condenser. Condensation on a soli d surface occurs in two different ways: dropwise and filmwise depending on the su rface wettability. Filmwise condensation occurs on a cooled surface easily wetted: th e liquid forms a film. On non-wetted surfaces the vapor condenses in drops which grow by further cond ensation and coalescence and then roll, fall, or detach from the surface unde r the action of gravity or other forces. New drops then form to take their place. A sche matic representation adopted from Tanasawa of the cycle of dropwise condens ation is illustrated in Figu re 3.10. The present review focuses on the mechanism of single drop let growth on a flat surface by direct condensation. Dropwise condensation offers th e attraction of increasing condensing heat transfer coefficients by an order of magnitude over filmwise condensation.
52 1)Formation of initial droplets 2) Growth due to direct condensation 3) Growth due to coalescence 4) Departure from surface Figure 3.10:Cycle of dropwise condensati on (adapted from Tanasawa, 1991). Two theories have been proposed to expl ain the mechanism of initial drop formation on a flat surface. One theory is based on the idea of â€œfilm ruptureâ€. This theory, first proposed by Jakob (1936) and later supported by Welch and Westwater (1961), states that vapor condenses on a dry subcooled surf ace and forms a thin film layer. When the film reaches a critical thic kness estimated to be about 1 m, it ruptures into droplets. The other theory first in troduced by Eucken (1937) and e xperimentally supported by Umur and Griffith (1965) and McCormick and We stwater (1965) is based on the idea of nucleation theory. The nucleation theory is be ing more accepted by researchers than the film rupture theory (Sadhal et al., 1996; Tanasawa, 1991). From a macroscopic point of view, however it does not matter, the way the droplet is originated. What matters is to know the critical radius rmin of the droplet for which cond ensation can be sustained. The smallest droplet, assuming hemispherical droplet ( = 90), size possible corresponds to the equilibrium radius of curvature for the specified wall subcooling s atwTT given in Carey (1992) is:
53 min2 lnwv vsatw lsatwRTp p pT vpTr (3.15) where is the surface tension of the liquid-vapo r interface, R is the gas constant for the vapor, lv is the specific volume of the liquid, wT is the wall temperature, v p is the farfield pressure of the vapor, and s atw p T, the equilibrium vapor pressure at the wall temperature. Here the foll owing approximation is made lnwv vsatw lsatwRTp ppT vpT (3.16) and Clausius-Clapeyron equation 11 lnv satwsatvwp p TRTpT (3.17) is used to yield min2satvl satvwTpv TpTr (3.18) The free energy of formation of a liquid em bryo on a flat surface is reduced by a shape factor () F that is function of the contact angle (Carey, 1992; Shyy, 1994; Collier and Thome, 1994): 22cos1cos 4() F (3.19) The growth of a single, surface-nucleate d droplet by condensation on the droplet surface is controlled by the heat transfer rate through the droplet and across the droplet solid interface (Kaviany, 1994). The transpor t process in a liquid droplet is often dominated by conduction, although there may be some convective transport in the liquid.
54 Since the conduction path is l onger through larger droplets than for smaller ones, larger droplets have greater resistance to heat tran sfer. The rate of condensation on the larger droplets is therefore less than on the smaller ones. Hence, the larger drops grow primarily by coalescence. This implies th at most of the heat transf er during dropwise condensation is transferred to the surface covere d with the smallest droplets. In the study of the heat transfer mechanism of the full cycle of dropwise condensation, two models have been proposed. The first one is the Le Fevre and Rose model (Le Fevre and Rose, 1966) and the second one is Tana kaâ€™s model (Tanaka, 1975). Le Fevre and Rose model derives the heat flux of dr opwise condensation by combining the heat transfer resistance of a single droplet of a given size and the average heat transfer coefficient on the surface to reflect the drop size distribution. Tana kaâ€™s model considers the transient change of local drop size distri bution, taking into acc ount the processes of growth and coalescence of drops. These models are reviewed by Rose et al. (1999) and in Tanasawa (1991) review paper on condensat ion. Computer simulation of dropwise condensation cycles have been done by Gose et al. (1967) and Burnside and Hadi (1999). Initially randomly distributed droplets are placed on a surface, then as time goes, the droplets growth, coalescence, and removal are simulated by using a simple droplet growth rate model and basic analysis. The expression obtained by Le Fevre and Rose (1966) for the heat tr ansfer resistance of a single droplet is derived by considering the following three aspects: (1) the influence of surface curvature on the phase equilibr ium temperature, (2) the mass transfer resistance at the liquid-vapor interface, a nd (3) the heat conduction resistance through the
55 droplet. The relation between the heat flux Bq from the base of a hemisphere drop to the condensing surface, and the surface subcooling T is expressed as: 2 1 221 12lsatvsatv v B lvTpRTp KvT Kr Tq rk (3.20) where s atvTp is the saturation te mperature of the vapor, r is the drop radius, is the surface tension, vl and vv are the specific volumes of th e liquid and vapor respectively, is the latent heat, kl is the liquid conductivity, is the ratio of the specific heat, K1 is a constant relating the temperature drop in the drop to Bq (K1=2/3 for hemispherical drop), and K2 is the ratio of the base area of a dr op to its curved surface area, i.e., for hemispherical drops K2=1/2. Other heat transfer resistan ce can be added to consider for example the resistance of a promoter layer if coated on the solid surface (Griffith, 1985; Abu-Orabi, 1998). Sadhal and Plesset (1979) have considered the heat transfer problem of a single sessile droplet of diameter less than 1 mm condensi ng on a semi-infinite flat surface. They have studied the effect of the so lid material. The droplet is assumed to have a spherical segment shape. They have solved the stead y heat conduction equati on for both the liquid and the solid and have considered vapor tran sport. The temperature distribution in both the liquid and solid phase is given by: 2 20 0l sT T (3.21) with the following boundary conditions: at the liquid-vapor interface
56 l llvT khTT n (3.22) at the solid-liquid interface ls ls lsTT TT kk nn (3.23) at the solid-vapor interface 0sT n (3.24) deep in the solid 0 sTT (3.25) where Tl is the liquid temperature, Ts is the solid temperature, and Tv is the vapor temperature, T0 is the temperature of the solid at large distances under the droplet, n the normal to the interface, kl liquid conductivity, ks solid conductivity, h the heat transfer coefficient. The equations are first transforme d in toroidal coordina tes, and then solved. Since the boundary condition given in Eq. (3.22) cannot be satisfied exactly, an approximate solution is sought. From this solution the Nusselt number Nu defined as 0 lvQ Nu krTT (3.26) is expressed as: 2 0 2sech 4 sintanhtanh tanhtanh, 2sinhl sd Nu k gBi kBi (3.27) where the correction function , g Bi is 2 14 3sin1sin ,1.8sin11ln 222 gBi BiBiBi (3.28)
57 where Bi is the Biot number defined as: lhr Bi k (3.29) and is the contact angle, r is the droplet base radius. Plots of the results of Eq. (3.27) (Nu versus =kl/ks) obtained by numerical integration are given in Sadhal and Plesset (1979) for =30, 60 and 90, 4101 , 41010Bi, and 1120Nu. It is observed that as th e liquid to soli d conductivity ratio approaches values greater than 1/Bi, Nu starts to decreases as the liquid to solid conductivity ratio increases. For values smaller than 1/Bi, the drop Nusselt number is constant for a given Biot number. Th e droplet growth rate is given by 3 0 2sin 1cos2coslv lkTTNu dr dtr (3.30) Because of the complexity of Nu function, no straightforward integration is possible. In Lai (1999), the growth of a single sphe rical droplet by condensation is studied theoretically. Initially a spherica l droplet having an initial radius larger than the critical radius is placed at the center of a spherica l domain filled with vapor. The vapor is initially at a uniform pressure and a subcooled temperature th at is less than the saturation temperature. The governing equations solved in the liquid and vapor phases are the unsteady, spherically symmetrical energy e quations. The energy equation in the liquid phase reduces to the heat conduction equati on since the convection terms are neglected. Three time scales are derived in the analysis. The first and second time scales reflect the thermal diffusion of the vapor and liquid. Th ese two time scales are small. The third time scale is the condensation time scale describi ng the droplet growth. The parametric range
58 investigates are: for the Jakob number (Ja) 10-7 < Ja < 10-5, for the Gibbs-Thomson effect number () 10-6 < < 10-3, and for the liquid to vapor density ratio 103 < l/g < 105. This parametric range is typical for the wate r vapor with droplet of radius ranging from 1 to 20 mm, subcooled temperat ure ranging from 1 to 10 C, and saturation pressure ranging from 0.1 to 1.0 atm. It is found that with a large Jakob number, that the subcooling effect is stronger and the droplet grows faster, that the Gibbs-Thomson effect suppresses the droplet growth and that induced convection due to large density difference between the liquid and vapor pha se contributes to the droplet growth. A summary of the analytical model of single droplet conde nsation is presented in Table 3.3. Table 3.3: Summary of the previous analytic al studies in single droplet condensation Authors Assumptions and Model Main Results Le Fevre and Rose (1966) hemispherical droplet on flat surface heat transfer resistance model taking account surface curvature effect, mass transfer resistance at interface, heat conduction in liquid phase relation between T subcooling and heat flux Eq. (3.20) Sadhal and Plesset (1979) spherical segmented sphere on semiinfinite flat surface, diameter < 1mm steady state heat conduction equation in liquid and solid phase =30, 60, 90; 4101l sk k , 10 < Bi < 104, 1< Nu < 120 approximate solution relation between Nu, , Bi and for values smaller than 1/Bi, Nu is constant for a given Bi. as approaches values greater than 1/Bi, Nu starts to decreases as increases. Lai (1999) spherical droplet surrounded by vapor droplet radius 1 to 20 mm heat conduction equation in liquid phase energy equation in vapor phase initially vapor is at uniform pressure and subcooled temperature 751010 Ja, 631010 351010l g condensation time scale is larger than the diffusion time scale droplet grows faster for stronger subcooling effect (large Ja) Gibbs-Thomson effect suppresses droplet growth large density difference between liquid and vapor phase contributes droplet growth
59 3.3.2 Dropwise Evaporation When a sessile droplet sits on a surface that is heated, different regimes of evaporation occur depending on the surface temperature, Tw (Zapalowicz, 1995; Sadhal et al., 1996; Aydin and Yang, 2000). In the study of evaporation, two curves, namely the droplet lifetime and the boiling curve are of fundament al interest. They ar e depicted in Figure 3.11, where the different evaporation regime s are indicated. These regimes are described in Table 3.4. In Figure 3.11, TONB stands for the temperat ure of the onset of bubble nucleation. Tmax is the temperature corresponding to the shortest droplet evaporation time and the critical heat flux, TLeid the temperature corresponding to the minimum heat flux. TLeid TONB Tmax tevap 3 1 2 4 Tw (a) 3 TLeid Tmax Critical Heat Flux Heat Flux Minimum Heat Flux 1 2 4 Tw (b) Figure 3.11: Schematic of the droplet evapor ation lifetime and boili ng curves. (a) Droplet evaporation time, tevap, versus surface temperature, Tw, and (b) Heat flux versus surface temperature.
60 Table 3.4: Classification and de scription of the dynamics and heat transfer of a droplet evaporation Zone Regime Superheat Description 1 Film evaporation Low When wONBTT , the droplet spreads and evaporates slowly Total contact of the droplet with the surface Conduction is the dominant mechanism inside the droplet 2 Nucleate boiling Intermediate When max ONBwTTT , vapor bubbles arise in the droplet Increase in momentum transport around the droplet causes increase of heat transfer Partial and intermittent contact of the droplets and solid surface 3 Transition boiling Intermediate When max wLeidTTT , the frequency of contact of the droplet and solid surface decreases, causing a decrease in heat transfer Formation of an unstable vapor film between droplet and solid surface 4 Film boiling High When wLeidTT , a thin and stable vapor film is present No contact between the liquid and solid surface The low superheat evaporation regime is reviewed in the following. In the low superheat evaporation regime, there is no nuc leation. The drop remains in contact with the surface during all the evaporation process. The heat transfer mechanism in the liquid phase is conduction dominated. The convect ion mechanism can be neglected. The smaller the drop is, the greater the cooling effect is. The droplet impact mechanism, reviewed in section 3.2, influences the h eat transfer and evaporation rates. When a droplet impacts on a flat surface, the droplet de forms to a shape that is more like a flat disk than a hemisphere. Such a shape allows wider conduction paths, and liquid convection parallel to the surface, improving the evaporation rate.
61 Sadhal and Plesset (1979) have obtained an exact analytical solution for an evaporating droplet. They have followed the same procedure that in the problem of a condensing droplet, except that the boundary condition at the liquid-vapor interface Eq. (3.22) is now for the case of an evaporating droplet: lvTT (3.31) For the evaporating droplet cas e, an exact solution is obtained. The Nusselt number defined in Eq. (3.26) is: 11 22 14tanss llkk Nu kk (3.32) where is the contact angle, kl is the liquid conductivity, and ks is the solid conductivity. The droplet depletion rate is given by Eq. (3 .30). Upon integration of Eq. (3.30), the life time or time for an initial droplet of radius r0 to be totally evaporated until its disappearance is: 2 2 0 31cos2cos 2sinl e lwsatr t kTTNu (3.33) where Tv is the vapor temperature, and T0 is the temperature of th e solid at large distances under the droplet. The Nusselt number Nu is plotted ve rsus the liquid to solid conductivity ratio for = 30 , 60 , 90 and 120 . For most liquid droplet evaporating from a solid metallic surface is of the order 10-2. For =10-2, the Nu is of the order 102 for each contact angle. The Nusselt number is larger for smaller contact angle for a given . In the dropwise evaporation studies of Bonacina et al . (1979), di Marzo and Evans (1989), Rizza (1981) and di Marz o et al. (1993), the liquid is assumed in full contact with
62 the solid surface, the only heat transfer mechanism considered is conduction, since convection is neglected. Bonacina et al. (1979) and di Marzo and Evans (1989) have developed simple analytical models, wher eas Rizza (1981) and di Marzo et al. (1993) solved numerically the conduction equation. In the model developed by Bonacina et al. (1979), a thin disk evaporating on a flat surface is considered. A linear temperature profile is assume d in the liquid. The following expression for the evaporation time of an initial disk of thickness 0 is obtained: 2 02llpsatl lwsatcTT t kTT (3.34) with l the density of the liquid, the latent heat, kl the thermal conductivity of the liquid, lpc the specific heat of the liquid, T the temperature. The subscrips l stands for liquid, sat for saturation, and w for wall or surface. They have extended their model to analyze the heat transfer in spray cooling by considering a coefficient that represents the fraction of area covered by the droplets. Thei r model is found to be in good agreement with their experiments. In the model of di Marzo and Evans (1989) , the solid surface is assumed isothermal, and is made of a highly c onductive material. The droplet shape is modeled by spherical segment. The temperature profile in the liqui d is assumed to be linear, and the droplet contact area with the solid remains constant during the evaporation process. Expressions of the instantaneous rate of evaporation and of the heat flux at the liquid-vapor interface are obtained by considering the energy bala nce for the droplet at the liquid-vapor interface. The Chilton and Colburn (1934) analogy at a liquid-air interface is used to express the heat flux q :
63 2 30.624 1aia conv pai x x D qh cx (3.35) and the rate of evaporation is deduced as: 2 2 1 3 020.624 1aconv ia lpairh xx dVD zdz dtcx (3.36) where 0.624 is the water-air molecular weight ratio, hconv is the convective heat transfer coefficient, r is the radius of the contact area, is the latent heat, cpa is the specific heat of the air, D is the diffusion coef ficient of the vapor in the air, a is the thermal diffusivity of air, i x is the vapor mole fraction at the liquid-air interface, a x is the air vapor mole fraction, and z is the nondimensional radial co ordinate. These expressions are found to be in reasonable good agreement with experimental data. It is observed that more heat transfer is taking place where the thickness of droplet is thinner. Rizza (1981) has numerically solved the heat conduction equation in two-dimensional cylindrical coordinates in the solid. The so lid surface temperature is estimated based on the local heat flux balance. The solid surf ace is made of a mate rial of high thermal conductivity and low heat capaci ty. The drop is assumed to ha ve a disk like shape. Heat flux lines are found to be highly concentrated at the outer radii. Ev aporation takes place from the droplet outer radius inward towa rd the droplet center and on the droplet upper surface. The conductance of the droplet is co ntinually increasing due to the reduction in the droplet diameter and thickness. Numeri cal solutions are performed for different droplet size and wall thickness. The nondimensional parameters b and ad e used in the investigation are defined as:
64 r bwettingratio R (3.37) where r is the impacted droplet radius, and R is the evaporation site radius. 2 r aspreadingratio d (3.38) 2 adr ee (3.39) where d is the droplet diameter before impact and e is the evaporator wall thickness. It is observed that a thicker wall gives shorter ev aporation time for a given impacted droplet radius and evaporation site radius, and that a larger evaporation radius site gives longer evaporation time for a given impacted radius droplet and wall thickness. In di Marzo et al. (1993), a coupled model is proposed where the transient conduction equation for the solid and the liquid are simu ltaneously numerically solved. The droplet shape is the one of segmented sphere of fixed base, i.e. th e contact area remains constant. The governing equations considered are: 2 2 l ll s s sT T t T T t (3.40) where l is the liquid thermal diffusivity, s is the solid thermal diffusivity, Tl is the liquid temperature, and Ts is the solid temperature with the following boundary conditions: at the liquid-air interface 2 3.0.624 1aia lllaconv pai x x D kTnhTTh cx (3.41) at the solid-liquid interface
65 ls ls lsTT TT kk nn (3.42) at the solid-air interface s s saT khTT n (3.43) deep in the solid 0 s sT kq n (3.44) where h is the overall heat transfer coefficient, Ta is the air temperature, 0.624 is the water-air molecular weight ratio, hconv is the convective heat transfer coefficient, r is the radius of the contact area, is the latent heat, cpa is the specific heat of the air, D is the diffusion coefficient of the vapor in the air, a is the thermal diffusivity of air, xi is the vapor mole fraction at the liquid-air interface, xa is the air vapor mole fraction, and q0 is the initial heat flux in the so lid. The transient c onduction equation in the solid is solved using boundary element method, and in the liqui d using finite volume element. Given the initial droplet volume, the droplet shape, th e initial surface temperat ure, the overall and convective heat transfer coefficients, thei r numerical code predicts the transient temperature profiles in the soli d and liquid, the heat fluxes, the total evaporation time, and transient droplet segmented sphere shape. Computations are performed for an initial droplet of 30 l on two different surface materials: al uminum and macor. The results of the simulations are found in good agreement w ith the experiments, especially at the droplet edge. A summary of the previous studi es on low superheat evaporation is given in Table 3.5.
66 Table 3.5: Previous studies on low superheat dropwise evaporation. Authors Study Assumptions and Model Main Results Sadhal and Plesset (1979) Analytical spherical segmented sphere on semi-infinite flat surface steady state heat conduction equation in liquid and solid phase =30 , 60 , 90 , 120 10-3< <1, 1< Nu < 500 exact solution Nusselt number as a function of contact angle and solid to liquid conductivity ratio Eq. (3.32) evaporation time Eq.(3.33) Bonacina, Del Giudice, and Comini (1979) Experimental and Analytical thin disk shape drop initial droplet diameter ~ 95 m disk diameter ~ 395 m linear temperature profile in droplet heat transfer only by conduction evaporation time of an initial thin disk Eq. (3.34) their model is in good agreement with their experiments Di Marzo and Evans (1989) Experimental and Analytical spherical segment drop isothermal solid surface linear temperature profile in drop constant contact area of drop and solid surface heat transfer only by conduction initial droplet volume=10, 30 50 l heat flux given by Eq. (3.35) rate of evaporation given by Eq. (3.36) high heat transfer where droplet thickness is thinner Rizza (1981) Numerical disk shape drop nonisothermal solid surface heat conduction in solid high thermal conductivity and low heat capacity solid material b=1/2, 1/4, 1/6, 1/8, 1/10 (ad/e)=1, 2, 3, 4, 5 heat flux lines highly concentrated under the contact line shorter evaporation time for thicker wall and smaller evaporation site radius for a given impacted droplet Di Marzo, Tartani, and Liao(1993) Experimental And Numerical spherical segment drop unsteady heat conduction equation in solid and liquid initial droplet volume=10, 20, 30 l aluminum and macor solid surface results in good agreement with experiment especially around droplet contact line
67 CHAPTER 4 4 THEORETICAL FORMULATION AND COMPUTATIONAL METHODS In this chapter, the theore tical formulation and the com putational methods adopted to solve multiphase thermo-fluid problems are presented. After describing the formulation, numerical assessment for single-phase flow computation will be made to establish the basis of the present approach. Then, in ch apter 5, multiphase computations will be assessed for rising bubble and drop dynamics. 4.1 Theoretical Formulation 4.1.1 Validity of Continuum Model First, the validity of a continuum model is discussed for modeling the thermo-fluid physical mechanisms involved in the microdevi ce. An important issue in modeling flows in micro-scale dimensions is the choice of the model to be employed to perform the simulation. The continuum model is no longe r valid when the Knudsen number becomes comparable to the characteristic length of the overall flow dime nsion. The Knudsen number is defined as: ...mfp Kn L (4.1) where m.f.p. denotes the mean free path length, and L a characteristic length scale of the flow. For 0.1 Kn , molecular models are required and for 3100.1 Kn, the continuum model can still be employed with a slip boundary condition at the wall (Lfdahl and Gad-El-Hak, 1999). The mean free path is defined as the average distance
68 that a molecule travels between successive collision and is for a single species gas (Vincenti and Krueger, 1975): 22 ..... 22molmolmM mfp dNd (4.2) where m is the mass of a molecule, dmol. is the diameter of a molecule, is the gas density, M is the molecular weight and N is the Avogadroâ€™s number N=6.02252 x 1023/mol. A simple calculation is now made for water to evaluate the Knudsen number in the vapor phase. The diameter of a molecule of steam is about dmol .=3.9 x 10-10m, the density of steam is 0.597 kg/m3, and the molecular weight of H2O is M=18.0156 gm/mol (Carey, 1992). Therefore, the mean free path of steam is m.f.p.=7.4 x 10-8m. In our case the characteristic length scale of the flow is taken to be the droplet diameter. Therefore, for a droplet diameter of 20 m, the Knudsen number is Kn=3.7 x 10-3m and for a droplet diameter of 1 mm, Kn=7.4 x 10-5m. Based on the above reasoning, it is clear that the continuum model is valid fo r droplet larger than 20 m in diameter. 4.1.2 Governing Equations The governing equations for a general multipha se thermo-fluid problems are the mass, momentum and energy conservation equations. For Newtonian and incompressible fluids with constant physical properties, the governing equations are: The mass conservation equation: .0u (4.3) where is the fluid density, u is the velocity vector. The momentum conservation equation with gravity: .. u uupu g t (4.4) where p is the pressure, is the dynamic viscosity, g is the gravitational acceleration vector.
69 The energy conservation equation: ..P PCT uCTkT t (4.5) where T is the temperature, Cp the specific heat and k the thermal conductivity. For the solid phase 0 u and 0 p ; there, the governing e quation represents the unsteady heat conduction process. The mass, momentum and energy governing equa tions in an axisymmetric coordinate system (r, , x) with /0 and no swirl velocity 0u are: 0 1 x ru x u r r r (4.6) 21 1 r u x u x r u r r r r p u u x u u r r r u tr r r r x r r r (4.7) x u x r u r r r x p u u x u u r r r u tx x x x x r x 1 1 (4.8) x T k x r T rk r r T C u x T C u r r r T C tP x P r P1 1 (4.9) 4.1.3 Interfacial Boundary Conditions The interfacial boundary conditions at the liquid-vapor interface with phase change (Shyy, 1994; Sadhal et al. 1996) are: The mass conservation condition is in the normal direction vvllmuVnuVn (4.10) in the tangential direction where the no-slip condition applies lvVtutut (4.11) where ( n) and ( t) are the unit normal and tang ential vector, respectively, V denotes the velocity of the interface, the subscripts v and l the vapor and liquid phase, respectively, and vu and lu the fluid velocity on the vapor and liquid side, respectively.
70 The conservation of momentum condition is in the normal direction 22vllvvv llVVnnnnppuu (4.12) or 211 lv vlv lnnnnppm (4.13) in the tangential direction .. .lvnnttt (4.14) where P is the pressure, is the surface tension, the curvature of the interface, is the viscous stress tensor. If the surface te nsion is considered constant, Eq. (4.14) simplifies to ..lvnntt (4.15) The conservation of energy conditions impose the continuity of temperature lvTT (4.16) the continuity of heat flux .llvv lluVnTTnmkk (4.17) where is the latent heat of conde nsation or evaporation, and k is the conductivity coefficient. For a curved liquidvapor interface, the su rface energy of the in terface causes a shift in the phase equilibrium. This effect is called the Gibbs-Thomson effect. The temperature of the interface is shif ted from the saturation temperature of a flat interface. For a heat transfer, the liquid pressure at the in terface is equal to the vapor pressure: lv s atppT (4.18) in this case, the Gibbs-Thomson equa tion in terms of temperature is: 1lv lsatTTT (4.19) where Tsat is the saturation temper ature on a flat surface, the surface tension, and the curvature of the interface.
71 When there is no phase change 0 m . Thus, Eq. (4.10), (4.13) and (4.17) simplify, respectively, to the following: vlVnunun (4.20) lvv lnnnnpp (4.21) vv llTnTnkk (4.22) For the solid-fluid interfaci al boundary conditions with no phase change, no slip and no penetration boundary conditions apply on th e solid surface. The energy conservation imposes the continuity of temperature and heat flux through the solid-liquid interface s lTT (4.23) ss llTnTnkk (4.24) and through the solid-vapor interface s vTT (4.25) ssvvTnTnkk (4.26) At the contact line, i.e., at the intersect ion between the liquidvapor-solid phases, a constant contact angle is considered. Therefore, the nor mal to the interf ace is defined as (Brackbill et al., 1992): cossinwallwallnnt (4.27) where is the contact angle, walln is the unit normal vector to the solid surface, wallt is the unit tangential vector to the wall.
72 4.1.4 Formulation with an Imm ersed Boundary Interface A single set of conservation equations valid for both phases that in clude the interface conditions can be derived. This formulati on considers the interface to be of small nonzero thickness within which the values of th e properties change smoothly. It is this formulation that we will numerically implement. In order to distinguish between the phases, a function , I xt named indicator function (Juric and Tryggvason, 1998; Brackbill et al., 1998) which ha s the value 1 in the liquid pha se and 0 in the gas-vapor phase, and that varies smoothly from 0 to 1 in the interface region is used. This function is similar to the color function used in the continuum surface force (CSF) method (Brackbill et al., 1992; Kothe et al., 1996). The fl uid property field at every location, such as density, viscosity, conductivity, specific heat, can then be given by: vlv x tIxt ,, (4.28) where denotes any fluid property, the subscripts l and v denote liquid and vapor-gas respectively. The indicator function , I xt has the properties of a heaviside function (Udaykumar et al., 1997; Sussman et al., 1994) . By definition, the gradient of the heaviside function , I xt is: k t I nxxds C (4.29) where k x x is the Dirac delta function that is zero except in the interface region, k x is the coordinate of marker points that represent the interface location, and t C is the interface contour.
73 The surface tension force of the interface is reformulat ed as a body force in the momentum equation. This model is referred as the continuum surface force model (CSF) (Brackbill et al., 1992). The momentum equation writes as: ..k tu uupunxxds t C (4.30) If phase change at the liquid-vapor interf ace is considered, additional source terms are added to the mass and energy conservation equa tions. This formulation has been used in the numerical simulation study of boiling fl ow by Juric and Tryggvason (1998), Brackbill et al. (1998) and is re viewed in Tryggvason et al . (2001). The mass and energy conservation equations then become: .lvk tuuunxxds C (4.31) ..P Pk tCT uCTkTmxxds t C (4.32) 4.1.5 Nondimensionalization of the Governing Equations The governing equations with no phase ch ange are nondimensionalized based on the following nondimensionalization: * x x L *u u U *tU t L * 2 ref p p U (4.33) * refg g g * ref * ref * ref (4.34) * refk k k *ref p p pC C C *cold hotcoldTT T TT (4.35) where L is the characteristic length scale, U is the reference speed, is the density, is the dynamic viscosity, is the surface tension, g is the gravitational acceleration, Cp is
74 the specific heat, k is the thermal conductivity and the thermal diffusivity, the subscripts ref is for the reference values. The re ference properties can be taken, for example, as the liquid properties. The nondime nsional governing equatio ns with no phase change are: .0 u (4.36) 111 .. Rek tu uupunxxdsg tWeFr C (4.37) 1 ..P PCT uCTkT tPe (4.38) in which the superscripts * are omitted for convenience, and where Re is the Reynolds number, We is the Weber number, Fr is the Fr oude number, and Pe is the Peclet number defined as: Reref refUL inertiaforce viscousforce (4.39) 2eref refUL inertiaforce W surfacetensionforce (4.40) 2 refUinertiaforce Fr gLgravitationalforce (4.41) refrefp refCUL convectioneffect Pe kconductioneffect (4.42) Based on Re, We, Fr and Pe, one can also de duce other dimensionless parameters such as the Capillary number (Ca), the Ohnesorg e number (Oh), the Bond number (Bo) and the Prandtl number (Pr).
75 Reref refU Weviscousforce Ca surfacetensionforce (4.43) Reref refCaviscousforce Oh surfacetensionforce L (4.44) 2 refref refgL Wegravitationalforce Bo F rsurfacetensionforce (4.45) Pr Rerefrefp ref refrefC Pemomentumdiffusivity kthermaldiffusivity (4.46) 4.1.6 Dimensionless Parameters within the Micro-Scale Cooling Device Here in the drop dynamics problem, we choose the droplet diameter d as the characteristic length scale, a nd the ejection or impact speed U as the characteristic velocity scale. The time scale is taken as td= d /U. First, a simple computation can be made to evaluate the impor tance of the surface tension force versus gravity within the microscale device. For example, let us consider two water droplet sizes, one of diameter 20 m and the other one of diameter 1 mm. With =997 kg/m3, g =9.81 m/s2, =0.072 N/m, the Bond number is 5.4 x 10-5 for d =20 m, and 0.135 for d =1 mm. This shows that the surf ace tension forces are dominant. Therefore, again, neglecting grav ity is a valid approximation. The liquid and vapor properties of water and R-12 at 1 atm are liste d in Table 4.1. The liquid to vapor ratio of flui d properties are shown in Tabl e 4.2. For the value of the contact angle, different values and models will be selected in chapter 6 to study the effect of the liquid wettability.
76 Table 4.1: Properties of saturated liq uid-vapor at 1 atm from Carey (1992) Fluid Water (H2O) R-12 (CCl2F2) Tsat (K) Tsat ( C) 373.15 100.00 243.20 -29.95 (kJ/kg) 2257 168.3 (N/m) 58.91 x 10-3 15.5 x 10-3 l (kg/m3) 957.9 1486 v (kg/m3) 0.597 6.33 l (N s/m2) 2.78 x 10-4 3.73 x 10-4 v (N s/m2) 1.26 x 10-5 1.03 x 10-5 (Cp)l (kJ/kg K) 4.22 0.896 (Cp)v (kJ/kg K) 2.03 0.569 kl (W/mK) 0.679 9.5 x 10-2 kv (W/mK) 2.5 x 10-2 6.9 x 10-3 Prl 1.72 3.51 Prv 1.02 0.85 Table 4.2: Liquid to vapor ratio of the fluid properties for water and R-12 at 1 atm and at saturation temperature. Properties ratios Water R-12 Density 1600 235 Viscosity 22 36 Specific heat 2.1 1.6 Conductivity 27 14 The dimensionless parameters ranges related to the drop dynamics within the cooling device are shown in Table 4.3 fo r drop diameter ranging from 20 m to 1 mm and for U=1m/s. Table 4.3: Range of dimensionless paramete rs in droplet ejection/impingement problem Water R-12 U=1 m/s U=1 m/s 20 m < d < 1 mm 20 m < d < 1 mm 75 < Re < 3400 80 < Re < 4000 0.32 < We < 16 1.92 < We < 96 Ca=0.005 Ca=0.024 O( 10-5) s < td < O( 10-3) s O( 10-5) s < td < O( 10-3) s
77 It is clear based on the above discussion, the interplay of capillarity, liquid and vapor phases, viscosity, conductiv ity, and convection forms a complex system with the following dimensionless paramete rs: Reynolds, Weber, Froude, Peclet, or others deduced from their combinations. Next , the computational methods fo r solving the field equations as well as tracking the interface are presented. 4.2 Computational Methods The interfacial dynamics is handled based on the immersed boundary technique originally developed by Peskin (1977). It is a mixed Eulerian-Lagrangian method, in the sense that the flow equations are solved using an Eulerian approach on a fixed Cartesian grid, and the interface is repr esented by a discrete set of points and is advected in a Lagrangian framework. In this method, a si ngle set of conservation equations valid for both phases is solved. The inte rfacial conditions are included in the conservation equation as source terms, by converting the surface tension force to a volume body force. The interface is small non-zero thic kness within which the values of the fluid change smoothly. Immersed boundary techniques or fr ont tracking methods have been employed for a variety of two-fluid problems by U nverdi and Tryggvason (1992), Shyy et al. (1996), Udaykumar et al. (1997), Juric and Tr yggvason (1998), and Ka n et al. (1998). For the present immersed boundary technique, th e methodology of Udaykumar et al. (1997) and Kan et al. (1998) is adopt ed. Unlike the approaches take n in these studies, the field equations are solved by a s econd order accurate projection method originally formulated by Chorin (1968) and recently discussed by Ye et al. (1999) instead of the iterative pressure-based algorithm, SIMPLE, of Pa tankar (1980). A time marching pressurecorrection algorithm such as the projection method is more time efficient to solve unsteady flow problem than th e SIMPLE algorithm. As all pressure-based algorithms, the
78 velocity and pressure are decoupled by intr oducing a pressure eq uation that is of Poissonâ€™s type. Because the pressure equa tion is computationally costly, a multigrid technique is adopted to accelerate the rate of convergence (Shyy, 1994; Udaykumar et al., 2001). It is noted that other transport equati ons, being of convectiv e-diffusive type are usually faster to converge in linearized form . In these equations, the nonlinearity from convection and coupling w ith the pressure field are issues which need due considerations. The finite-volume method is employed to disc retize the flow equations. In this study, a computer code written in FORTRAN, specifi cally dedicated to so lve two-dimensional planar and axisymmetric two-flui d flow with heat transfer, is developed. In this chapter, the technical aspects, related to the projection method and the immersed boundary technique, are presented along with validation exercises. 4.2.1 Fluid Flow Field Equation Solver 126.96.36.199 The projection method The projection method, or the fractional-step method, originally formulated by Chorin (1968) is a time-marching pressure-correc tion algorithm for solving the unsteady, incompressible Navier-Stokes equations. The projection method consists of splitting the solution procedure into distinct steps. The present approach follows the work of Ye et al. (1999), which is based on the previous work s of Kim and Moin (1985) and Zang et al. (1994). The continuity and momentum equations are discretized using the finite-volume formulation on a Cartesian, a cell-centered, collocated grid arrangement for primitives variables (u, p ). A face-cell velocity variable U is introduced for computing the volume flux in each cell and to facilitate the enfo rcement of mass conservation. The current cell
79 arrangement is illustrated in Figure 4.1. Th e advantage of using such arrangement is discussed in Ye et al. (1999). V u p v U Figure 4.1: Schematic of a computati onal cell arrangeme nt. The variables p , u, v , are defined at the cell-center, and U , V are the face-c ell velocities. In the first, or projector, step, the moment um equation without the pressure term is solved for an intermediate solution of the velocity *u . A second order accurate AdamsBashforth scheme is employed for the convec tion terms, and Crank-Nicolson scheme for the viscous terms. The semi-discrete projector equation is: 1 223 22* * nn n nUuUu uuuu t (4.47) where is the fluid density, u is the velocity at face-center, and U is the velocity at face-cell, is the dynamic viscosity, t is the time step size, the superscripts n, n-1 denote the current and previous time st ep, respectively, and the superscript * denotes the intermediate step. The veloci ty boundary condition applied at the projection step is 1 nuu*. Once the intermediate face-centered velocity *u is obtained, the intermediate face-center velocity U* is computed by interpolating *u . The second step is the pressu re-correction step. Using:
80 u p t (4.48) and invoking the mass conservation: 10nU (4.49) the following pressure Poisson equation is derived and solved for p : 111n p U t * (4.50) For the problem to be well posed, the velocity component normal to the boundary should be specified which imposes 10npn . Once the pressure is obtained, both the cel l centered and face cell velocities are corrected as: 111nn ccuutp * 1*11nn f cUUtp (4.51) where the subscripts cc and fc denote cel l-centered and face-cell, respectively. The complete description of the above form ulation is given in Ye et al. (1999). 188.8.131.52 Finite-volume discretization The finite-volume method is employed to discretize the governing equations. The detailed discretization of a convection-diffusi on equation is presented here as an example (Patankar, 1980; Versteeg and Malalasekera, 1995). Note that the pr edictor equation (Eq. 4.47) and the energy equatio n are both an unsteady conve ction-diffusion problem and that the pressure Poisson equation (Eq. 4.50) is a diffusion problem. The general unsteady convection-diffusion equation can be written for any field variable as: uS t (4.52)
81 and in integral form as: ..AAdundAndASds t (4.53) where is the volume, A is the surface, n is a unit vector normal to the surface of the control volume, and is the diffusion coefficient. The discretization notation is illustrated in Figure 4.2 to help understand the following discussion. The cell centered at point P, is the reference cell on whic h the field equations are discretized. x )e x )w y )n y )s x y w s e P n E W N S Figure 4.2: Illustration of the discretization notations. The discretized form of Eq. (4.52) using the Adams-Bashforth scheme for the convection terms and the Crank-Nicolson scheme for the diffusion terms is: 131 22nnn nHH x yDDSxy t (4.54) where H represents the convection terms and D the diffusion terms. H is expressed as: eewwnnssHFFFF (4.55)
82 where Fe, Fw, Fn, Fs are the convective normal fluxes through each cell face of the control volume respectively in the east , west, north and south direction: eeeFuy wwwFuy (4.56) nnnFvx sssFvx (4.57) where ue, uw, vn, vs are the velocity at the face cell: 1 2eijEPuUuu , 11 2wijPWuUuu, (4.58) 1 2nijNPvVvv , 11 2 s ijPSvVvv, (4.59) The diffusion term D is: eEPwPWnNPsPSDDDDD (4.60) where e e ey D x w w wy D x (4.61) n n n x D y s s s x D y (4.62) The source term is CPPSSS (4.63) The final discretizated form is cast in the following form: PPEEWWNNSSaaaaab (4.64) where 1 2 E eaD 1 2WwaD (4.65)
83 1 2 N naD 1 2SsaD (4.66) 1 031 22nn nn PPCHH bDaSxy (4.67) 0 PEWNSPPaaaaaaSxy (4.68) 0 n P Paxy t (4.69) Equation (4.64) consists of a five-point stencil and s econd-order central difference schemes are adopted for the convection and diffusion terms. In matrix form, Equation (4.64) reads:  ab (4.70) where [a] is a pentadiagonal coefficient matrix, is the solution vector and b is the source vector. The system is solved in an iterative procedure using line Successive Over Relaxation (LSOR) technique that decompos es the system in two tridiagonal matrix systems. The tridiagonal systems are so lved using the Thomas algorithm. The axisymmetric formulation requires s light modifications in the discretization (Ferziger and Peric, 1996). The x-axis is the axial direction and y-axis is the radial direction. To consider the radius effect, the radii at the node P, east, west, north and south faces are defined as: 11 2Pewjjrrryy (4.71) njrr 1 s jrr (4.72) The axisymmetric volume of the cell is Pdrxy (4.73)
84 and the expressions of the face area are: eedAry wwdAry (4.74) nndArx ssdArx (4.75) In the momentum equation in the radial direction (or v-momentum) an additional term appears: 2Pv d r (4.76) which is treated as a source term. 184.108.40.206 Multigrid technique for the pressure equation A multigrid technique is employed to accelerate the rate of convergence of the pressure Poisson equation. In multigrid tec hniques, computations are performed on fine and coarser grid sizes following a certain schedule or cycle order. The idea behind multigrid is to take advantage of the fact that high wave number components decay faster than low wave number components. Whether we consider a wave number high or low depends on the mesh size. Fo r example in the course of iterations, a low wave number component in a fine mesh system behave s like a high wave number component in a coarse mesh system (Shyy, 1994). Consequently , it is advantageous to treat various wave number components on differe nt grid levels to acceler ate the convergence rate. The multigrid procedure adopted in this work is based on a standard multigrid Wcycle and follows that adopted in Udaykuma r et al. (2001), summarized below. Consider a sequence of g=1,,G grids, where the grid spacing hg-1 on grid g-1 is twice of the grid spacing hg on the grid of level g, i.e. hg-1=2hg. The final discretized form of the pressure equation simplifies to the following system:
85 GGGab (4.77) where [a]G is the coefficient matrix, G is the solution vector and Gb is the source vector on G, the finest grid level. The final converged solution is of course obtained on the finest level G, but intermediate computations are performed on coarser grid levels. First, few calculations are performed on the finest grid G. g g ggabR (4.78) where g R represents the residual vector. The residuals are then transferred by a restriction operator to the next coarser grid: 1 g gM M g RR (4.79) where M denotes the four surrounding points on th e fine grid. The following system is then solved on the coarse grid: g-111 g g gaR (4.80) The solution vector is then transferred by a prolongation operator to the fine grid: 1 1 g gMM M g (4.81) where M in this case denotes the four su rrounding points on the coarse grid and M denotes the interpolation we ights. Finally the solution vector is corrected as: 1 g ggg (4.82) To illustrate the acceleration of the conve rgence rate and theref ore the overall speedup of the computations, a 2D lid driven square cavity flow for Re=100 is considered. The velocity is specified on the four walls: on the top u=1.0, v=0.0 and on the bottom, left and
86 right walls u=v=0.0. The present results are in agreement with those obtained by Ghia et al. (1982) using a vorticity-stream function fo rmulation. The streamlines at steady state are shown in Figure 4.3 to illustrate the pr oblem. Figure 4.4 shows the convergence path on the pressure equation. The pressure equa tion maximum residual versus the number of fine grid iterations to reach a convergence of 10-8 at the first time st ep iteration for a 102 x 102 grid is plotted in Figure 4.4 for one, tw o and three multigrid level. Utilizing twoand three-grid levels clearly enhances the convergence rate. Table 4.4 presents the work units to reach the residual level of 10-8 for the pressure equation and total CPU time to solve the flow equations normalized by the tota l CPU time of single grid computations on the 102 x 102 grid, All data are based on the comput ations conducted at the first time step of the method, with zero values assigned as the initial guess. A 102 x 102 grid with 3 levels and a 202 x 202 grid with 4 levels are consider ed. With three levels on the 102 x 102 grid, it takes 10% of the CPU time needed for the single grid computations to reach a convergence level of 10-8. As the mesh is refined, the grid size increases and the advantage of the multigrid technique grows.
87 X Y 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Figure 4.3: Streamlines plot at steady state of the cavity li d driven flow for Re=100 on a 102 x 102 grid. Figure 4.4: Convergence path of the pressure e quation for a cavity lid driven flow for Re=100 at the first time step. The finest grid is 102 x 102. 0 500 1000 1500 10 -10 10 -8 10 -6 10 -4 10 -2 10 0 Fine Grid Iterations Residual 2 levels 1 level 3 levels
88 Table 4.4: Multigrid performance for a cavity lid driven flow for Re=100. 102 x 102 202 x 202 Levels Work units Relative CPU time Work units Relative CPU time 1 1299 1.0 4702 15.78 2 383 0.32 1364 4.82 3 91 0.10 288 1.12 4 89 0.43 220.127.116.11 Single phase flow examples To validate the single phase flow solver two benchmark problems are considered. The first problem is the buoyancy-induced flow in a square cavity w ith vertical walls differentially heated. This case demonstrates the coupling of the energy equation with the flow equations. The second test case is the Po isseuille flow, i.e., a flow in a tube of circular cross section to validat e the axisymmetric formulation. Buoyancy-induced convection. The problem of buoyancy-driven flow in a square cavity with vertical walls differentially he ated has been presented by Jones (De Vahl Davis and Jones, 1983) as a benc hmark case for code validati on. In addition to the study by De Vahl Davis and Jones ( 1983) additional results have b een reported, e.g. Le Quere (1991), Shyy (1994), Ferziger and Peric (1996) and Wan et al. (2001). The governing equations are the coupled continuity, mome ntum and energy equations. The fluid is assumed incompressible and the Boussinesq approximation is adopted for the density treatment. The nondimensional parameters are: * x x L *y y L *uL u *vL v * 2t t L 2 * 2 p L p *cold hotcoldTT T TT (4.83) 3hotcoldgTTL Ra Pr (4.84)
89 The nondimensional equations, omitting the superscript * for convenience, are (Wan et al. 2001): 0uv xy (4.85) 22 22Pruuupuu uv txyxxy (4.86) 22 22PrPrvvvpvv uvRaT txyyxy (4.87) 22 22TTTTT uv txyxy (4.88) The top and bottom walls are adiaba tic. The left wall is held at Thot=1.0 and the right wall at Tcold=0.0. The boundary condition for th e velocity is no-slip on the four walls. The computations are made for Pr=0.71 and Ra=104, 105, and 106, all in the laminar flow regime. The following runs are performed: Ra=104 5,000 time steps of t=10-3 on uniform grid 82 x 82, and non-uniform grid 114 x 114 10,000 time steps of t=5 x 10-4 on uniform grid 162 x 162 Ra=105 50,000 time steps of t=10-4 on uniform grid 82 x 82, and non-uniform grid 114 x 114 100,000 time steps of t=5 x 10-5 on uniform grid 162 x 162 Ra=106 500,000 time steps of t=10-5 on uniform grid 82 x 82, and non-uniform grid 114 x 114 The time step is chosen smaller for larger Ra since the velocity is known to increase with Ra from previous results. The initial solution for the computations for Ra=105 is based on the final solution of Ra=104, and similarly, the initial solution for Ra=106 is based on the
90 solution of Ra=105. The non-uniform grid 114 x 114 is shown in Figure 4.5. The mesh is refined along the wall in order to improve the solution accuracy without having to compute on a very fine mesh everywhere. The streamlines and isotherms are shown in Figure 4.6 for the computations made with the grid 114 x 114. As Ra increases, one observes (i) the skew symmetry solutions of the velocity and temperature field with respect to the cavity cen ter, and (ii) higher heat transfer rate along the vertical walls as deduced from closer isotherms. The followi ng results are computed from the solution: the maximum and mid-point valu es of the streamfunction, |max|, |0.5|, where u y the maximum vertical velocity on the horizontal mid-plane (on y=0.5), Vmax the average Nusselt number along th e vertical plane of coordinate x, 1 0 xxPNuxNudy where PrxPT NuRauT x the maximum and minimum local Nu sselt numbers on the hot wall (x=0, u=0), 0 xT Nu x The Nusselt number is calculated using the trapezoidal integration scheme. The results are shown in Table 4.5 for Ra=104, 105, 106 and agrees with the benchmark results of De Vahl Davis and Jones (1983). As one observes, the present solution approach the benchmark ones as the mesh is refined. Fu rthermore, the results obtained on the nonuniform grid are consistent with the ones obtained on the uniform fine mesh, demonstrating that mesh refinement along the wall improve the results accuracy as expected.
91 NonUniformGrid (114x114) Figure 4.5: Illustration of the non-uniform grid 114 x 114. Ra=104 Ra=105 Ra=106 Figure 4.6: Natural convection patter ns for Pr=0.71 obtained with a 114 x 114 grid. Top row: streamlines. Bottom row: isotherms.
92 Table 4.5: Natural convection results for Ra=104, 105, 106 compared to benchmark results. Case Ra=104, Pr=0.71 Method Navier-Stokes equa tions and energy equation with projection step method Stream functionvorticity Present Grid or Reference Uniform 82 x 82 Non-uniform 114 x 114 Uniform 162 x 162 De Vahl Davis and Jones (1983) |max| 5.074 5.075 5.075 5.071 |0.5| 5.074 5.075 5.075 5.071 Vmax 19.630 19.626 19.624 19.617 Nux0.5 2.247 2.246 2.246 2.243 Nux0 2.250 2.247 2.247 2.238 Nux0 max 3.542 3.536 3.536 3.528 Nux0 min 0.496 0.550 0.557 0.586 Case Ra=105, Pr=0.71 Method Navier-Stokes equa tions and energy equation with projection step method Stream functionvorticity Present Grid or Reference Uniform 82 x 82 Non-uniform 114 x 114 Uniform 162 x 162 De Vahl Davis and Jones (1983) |max| 9.625 9.618 9.619 9.612 |0.5| 9.127 9.119 9.119 9.111 Vmax 68.580 68.654 68.653 68.59 Nux0.5 4.540 4.528 4.526 4.519 Nux0 4.546 4.528 4.528 4.509 Nux0 max 7.834 7.747 7.747 7.717 Nux0 min 0.634 0.659 0.688 0.729 Case Ra=106, Pr=0.71 Method Navier-Stokes equations with projection step method Stream functionvorticity Pseudo-spectral Chebyshev algorithm Present Grid or Reference Uniform 82 x 82 Non-uniform 114 x 114 De Vahl Davis and Jones (1983) Le Quere (1991) |max| 16.980 16.829 16.750 16.81 |0.5| 16.672 16.333 16.32 16.386 Vmax 218.70 220.15 219.36 220.559 Nux0.5 9.021 8.788 8.799 8.825 Nux0 8.986 8.882 8.817 8.825 Nux0 max 18.659 17.798 17.925 17.536 Nux0 min 0.930 0.941 0.989 0.979
93 The above results validate the ability of th e present numerical method to solve coupled fluid dynamics heat transfer pr oblem. Next, the Poisseuille flow is considered, to verify the axisymmetric formulation. Poisseuille flow. Posseuille flow is a flow in a tube of circular cross section under the action of a pressure gradient. This problem is a classical problem in fluid mechanics and was first studied by Hagen in 1839 and by Poiseuille in 1840 (B atchelor, 1967). The exact solution to the Navier-Stokes equations exists for this flow problem. It is well suited for validating the axisymmetric formula tion of the code and to estimate the global order of accuracy of the single phase flow solver. Because of the axial symmetry, the flow is unidirectional in cylindr ical coordinates. Therefore, the velocity is only a function of the radial direction. The ex act solution for the velocity is: 224 R r p ur L () (4.89) where p is the pressure difference between two point s of distance of L, R is the radius of the tube cross section, and is the fluid viscosity. The tube considered for the present computa tions is 4 units in length and 0.5 unit in radius. At the inlet, the velocity profile is specified based on the exact solution. A pressure gradient, p/L, of 3.2 is selected. At r=0, i. e., on the centerline of the tube, the symmetry boundary condition is imposed, and at r=R, i.e., on the tube wall, the no-slip condition is applied. At the outlet, the velocity is corrected by en forcing global mass conservation. The Reynolds number (Re=u(0)R/) is 100. Initially, the velocity and pressure are set to zero. Four computations are performed for four different grid spacing h=x=y=0.1, 0.05, 0.025, 0.0125 to determine the order of accuracy based on the 2-
94 norm error. The computations are performed for 3000 time step s with a sufficiently small time step t=0.005, to ensure a fully developed flow everywhere in the tube. The velocity profiles at the outlet are compared to the exact solution and are shown in Figure 4.7. Only results obtained with three grid spacing ar e presented, since it becomes difficult to observe difference in the velocity profiles fo r h smaller than 0.05. The 2-norm error is defined as 2 2 11xyNN numericalexact jj j xyEuu NN (4.90) where Nx, Ny are the number of grid points in the axial and radi al direction respectively, and j is the index of the grid cell. The results are plotted on a logarith mic scale in Figure 4.8 with a reference line of slope close to 2, demonstrating a second order accuracy. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 u(r)rVelocity Profiles exact h=0.1 h=0.05 h=0.025 Figure 4.7: Velocity profile for Poisseuill e flow, comparison of exact solution with numerical solution.
95 10-2 10-1 10-4 10-3 10-2 10-1 Grid spacingError 2-norm 2nd order Figure 4.8: Error 2-norm of th e radial velocity versus grid spacing for Poisseuille flow. The above investigations show that th e coupled thermo-flu id problem and the axisymmetric formulation have been correc tly implemented in the code and that the single flow solver is second-order accurate. Next, the two-phase flow formulation is presented with an emphasis on the trea tment of interfacial characteristics. 4.2.2 Immersed Boundary Technique An illustration of a typical computational domain composed of two immiscible fluid is given in Figure 4.9. The squa re fluid domain denoted by is covered with a fixed Cartesian grid. The interface between Flui d 1 and Fluid 2 represented by the curve C is present in the fluid domain and is marked by particles (dot s) that do not coincide with the grid points.
96 Figure 4.9: Illustration of a computational domain composed of two immiscible fluids. 18.104.22.168 Interface tracking The interface tracking procedure employed fo llows the work of Udaykumar et al. (1997). The immersed boundary or interface denoted by tC (a curve in 2D or a surface in 3D) is represented by K marker or interfacial points of coordinates k x s with k=1,2,,K. The markers are uniformly distributed along tC at some fraction of the grid spacing, 0.5h
97 2222ss ssssyx n x yxy (4.91) The curvature is then obtained by taki ng the divergence of the normal vector: . n (4.92) The curvature for the 2D planar case is: 3/2 22 s sssss ssyxyx xy (4.93) and for the axisymmetric case is: 3/23/2 2222 s ssssss ssss x yxyx yxyxy (4.94) The derivatives xs,xss,ys, and yss are computed by taking the a ppropriate derivatives of the piecewise polynomials. The normal and curvatur e are evaluated for each marker points. Assignment of the material properties in each phase. With the interface location known with respect to the grid, the material properties are a ssigned in each fluid based on a Heaviside step function: 121k x x H (4.95) where is any material property such as density, dynamic viscosity, k thermal conductivity, and Cp specific heat. The subscripts 1 and 2 denote Fluid 1 and Fluid 2, respectively. And k x xHis the discrete Heaviside step function defined as: 111 1 2 1 0mmmm kk k m k k kxxxx ifxxd dd xx ifxxd ifxxd dimsin H (4.96)
98 where dim is the space dimension, d=2h with h is the grid spacing, x is the grid coordinate, and k x is the interfacial marker coordinate. The discrete Heaviside function in 1D is plotted for different gr id spacing in Figure 4.10. The Heaviside step function allows a smooth distribution of the material properties over a transition zone d. The transition zone here is taken to be twice the grid spacing (d=2h). Thus the interface is of finite thickness d. This approach provides numerical st ability and smoothness; and with a fine mesh around the interface, sati sfactory accuracy can be attained. The above discrete Heaviside step function has been previo usly used by Sussman et al. (1994) and Udaykumar et al. (1997). Once the material properties have been assigned to each fluid, the fluid flow equations are solved. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Heaviside function in 1D xH(x) h=0.05 h=0.01 h=0.005 Figure 4.10: Discrete Heaviside step function in 1D for three different grid spacing with the interface located at x=0.5.
99 22.214.171.124 Immersed boundary treatment As already mentioned, the immersed boundary technique, originally formulated by Peskin (1977) has been empl oyed for a variety of two-fluid or twophase flow problems by Unverdi and Tryggvason (1992), Udayku mar et al. (1997), Juric and Tryggvason (1998). The immersed boundary technique is a technique to communicate the information between the moving Lagrangian interface and the fixed Eulerian grid. In particular, since the locations of the marker po ints do not coincide with the underlying grid points, the velocity field stored at the ce ll-center of the Cartesian grid system is interpolated to obtain the velocity of the interface points. As will be detailed next, the interface force acting on the marker points is spread to th e nearby grid points vi a a discrete Delta function. The discrete Delta function is given by: dim 11 1cos 2 0mm k k m kxx ifxxd dd xx otherwise (4.97) where dim is the space dimension, d=2h with h the grid spacing. In other words, the delta function spreads over 4h. The delta function is nothing but that the derivative of the Heaviside step function. In Figure 4.11, the discrete Delta function is plotted for three different grid spacing to illustrate the eff ect of grid refinement on the interfacial treatment.
100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Delta function in 1D xdelta(x) h=0.05 h=0.01 h=0.005 Figure 4.11: Discrete delta f unction in 1D for three diffe rent grid spacing with the interface located at x=0.5. Estimation of the interfacial force. The interfacial force F exerted by the interface on the flow is incorporated in the momentum equation by means of an integral source term and is expressed as: k tFnxxds C (4.98) where the integral is taken over t C, is the surface tension, is the curvature, n is the normal unit vector, k x x is the â€˜Diracâ€™ Delta function, x gives the location of the grid points, k x gives the position of the marker s as a function of the arclength s. In the discretized form, the force is estimated by: P kkkk k F nxxs (4.99)
101 The force at point P, P F , is simply the sum of the interf acial force of the marker points located inside a ci rcle of radius 2h weighted by the Delta f unction. This procedure is illustrated in Figure 4.12. Circle of radius 2 h around P P h Circle of radius 2 h around P P h Figure 4.12: Marker points considered for the estimation of the force at point P. Estimation of the interface velocity. The velocity of a marker point or interface velocity denoted by kV should satisfy the c ontinuity condition: kkVuxxdx (4.100) where the integral is over the entire fluid domain , u is the fluid velocity vector, and k x x is the Delta function. And in discretized form, the velocity of a marker point is: 2kijk ijVuxxh (4.101) where h is the grid spacing, i and j are the grid location indices, and u the fluid velocity. The velocity at the kth marker, kV , is computed as the sum of th e fluid velocity of the grid points located inside a circle of radius 2h weighted by the De lta function. This interpolation is illustrated in Figure 4.13. The interpolation may not result in a smooth
102 velocity distribution; therefor e a standard Fast Fourier Transform (FFT) is performed to smooth the distribution (Press et al., 1992). Circle of radius 2 h around marker point hX Circle of radius 2 h around marker point hX Figure 4.13: Fluid points consid ered for the evaluation of th e velocity at marker point X. Advection of the interface. The interface velocity kV is by definition the time derivative of the marker point position: k k x V t (4.102) Therefore, the interface is advected as: 11 nnn kkkxxtV (4.103) 126.96.36.199 Mass conservation of the object The immersed boundary technique is not sufficient by itself to ensure mass conservation of the two phases. This is because the interface is advected with the marker velocity that is computed thr ough interpolation of the velocity field, which can results in loss or gain of mass of an enclosed object. The mass conservation technique of Kan et al. (1998) is adopted here. It is based on the following concept: since there is no phase change, the initial mass of the object s hould remains identic al throughout the time
103 marching computations. The mass of object (drop or bubble) is evaluated as follows for the 2D and axisymmetric computations: 22 Ds M yxdsC (4.104) 2 .2 Axisyms M yxdsC (4.105) Before proceeding to the next time iteration, the mass of the new object is calculated and compared to the initial value. If ther e is a difference with the initial mass, a correction procedure is invoked and the iter ation continues. Specifically, the shape correction in the normal direc tion is updated iteratively by m eans of a bisection method. For details, see Udaykumar et al . (1997) and Kan et al. (1998). 188.8.131.52 Overall procedure The overall procedure is as follows: 1Assignment of the material properties 2Evaluation of the normal and curvature of the interface to estimate the surface tension force 3Computation of the velocity and pre ssure (momentum + pressure equations solved) 4Estimation of the interface velocity a nd advection of th e interface with correction to enforce mass conservation 5March in time, return to step 1 184.108.40.206 Spherical drop in static equilibrium Spherical droplet in static e quilibrium is an ideal test fo r evaluating a two-phase flow solver because it has an exact solution sinc e there is no flow fi eld in the surrounded region. This test has become a standard test to examine the pressure jump and spurious velocity field resulting from the interfacial treatment (Brackbill et al., 1992; Kothe et al.,
104 1996; Scardovelli and Zaleski, 1999; Williams et al., 1999; Bussmann et al., 1999). The exact pressure difference across the interface, p, is 2exact dp R (4.106) where is the surface tension and Rd is the droplet radius. The computational domain considered here is a closed cylinder of height H=1.0, and diameter D=1.0. The droplet is placed at the center of the domain, and has a radius Rd=0.25. The surface tension taken in th e calculation is 1.0, the density ratio is 2.0, and both fluids are assumed inviscid . The exact pressure difference p is 8.0. In the numerical solution, the pressure difference is evaluated as: 1111inoutNN numnn nn inoutppp NN (4.107) where Nin is the number of computationa l cell lying inside the drop, Nout is the number of computational cell lyin g outside the drop, and pn is the pressure at the center of each computational cell. The relative error between the theoretical and computed pressure difference can be expressed in terms of the L2-norm error: 2 2 2 1() 1inN numexact B n inexactpp L Np (4.108) 2 2 2 2 1()inN numexact Kn n exactpp L p (4.109) L2B and L2K are based on the definition used in Br ackbill et al. (1992) and Kothe et al. (1996), respectively. n is the volume of cell n and Nin is the total number of interior
105 cells of the drop. Note by the above definiti ons, that the error norm of Brackbill et al. (1992) does not take into considera tion the axisymmetric formulation. The mesh used in the computation is uniform, i.e., x=y. Five different mesh sizes are considered. The error study is presen ted in the logarithmic plot of the L2 norm error versus the grid size of Figure 4.14. In addi tion to the baseline case, the effect of the density ratio is also inve stigated and the outcome is presented in Table 4.6. As expected, the error decreases as the mesh is refined as shown in Figure 4.14. The order of accuracy is about 0. 3 using the error norm proposed in Brackbill et al. (1992); it is about 1.0 using the error norm proposed in Kothe et al. (1996). As the density ratio increases the error in the pressure evaluation increases as shown in Table 4.6. The error in the pressure evaluation for a density rati o of 1000 is 6%. The pr essure contours and velocity field obtained for a computation with drop/fluid=10, t=10-4, R/x=10 are shown in Figure 4.15. The velocity field is zero everywhere excep t around the interface, because of the continuum su rface tension force model. Table 4.6: Effect of the density ratios fo r a spherical drop in static equilibrium. R x drop f luid num exactp p L2B Norm Error 100 2 0.983 24.8410 100 10 0.962 26.0310 100 100 0.943 27.3810 100 1000 0.939 27.7110
106 100 101 102 103 10-5 10-4 10-3 10-2 10-1 100 Drop in Static Equilibrium density ratio = 2 log(R/Dx)log(L2) 1rst order L2K L2B Figure 4.14: Effect of grid refi nement on the error fo r the validation test case of a droplet in static equilibrium. 0 0 1 1 2 2 3 3 4 5 5 6 6 7 7 X r 0 0.25 0.5 0.75 1 0 0.1 0.2 0.3 0.4 0.5 (a) Pressure Contours X r 0 0.25 0.5 0.75 1 0 0.1 0.2 0.3 0.4 0.5 (b)Velocity Field Figure 4.15: Pressure contour s and velocity field of a drop in static equilibrium. drop/fluid=10, t=10-4, R/x=10. (a) Pressure contours and (b) velocity field.
107 4.3 Summary In this chapter, the theo retical formulation and th e various aspects of the computational methods employed have been described. The latest summarize as: A second order projection method is empl oyed to solve the mass, momentum and energy conservation equations. The equations are discretized using the finite volume formulation. A multigrid technique is implem ented to accelerate the calculation of the pressure equation. The coupling of the ma ss, momentum, and energy equations with the energy equations is validated by the st udy of the natural convection in a square cavity having its vertical walls differentially heated. The axisymmetric formulation and the accuracy using the Poisseuille flow problem. The two-fluid flow problem is numerica lly solved using the immersed boundary technique. The interface is tracked by ma rker points. The surface tension force is modeled as a continuum body force in the momentum equations. The accuracy of the method is investigated for a drop in static e quilibrium. The errors as functions of grid size and the density ratio are assessed. As expected, the present interfacial treatment is essentially a first order technique. The immersed boundary method has the advant age of solving a single set of equation which is simpler to code and less demanding in computation time. The interface is no longer a discontinuity, but is of thin finite thickness. This has the advantage to provide numerical stability and smoothness. On the other hand, the immersed boundary method is limited by a modest order of accuracy and of handling flows with a large density ratio. However, techniques such as adaptive mesh refinement are available to improve spatial resolution. In the following, a brief di scussion of this approach is given. Adaptive Cartesian grid refinement algorithm s have been under development for more than 15 years. (Berger, 1984; Melton, 1996; Aftomis et al., 2000). The algorithm has been included to a front tracking method by Agresar et al. (1998) to study the deformation of cells. In an adaptive grid refinement algorithm, first a coarse grid is generated. Then cells that contain the interface are subdivided into 2m sub-cells with m the space dimension. The parent an d child cells are stored in a 2m-tree data structure. The
108 subdivision can continue to a specified resoluti on level. Once the grid is locally refined special treatment (interpolati on) is required for the discre tization of the flow equation around discontinuous grids. An illustration of a one level grid refinement around an interface is given in Figure 4.16. A main issue of this technique is the enforcement of the conservation laws, i.e., mass, momentum, a nd energy fluxes at the interfaces between fine and coarse meshes. As discussed in Shyy et al. (1998) and Burton and Eaton (2002), various strategies need to be devised in orde r to avoid creating arti ficial sources at the discontinuous grid boundaries. h Figure 4.16: Illustration of one level local grid refinement.
109 CHAPTER 5 5 DYNAMICS OF RISING DROP AND BUBBLE 5.1 Rising Bubble and Drop by Buoyancy The study of rising drop and bubble by buoyanc y in viscous liquid has received considerable attention over the years. Th e bubble as it rises can deform to spherical, ellipsoidal, skirted, spherical cap shap es (Clift et al., 1978) depending on three dimensionless parameters: the Etvs number, Eo 2gd Eo (5.1) the Morton number, Mo 4 23g Mo (5.2) the Reynolds number, Re ReUd (5.3) or combinations of these parameters. The different bubble shape regimes have been classified on a map in terms of Eo, Mo, and Re (Grace, 1973; Clift et al. 1978). Among the experimental studies, insi ghtful visualization can be f ound in Batchelor (1967), Hnat and Buckmaster (1976), and Bhaga and Weber (1981). One of the first numerical st udies is reported by Ryskin a nd Leal (1984) for the steady motion of an axisymmetric bubble rising in a liquid. Their method employs an orthogonal curvilinear grid that conforms to the bubble shape. They have considered the bubble to be a void. The bubble shape is determined based on the norma l stress balance at the bubble
110 interface. Good agreement is obtained be tween their numerical solution and the experiments of Hnat and Buck master (1976). Later, Dandy and Leal (1989) have further developed the method to consid er the bubble/drop flui d flows. They have investigated the effect of the density and visc osity ratio on the bubble/drop shap es and the associated flow structure. Their results are in good agreement w ith the experimental results of Thorsen et al. (1968). In the method of Ryskin and Leal (1984) and Dandy and Leal (1989), the surface of the bubble/drop is considered as a sharp interface. In th e last decade, the motion of bubble and drop due to gravity have been simulated successfully using front tracking methods (Unverdi and Tryggvas on, 1992; Han and Tryggvason, 1999), VOF methods (Gueyffier et al., 1999; Chen et al . 1999; Lrstad et al ., 2002) and level-set methods (Sussman and Smereka, 1997; S on 2001). Recently, Ye et al. (2001) have successfully simulated rising bubble by buoyancy using a sharp interface method, referred as cut-cell approach, in which the governing equations for each phase are solved simultaneously on a fixed Cartesian grid. In this chapter, we demonstrate the present numerical abilities to simulate axisymmetric two-fluid flows for different fl uid properties and flow parameters. In order to do so, the results obtained with the present numerical method are compared with reported numerical and experimental results , including a rising bubble in viscous liquid and a rising drop of oil in water. 5.1.1 Rising Bubble in a Viscous Liquid In this section, the numerical study of the rise by buoyancy of a bubble of lower density (fluid 1) in a conti nuous phase of higher de nsity (fluid 2) is presented. Several computations are performed fo r different density ratios (1/2) and viscosity ratios
111 (1/2). A main purpose is to in vestigate the effect of thes e ratios on the convergence rate of the computations and on the bubble dynami cs. To assess the present computational method, the results are compared with those of Ryskin and Leal (1984) and Ye et al. (2001). The present computational setup is illustrated in Figure 5.1. Fluid 2 is initially a spherical bubble of diameter d=0.2 centered at x=0.2 on the axis of symmetry. The computational domain is a cylinder of height H=10d and diameter D=5d. Because of the axisymmetric geometry, only half of the domain is considered in the computation. On the wall, the no-slip boundary condi tion is applied. The Reynolds (Re), Weber (We) and Froude (Fr) numbers of this se ction are based on the properties of the continuous phase, i.e., Fluid 1: 1 1ReUd (5.4) 2 1Ud We (5.5) 2U Fr gd (5.6) where 1 is the density, 1 is the dynamic viscosity of Fluid 1, is the surface tension, and g is the gravitational acceleration. The charac teristic length scale is the bubble initial diameter d, and the characteristic speed is the taken as Ugd . In the following the characteristic speed is taken as 1.
112 H = 2 . 0 H=2 D=1 axis Fluid 1 Fluid 2 g H = 2 . 0 H=10d D=5d axis Fluid 1 Fluid 2 g Figure 5.1: Schematic of the computationa l setup for a single bubble rising by buoyancy. Initially the bubble (Fluid 2) is a sphere of diameter d=0.2 centered on the axis at x=0.2. Multigrid performance. First, the multigrid performance is tested for a case with Re=100, We=4, density ratio of 10 and viscosity ratio of 1 on a 402 x 72 grid. One, two and three levels are considered. The maximum residual of the pressure equation is plotted versus the number of fine grid iterations in Figure 5.2. Si milar performan ce is observed than for a single phase flow. Therefore, three grid levels are chosen in the computations. Grid resolution effect. The grid resolution is investig ated for a case with Re=100, We=4, Fr=1, density ratio of 100, and visc osity ratio of 1. A coarser grid of 202 x 42 meshes and a finer grid of 402 x 72 meshes are considered. Both grids are uniform in the axial direction and non-uniform in the radial di rection with the mesh being finer closer to the axis of symmetry. Around the bubble, the grid spacing, h, is equal in the axial and radial direction. For the 202 x 42 grid, the number of cells per diameter is d/h=20, and for the 402 x 72 grid, d/h=40. The effect of the grid refinement on rising speed and aspect ratio of the bubble are plotted vers us time in Figure 5.3. The as pect ratio is defined as the
113 bubble height on the axis of sy mmetry divided by the bubble widt h in the radial direction. The difference between the result s on the coarse and fine grid is small. A lthough there is more difference in the aspect ratio result s than in the rising speed, asymptotically, solutions on the two grids converge to each other. The fine non-uniform grid with 402 x 72 meshes is used to study the e ffect of the properties ratios. Figure 5.2: Multigrid performance on the co mputations of a bubble rising by buoyancy. Figure 5.3: Rising speed and as pect ratio versus time for Re=100, We=4, Fr=1, density ratio=100, viscosity ratio=1, on 202 x 42 (d/h=20) and 402 x 72 grids (d/h= 40). x 10 0 1 1.5 2 2.5 4 10 7 10 6 10 5 10 4 10 3 10 2 10 1 0 0.5 7 Fine Grid IterationsResidual 1 level 2 levels 3 levels time 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 rising speed 202 x 42 402 x 72 time 0 0.5 1 1.5 0.4 0.5 0.6 0.7 0.8 0.9 1 aspect ratio 202 x 42 402 x 72
114 Effects of transport properties. The effects of the density and viscosity ratios on both bubble shape and rising spee d are investigated. Comput ations are performed at Re=100, We=4, and Fr=1 for three density ra tios: 10, 100, and 1000 w ith a viscosity ratio of 1, and for three viscosity ratios 1, 10, 100 with a density ratio of 10. The computation of the pressure Poisson equation is most de manding especially with large density ratios. As shown in Figure 5.4, more iterations are re quired for a larger density ratio to achieve the same convergence level. Figure 5.4: Effect of the density ratio (d.r .) on the convergence path of the pressure equation for Re=100, We=4, Fr=1 at the first time iteration on 402 x 72. To present the effects of density and viscosity ratio on the bubble dynamics, the instantaneous bubble shapes with varying dens ity and viscosity ratios are shown in Figure 5.5 and Figure 5.6, respectively. From Figure 5. 5 and Figure 5.6, one notices that with the selected combinations of Re, We and Fr, as the initially spherical bubble rises in the cylinder, it deforms toward a steady-state shap e. The time evolutions of the rising speed and aspect ratio of the rising bubble are plotted in Figure 5.7 and Figure 5.8. As expected 0 500 1000 1500 2000 2500 3000 3500 4000 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 IterationsResidual d.r. = 1000 d.r. = 100 d.r.=10
115 and shown in Figure 5.5 and Figure 5.7 when the density ratio increases, the bubble rises faster. For the three viscosity ratios consid ered, the rising speed of the bubble eventually converges to the same value, as shown in Figur e 5.8. For the cases with viscosity ratio of 10 and 100, the bubble reaches a constant rising speed at an earlier time and deforms more compared to the case with a viscosity ratio of 1. Also, for the cases with viscosity ratio of 10 and 100, the bubble approaches the same shape as shown in the aspect ratio plot. densityratio=10 densityratio=1000 densityratio=100 Figure 5.5: Instantaneous bubble sh apes at time t=0, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5 for Re=100, We=4 and Fr=1 for density ratios of 10, 100 and 1000, and viscosity ratio of 1.
116 viscosityratio=1 viscosityratio=100 viscosityratio=10 Figure 5.6: Instantaneous bubble sh apes at time t=0, 0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5 for Re=100, We=4 and Fr=1 for viscosity ratios of 1, 10 and 100, and density ratio of 10. Figure 5.7: Effect of the density ratio on the rising speed and aspect ratio versus time for Re=100, We=4, Fr=1, viscosity ratio=1. time 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 rising speed density ratio = 10 density ratio = 100 density ratio = 1000 time density ratio = 10 0 0.5 1 1.5 0.4 0.5 0.6 0.7 0.8 0.9 1 aspect ratio density ratio = 100 density ratio = 1000
117 Figure 5.8: Effect of the visc osity ratio on the rising speed a nd aspect ratio versus time for Re=100, We=4, Fr=1, density ratio=10. The streamlines in a referen ce frame moving with the bubble are plotte d in Figure 5.9. A wake structure behind the bubble is clearly observed in Figure 5.9, which is in accordance with experimental observations (Hnat and Buckmaster, 1976; Bhaga and Weber, 1981) and numerical simulations (R yskin and Leal, 1984; Ye et al., 2001). The wake length increases as the density ratio is raised, which can be explained by the shape of the bubble. For a larger density ratio, the steady bubble shape takes a higher aspect ratio. The wake structure behind a bubble is different from the one observed for flow over a solid object with the same shape, becau se the Young-Laplace condition dictates the distribution of the pressure field. time 0 0.5 1 1.5 0 0.2 0.4 0.6 0.8 1 1.2 1.4 rising speed viscosity ratio = 1 viscosity ratio = 10 viscosity ratio = 100 time 0 0.5 1 1.5 0.4 0.5 0.6 0.7 0.8 0.9 1 aspect ratio viscosity ratio = 1 viscosity ratio = 10 viscosity ratio = 100
118 Figure 5.9: Streamlines in a reference fram e attached to the bubble. Re=100, We=4 and Fr=1 for the three density ratios cases at t=1.5. Effects of Re and We. The effects of the Reynolds and Weber numbers on the bubble steady state shape are investig ated next. In Figure 5.10, the computed axisymmetric rising bubble shapes at t=1.5 are represented as functions of Re and We, with Fr=1, the density ratio of 100 and viscosity ratio of 1. Spherical to ellipsoidal shapes are observed as Re and We increases. The computed as pect ratios for the corresponding cases are given in Table 5.1. Since the present numerical method is based on the immersed boundary method, the interface is considered to be of finite thickness. In the present simulations, the pressure jump at the interf ace is not sharp but spr eads over 4 grid cells. Therefore, it is difficult to calculate accura tely the drag coefficient from the pressure distribution on the interface. Another way to evaluate the drag coefficient, CD, is to consider the balance between the buoyancy force and hydrodynamic drag, which leads to the following relation (Bhaga and Weber, 1981; Ryskin and Leal, 1984): 244 33Dgd C UFr (5.7)
119 The drag coefficient is calculated based on the rising speed U at t=1.5. The present results are compared in Table 5.2 to the dr ag coefficient calculated in the numerical studies of Ryskin and Leal (1984) and Ye et al. (2001). Despite the fact that different density ratio and viscosity ratio are employe d in the present simulations and the drag coefficient is evaluated t=1.5, the present results are in reasonabl e agreement with the previous numerical studies of Ryskin a nd Leal (1984) and Ye et al. (2001). We Re 4 10 15 20 20 50 ----100 200 ----Figure 5.10: Computed rising bubble shapes at time t=1.5 as a function of Re and We with a density ratio of 100 and viscosity ra tio of 1. Re and We are based on U=1. Table 5.1: Bubble aspect ratio at t= 1.5 as a function of Re and We. Re We 4 10 15 20 20 0.82 0.65 0.58 0.52 50 0.66 0.45 ----100 0.54 0.39 0.29 0.22 200 0.46 0.30 ----
120 Table 5.2: Comparison of drag coefficients as a function of Re and We. First line: CD computed from present simulation with density ratio of 100 and viscosity ratio of 1. Second line: CD from Ryskin and Leal (1984) fo r a void bubble. Third line: CD from Ye et al. (2001) for density ratio of 1600 and viscosity ratio 22. Re We 4 10 15 20 20 2.32 2.16 2.20 2.88 3.22 3.00 3.24 3.55 3.70 3.61 3.60 3.70 50 1.19 1.23 1.20 1.95 2.77 2.30 ----100 0.95 0.81 0.80 1.77 2.26 --2.44 --3.15 --200 0.89 0.59 --1.91 1.64 ------5.1.2 Rising Oil Drop in Water Next, the rise of S.A.E. 10W oil drops in water by buoyancy is simulated for different drop diameter. The results of the present simu lation are compared with the experimental results of Klee and Treybal ( 1956) and with the results of the numerical simulation of Lrstad et al. (2002). The numerical method em ployed by Lrstad et al. (2002) is a three dimensional VOF technique. Th e fluid properties selected for the simulation are presented in Table 5.3. The f our cases investigated are shown in Table 5.4 along with the corresponding dimensionless numbers, and the estimated terminal velocity using the empirical relation of Klee and Treybal (1956). Table 5.3: Physical properties of S. A.E. 10W oil drops and water at T=23oC from Klee and Treybal (1956) Fluid Properties Water S.A.E. 10W OilRatio Density () [kg/m3] 997.5 865 1.153 Viscosity () [kg/(sm)] 1.06 x 10-372.1 x 10-3 0.0147 Surface Tension () [N/m] 18.5 x 10-3
121 Table 5.4: Computational cases fo r the rising oil drop in water. Case d (mm) U (cm/s) Re wUd/w Eo= d2g/ Fr= U/(gd)0.5 We= wU2d/ 1 2.5 7.42 174.6 0.439 0.4738 0.7421 2 5 10.74 505.3 1.757 0.4849 3.110 3 10 10.74 1011 7.026 0.3429 6.219 4 15 10.74 1516 15.81 0.2800 9.329 Wall effect. Before investigating the effect of the drop diameter, the effect of the computational domain size and the grid size are investigated for a drop of diameter d=5 mm. To investigate the effect of the comput ational domain size, i.e. the effect of the position of the cylinder wall, four cases are considered: (1) D=6d x H=16d, (2) D=8d x H=16d, (3) D=10d x H=16d, (4) D=12d x H=16d. The computed instantaneous drop shapes in two computational configurations D=6d x H=16d and D=12d x H=16d are shown in Figure 5.11. At first glance, th e drop shapes are not affected by the wall positions. A closer investigation on the drop rising speed is considered. The wall effect on the terminal rising velocity can be determ ined with the semi-empirical relation of Harmathy (1960): 210.5udd for UDD (5.8) where u is the terminal rising speed, U is the terminal rising sp eed if the drop is rising in an infinite domain, d is the drop diameter and D is the cylinder diameter. The wall effect predicted with Harmathyâ€™s correlation and the present simulation is compared in Figure 5.12 versus D/d. According to Harmathyâ€™s rela tion, a computational domain of D/d=12 yields a wall effect of around 0.7%. In the com putations, the domain with D/d=12 is taken as the reference. The rising sp eed at t=10 is used for comparison. Due to
122 the viscous effect at the wall, a smaller dom ain reduces the rising speed. This is mainly due to the viscous effect at the wall. Figure 5.12 shows good agreement between Harmathyâ€™s correlation and th e present simulation. D=6dxH=16d D=12dxH=16d Figure 5.11: Effect of wall pr oximity on instantaneous risi ng drop shapes at the same nondimensional time t=0, 2, 4, 6, 8, 10 of a 5 mm oil drop in water. Figure 5.12: Effect of wall proximity on ri sing speed in percentage. Comparison of present simulation with Hama thyâ€™s relation Eq. (5.8). D/d 6 7 8 9 10 11 12 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 3 3.25 Wall effect in % Harmathy's relation Present simulation
123 Grid refinement study. A grid refinement study is now performed. Four grid sizes are considered having respectively 20, 32, 40, and 64 cells per diameter (d/h). The rising speed and the aspect ratio of the drop are pl otted versus time in Figure 5.13 for different grid sizes. The rising speed and aspect ratio at t=10 versus grid si ze is shown in Figure 5.14. From Figure 5.13, one observes that the transient results are not grid independent and that a grid with d/h=20 is too coarse as seen in the aspect ratio plot. The results suggest that for grid size with d/h 32, the rising speed and aspect ratio will eventually converge to the same steady state value. As th e grid is refined, the difference in results decreases as seen in Figure 5.14. Figure 5.13: Drop rising speed and aspect ratio versus time for d=5 mm. Figure 5.14: Effect of grid refinement on ri sing speed and aspect ratio at nondimensional time t=10. 15 25 35 45 55 65 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 d/h rising speed 15 25 35 45 55 65 0.75 0.8 0.85 0.9 d/h aspect ratio 0 2 4 6 8 10 0 0.2 0.4 0.6 0.8 1 time rising speed d/h=20 d/h=32 d/h=40 d/h=64 time 0 2 4 6 8 10 0.75 0.8 0.85 0.9 0.95 1 1.05 aspect ratio d/h=20 d/h=32 d/h=40 d/h=64
124 Effect of the drop size. The instantaneous drop shapes for the different cases investigated are shown in Figure 5.15. As th e drop diameter increases, the Reynolds and Weber numbers increase. Like in the case for bubble, as Re and We increase, spherical to ellipsoidal drop shapes are pr edicted. The terminal rising speed and aspect ratio are plotted versus the drop diameter in Figure 5.16 along with the experimental results of Klee and Treybal (1956) and the numerical results of Lrstad et al. ( 2002). The results are in good agreement for small drop of diameter (d = 2.5 and 5 mm); however there are noticeable differences for larger drops (d=10 and 15 mm). No data for d=15 mm for the experimental data of Klee a nd Treybal (1956) and the numerical simulation of Lrstad et al. (2002) since drop break-up occurs. The pr esent immersed boundary technique predicts greater aspect ratio for the drop diameter d=10 and 15 mm compared to the VOF method of Lrstad et al. (2002) for the same grid spacing d/h=32.
125 d=2.5mm Re=175 We=0.74 Fr=0.47 d=5mm Re=505 We=3.1 Fr=0.48 d=10mm Re=1011 We=6.2 Fr=0.34 d=15mm Re=1516 We=9.3 Fr=0.28 Figure 5.15: Effect of the drop diameter on th e instantaneous shapes of a rising oil drop in water at t=0, 2, 4, 6, 8, 10, 12.
126 Figure 5.16: Rising speed and aspect ratio ve rsus drop diameter. Comparison of present simulation using immersed Boundary method w ith the experimental results from Klee and Treybal (1956) and with the numerical si mulation of Lrstad et al. (2002) using VOF method. The numerical si mulations are performed on a grid with d/h=32. 5.2 Drop Motion within the Micro-Scale Cooling Device To simulate the droplet dynamics within the micro-scale device describe in chapter 2, the problem illustrated in Figure 5.17 is cons idered. The droplet is represented by Fluid 2 and is given an initial velocity U. There is no gravity. The nondimensional parameters are the Reynolds (Re), Weber (We) numbers ba sed on the liquid phase or fluid 2: 2 2ReUd (5.9) 2 2Ud We (5.10) where 2 is the droplet density, 2 is the droplet dynamic viscosity, is the surface tension, d is the droplet diameter and U is the initial velocity. The two above nondimensional parameters can be combined to define the Capillary (Ca) number: 2ReWeU Ca (5.11) The properties of saturated liquid-vapor of water and R-12 are listed chapter 4. drop diameter (mm) 0 2 4 6 8 10 12 14 16 0 5 10 15 rising speed (cm/s) Experiment (Klee and Treybal) VOF (Lrstad et al.) IB (present simulation) drop diameter (mm) 0 5 10 15 0 0.2 0.4 0.6 0.8 1 aspect ratio Experiment (Klee and Treybal) VOF (Lrstad et al. ) IB (present simulation)
127 Numerical simulations are performe d on a non-uniform grid with 202 x 82 meshes. The effect of the droplet diameter and the ini tial velocity on the flow field is investigated for a density ratio of 100 and viscosity ratio of 10. The droplet diam eters considered are: (a) d=0.1, (b) d=0.15 and (c) d=0.2 for constant initial velocity U=1.0 and surface tension =1.0. The corresponding dimensionless parame ters are (a) Re=50, We=10, Ca=0.2; (b) Re=75, We=15, Ca=0.2; (c) Re=100, We=20, Ca=0.2. The instantaneous droplet shapes are shown in Figure 5.18. The larger the droplet is the faster it travels, due to its higher inertia. The initial velocities U considered are: U=1.0 and U=2.0 for constant droplet diameter d=0.1 and surface tension =1.0. The corresponding dimensionless parameters are Re=50, We=10, Ca=0.2, and Re=100, We =40, Ca=0.4. The instantaneous droplet shapes are shown in Figure 5.19. Finally, numerical results for Re=278, We =10, Ca=0.036, density ratio of 235 and viscosity ratio of 36, which corresponds to the case of R-12 , are shown in Figure 5.20. Two cases with a density ratio of 1600 and vi scosity ratio of 22, which correspond to water, are shown in Figure 5.21, with two parameter sets: (a) Re=45, We=10, Ca=0.22 and (b) Re=455, We=10, Ca=0.022. Fluid 2 H=5d D=5d axis Fluid 1 Fluid 2 Figure 5.17: Schematic of the computational setup for the droplet dynamics. Initially the droplet (Fluid 2) is a sphere of diam eter d centered on the axis at 0.2.
128 X Y -0.5 0 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1Re=50We=10Ca=0.2 d=0.10 X Y -0.5 0 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 d=0.15 R e= 7 5 W e= 1 5 C a= 0 . 2 X Y -0.5 0 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1Re=100We=20Ca=0.2 d=0.20 (a) (b) (c) Figure 5.18: Effect of droplet diameter on droplet dynamics. Instantaneous droplet shapes at same instants: t = 0, 0.2, 0.4, 0.6, 0.8, 1.0. (a) d=0.1, Re=50, We=10, Ca=0.2; (b) d=0.15, Re=75, We=15, Ca=0.2; (c) d=0.2, Re=100, We=20, Ca=0.2. Density ratio = 100, viscosity ratio = 10, U=1.0, =1.0. X Y -0.5 0 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1Re=50We=10Ca=0.2 U=1.0 X Y -0.5 0 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1Re=100We=40Ca=0.4 U=2.0 (a) (b) Figure 5.19: Effect of initial velocity on droplet dynamics. Instantaneous droplet shapes for (a) Re=50, We=10, Ca=0.2 at t = 0, 0.2, 0.4, 0.6, 0.8, 1.0; (b) Re=100, We=40, Ca=0.4 at t = 0, 0.2, 0.4, 0.6.
129 X Y -0.5 0 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1Re=278We=10Ca=0.036 densityratio=235 viscosityratio=36 Figure 5.20: Instantaneous shapes for a R-12 dr oplet. Computations with a density ratio of 235, and viscosity ratio of 36, at t=0, 0.2, 0.4, 0.6, 0.8 for Re=278, We=10, Ca=0.036, d=0.1, U=1.0, corresponding to a R-12 droplet of diameter d=46 m with initial velocity U =1.5 m/s, at t=0, 6, 12, 36, 48 s. X Y -0.5 0 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1Re=45We=10Ca=0.22 densityratio=1600 viscosityratio=22 X Y -0.5 0 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1Re=455We=10Ca=0.022 densityratio=1600 viscosityratio=22 (a) (b) Figure 5.21: Instantaneous shapes for water dr oplets. Computations with a density ratio of 1600 and viscosity ratio of 22 for (a) Re=45, We=10, Ca=0.22, d=0.1, U=1.0, at t=0, 0.2, 0.4, 0.6, 0.8, 1.0, corresponding to a water droplet of diameter d=0.28 m with initial velocity 46.6 m/s, at t = 0, 0.0012, 0.0024, 0.0036, 0.0048, 0.006 s, (b) Re=455, We=10, Ca=0.022, d=0.1, U=1.0, at t=0, 0.2, 0.4, 0.6, corresponding to a water droplet of diameter d=28m with initial velocity 4.66 m/s, at t=0, 1.2, 2.4, 3.6 s.
130 5.3 Summary Numerical simulations of rising bubble and dr op are performed to assess the capability and accuracy of the present method. Comparis ons are made with the available numerical and experimental studies. The effects of fl uid properties ratio, grid refinement, and dimensionless parameters are investigated is some detail. Overall, the present numerical method performs adequately for bubble and drop dynamics. Simulations of drop dynamics relevant to the micro-scale de vice application are performed. Within the present flow regime, the drop remains larg ely spherical due to high surface tension effect.
131 CHAPTER 6 6 DROP IMPACT DYNAMICS WITH HEAT TRANSFER In this chapter the dynamics of a drop impacting on a flat solid surface and the associated heat transfer characteristic s are numerically investigated using the computational techniques described in chapter 4. The literature on this subject has been reviewed in chapter 3. As already mentione d, previous numerical studies have mostly treated the problem as a free surface fl ow without accounting for the physical mechanisms of the gas phase. Here, the liquid and gas phases are modeled simultaneously. The governing equations are the mass, momentum and energy conservation equations. The problem is considered to be axisymmetric, and the fluids are Newtonian with constant physical properties in each phase. Drop splashing from one to multiple objects is not studied here. Furthermore, to eliminate the possibility of bubble entrapment, the simulation is ini tiated at the moment when the drop is in contact with the solid surface. To validate the computati ons, selected results are compared with experimental and numerical st udies published in the literatu re. The effects of the liquid transport properties, such as the surface tensio n, viscosity, contact angle and the impact velocity on the drop dynamics are investigated. The heat tr ansport process in the liquid and gas phases during imp act is also studied. 6.1 Initial and Boundary Conditions The computations are performed in a cylindric al domain of diameter and height 5 drop diameter long. An illustration of the initia l and boundary conditions is shown in Figure 6.1. Initially, a spherical droplet of diameter d0 is placed on the surface and given an
132 impact velocity U. The symmetry condition is applied on the cylinder axis and the no-slip condition for the velocity on the cylinder walls is applied. The boundary condition for the contact line position is discussed next. -0.4 -0.2 0 0.2 0.4Y -0.4 -0.2 0 0.2 0.4Z -0.4 -0.2 0 0.2 0.4Y -0.4 -0.2 0 0.2 0.4Z (a) Initial Configuration r x Immersed Interface Fluid 1 = Vapor Fluid 2 = Liquid Cylinder axis n 0 0rxuu 0rxuu 0 T r 0 T r 0rxuu 0 T x nwTT 0x ru u rr (b) Boundary Conditions Figure 6.1: Initial and boundary conditions fo r the droplet impingement problem with heat transfer. (a) Initial confi guration: a droplet of diameter d0 is placed on the surface with given velocity U and (b) boundary conditions. The contact line is located at the intersection between the liquid, gas and solid phase. From a continuum point of view, the contac t line position is described by the contact
133 angle . This angle represents the wettability of a liquid drop on a flat solid surface. This angle varies depending on the fl uid and solid properties, surface characteristics, such as roughness and temperature, and fluid dynamics . A discussion on the contact angle has been given in chapter 3. Analytical me thods based on doma in perturbation and lubrication theory have been developed. These models lead to a singularity at the contact line for Newtonian and incompressible fluids , when the no-slip boundary condition at the wall is applied. The singularity can be rem oved if a slip boundary condition is applied (Dussan V., 1979). In direct numerical simulation of the flow equations involving a contact line, the contact angle is typically prescribed. Either a single contact angle value corresponding to the static contact angl e (Kothe et al., 1991; Brackbill et al., 1992) or a range of dynamic contact angles obtained from measuremen t can be imposed as boundary conditions (Fukai et al., 1995; Bussmann, 1999). Another way to treat the cont act line motion is to consider the classical Young-Dupr relation (E q. 3.2). However, this treatment requires the knowledge of the surface tension coefficien t for each phase, such coefficients are not always available. Renardy et al. (2001) have implemented two different methods to treat the contact line in a two-dimensional VOF numerical method. The first one employs a contact angle value and the second one empl oys the Young-Dupr rela tion. It is found that the second method can introdu ce an artificial flow near the contact point; hence, the first approach is preferable. Renardy et al . (2001) have only performed computational exercises with a density ratio of 1 and 10 and a viscosity ratio of 1 between the liquid and gas, to demonstrate the difference betw een the methods. Their simulation does not represent any specific liquids.
134 The present numerical method does not solve fo r the contact line velocity directly. The contact point position is extrap olated based on prescribed cont act angle, with the no-slip condition on the wall for the fluid velocity imposed. Specifically, the position of the contact line, which is represented by the positio n of the first marker point, is extrapolated from the position of the second marker point and with a specified contact angle. Both static and dynamic contact angle models are considered. For the static contact angle model, a constant value is imposed in the course of the computations. For the dynamic contact angle model, the contact angle va lue is assigned based on the contact line velocity. 6.2 Results with a Static Contact Angle Model 6.2.1 Impact of Water and Ink Drop on Polycarbonate Surface The computations of a water droplet impinging on a polycarbonate surface and the computations of an ink dr oplet impinging on a polycarbonat e surface are now compared with the experiment of Kim and Chun (2001) . A direct comparison between numerical simulation and physical observation of the dropl et shapes at different time instants is shown in Figure 6.2 and Figure 6.3, for water a nd ink, respectively. Despite the fact that a static contact angle model is employed in both cases the results ar e in reasonably good agreement with the experiment. To understand the effect of the different parameters on the drop dynamics and to gain insight on the physical pheno mena involved, a parametric study is presented next.
135 -0.4 ms 0.6 ms 1.6 ms 2.6 ms 3.6 ms 5.6 ms 11.6 ms 14.6 ms 17.6 ms 21.6 ms 0.6 ms 1.7 ms 2.8 ms 3.4 ms 5.7 ms 11.4 ms 14.8 ms 17.6 ms 21.6 ms 0 ms Figure 6.2: Water droplet impinging on pol ycarbonate surface for Re=3200 and We=30. Left column: images from Kim and Chun (200 1) experimental results. Reprinted with permission. Right column: pr esent numerical simulation.
136 -0.1 ms 0.9 ms 1.9 ms 3.9 ms 6.9 ms 7.9 ms 9.9 ms 10.9 ms 13.9 ms 0 ms 0.84 ms 1.96 ms 3.92 ms 7.0 ms 7.8 ms 9.8 ms 10.9 ms 14.0 ms Figure 6.3: Ink droplet impi nging on polycarbonate surface for Re=2300 and We=190. Left column: images from Kim and Chun (200 1) experimental results. Reprinted with permission. Right column: pr esent numerical simulation. 6.2.2 Parametric Variation Effects on the Drop Dynamics The dimensionless parameters for the ba se case of Re=100 and We=20 are based on the following nondimensional scales: d0=0.2, U=1.0, l=100, g=1, l=0.2, g=0.02, =1.0, which correspond for example to the following set of physical parameters: d0=100 m, U=1.0 m/s, l=100 kg/m3, g=1 kg/m3, l=10-4 Ns/m2, g=10-5 Ns/m2, =5 x 10-4
137 N/m. The computational time t=1.0 correspond s to a physical time of 0.5 ms. For the contact angle we have used 60 for our reference computation, which corresponds to the equilibrium contact angle of deionized wa ter on silicon oxide measured by Kim and Chun (2001). However, experimental images of the impingement of a deionized water droplet on silicon oxide were not pr esented in Kim and Chun (2001). 220.127.116.11 Grid size, background fl uid properties effects To demonstrate the effects of the grid size, density and viscosity ratios on the results, three grid sizes with 20, 40, 80 cells per drop diameter, two density ratios of 100 and 1000, and two viscosity ratios of 10 and 100 for constant liquid properties are considered for Re=100, We=20, =60. The effects of the grid size, density and viscosity ratios on the spreading coefficient at t=1.0 and 5.0 are relatively small as shown in Table 6.1, which indicates that the numerical solutions ar e satisfactory. In the following, the effects of the contact angle, impact velocity, droplet viscosity, and surface tension are assessed on grid having 40 cells per drop diameter, a dens ity ratio of 100, and a viscosity ratio of 10. Table 6.1: Effect of the gr id size and background fluid, i.e. gas properties on the drop spreading coefficient at t=1.0 and 5.0 for Re=200, We=20 and =60. Grid size d/h Density ratio (l/g) Viscosity ratio (l/g) at t=1.0 at t=5.0 20 100 10 1.484 1.632 40 100 10 1.513 1.619 80 100 10 1.558 1.613 40 1000 10 1.507 1.618 40 100 100 1.488 1.619
138 18.104.22.168 Effect of the contact angle The contact angle values considered are 60, 80, 100, and 120. These computations are all done for Re=100 and We=20. For va lues of the contact angle between 0 and 90, the liquid is said to wet the surface, and for values of the contact angle between 90 and 180, the liquid is said non-wetting (Carey, 1992). The instantaneous droplet shapes are presented in Figure 6.4 for two contact angles, 60 and 120. For larger contact angle values, the wetting property of the surface di minishes, therefore reducing the contact area or wetting area. From Figure 6.8 (a) one notices for larger contact angle values, both maximum and equilibrium values of the spre ading coefficient decrease and the amplitude of the recoiling oscillations becomes great er, eventually reaching a point beyond which the droplet bounces back, as it happens for computations with =120. 22.214.171.124 Effect of the impact velocity Three impact velocity values U= (i) 1.0, (ii) 2.0 a nd (iii) 3.0 are considered, corresponding to (i) Re=100 We=20, (ii) Re=200 We=80, and (iii) Re=300 We=180, respectively. The corresponding droplet shapes and the spreading co efficient are shown in Figure 6.5 and Figure 6.8 (b) respectively. As expected, for a larger impact velocity, the droplet attains a larger maximum spreads. However, within the range of parameters investigated, the impact velocity has no eff ect on the frequency and the amplitude of the recoiling process. 126.96.36.199 Effect of the droplet viscosity The Reynolds number is varied from 100 to 400 while holding the Weber number unchanged at We=20. The instantaneous drople t shapes are shown in Figure 6.6. As the Reynolds number increases, the maximum spreading increases and the amplitude of the
139 recoiling oscillations increases as can be s een in Figure 6.8 (c). On the other hand, the frequency of the oscillations is not change d. The larger the drop let viscosity is, the greater the recoiling energy dissipates. The vi scous effects clearly dampen the spreading and recoiling of a droplet impinging on a flat surface. 188.8.131.52 Effect of the surface tension Three computations are performed with We=2, 20 and 200 for Re=100 constant, corresponding to =10.0, 1.0, 0.1 respectively. From Figure 6.7 and Figure 6.8 (d), one notices that as the Weber number increases, the maximum spreading coefficient increases while the recoiling frequency decreases. For la rge surface tension, the recoiling process is fast and oscillatory. For lower surface tensi on, the recoiling process is slower and less oscillatory. The larger the surface tension is (lower We), the higher the frequency of the recoiling process is. The qualitat ive trend of the recoiling frequency given in Eq. (3.3) is observed. Surface tension resists drop deforma tion. For higher surface tension, there is more resistance for the drop to deform; theref ore, the drop has a tendency to recoil fast towards successive hemispherical shapes. In contrast, for lower surface tension, a flat disk shape with a ring-like structure in the outer region is maintained for a noticeable period of time. 184.108.40.206 Maximum spreading coefficient A comparison between the values of the maximum spreading coefficient obtained from the present computations and from the energy correlation of Eq. (3.13) is made. As shown in Figure 6.9, the computed maximu m spreading coefficients follow the same trend as the one obtained from the correla tion and are always slightly lower. The maximum spreading coefficient increases, for decreasing contact angl e (Figure 6.8 (a)),
140 for increasing impact speed (Figure 6.8 (b), for increasing Re (Figure 6.8 (c)) and for increasing We (Figure 6.8 (d)). The gap betw een the simulation and correlation results is increasing as the contact angle increases Figure 6.8 (a)), and as the Re increases, i.e. as the viscosity decreases (Figure 6.8 (c)). Th is suggests that the viscous effect term appearing in the derivation of correlation of Eq. (3.13) is no t approximated correctly as also discussed in Middleman (1995).
141 t=0.1 t=0.2 =60 Re=100 We=20 t=0.3 t=0.4 t=0.5 t=0.7 t=1.0 =120 Re=100 We=20 t=0.1 t=0.2 t=0.3 t=0.4 t=0.5 t=0.7 t=1.0 Figure 6.4: Droplet shapes at different time instants for va rying static contact angle.
142 t=0.1 t=0.2 U=1.0 Re=100 We=20 t=0.3 t=0.4 t=0.5 t=0.7 t=1.0 U=2.0 Re=200 We=80 t=0.1 t=0.2 t=0.3 t=0.4 t=0.5 t=0.7 t=1.0 U=3.0 Re=300 We=180 t=0.1 t=0.2 t=0.3 t=0.4 t=0.5 t=0.7 t=1.0 Figure 6.5: Droplet shapes at different time instants for varying impact velocity.
143 t=0.1 t=0.2 =0.2 Re=100 We=20 t=0.3 t=0.4 t=0.5 t=0.7 t=1.0 =0.1 Re=200 We=20 t=0.1 t=0.2 t=0.3 t=0.4 t=0.5 t=0.7 t=1.0 =0.05 Re=400 We=20 t=0.1 t=0.2 t=0.3 t=0.4 t=0.5 t=0.7 t=1.0 Figure 6.6: Droplet shapes at different time instants fo r varying droplet viscosity.
144 t=0.1 t=0.2 =1.0 Re=100 We=20 t=0.3 t=0.4 t=0.5 t=0.7 t=1.0 =10.0 Re=100 We=2 t=0.1 t=0.2 t=0.3 t=0.4 t=0.5 t=0.7 t=1.0 =0.1 Re=100 We=200 t=0.1 t=0.2 t=0.3 t=0.4 t=0.5 t=0.7 t=1.0 Figure 6.7: Droplet shapes at different time instants for varying surface tension.
145 t i m e 0 1 2 3 4 5 6 7 8 9 1 0 0 0 . 5 1 1 . 5 2 2 . 5 " b o u n c e b a c k " E f f e c t o f t h e C o n t a c t A n g l e R e = 1 0 0 W e = 2 0 t i m e 0 1 2 3 4 0 0 . 5 1 1 . 5 2 2 . 5 3 E f f e c t o f t h e I m p a c t V e l o c i t y U=2 U=1 U=3 (a) (b) t i m e 0 1 2 3 4 5 6 7 8 9 1 0 0 0 . 5 1 1 . 5 2 2 . 5 E f f e c t o f t h e D r o p l e t V i s c o s i t y R e = 1 0 0 R e = 4 0 0 We=20 R e = 2 0 0 t i m e 0 1 2 3 4 5 6 7 8 9 1 0 0 0 . 5 1 1 . 5 2 2 . 5 E f f e c t o f t h e S u r f a c e T e n s i o n R e = 1 0 0 We=200 We=20 We=2 (c) (d) Figure 6.8: Spreading coefficient versus time. (a) Effect of the static contact angle for Re=100 and We=20, (b) effect of the impact velocity U for =60, (c) effect of the Reynolds number for We=20 and =60, (d) effect of the surface tension for Re=100 and =60.
146 Contact Angle Effect1 1.5 2 2.5 5060708090100110120130b max Simulation Correlation Impact Velocity Effect1.5 2 2.5 3 0.511.522.533.5 Ub max Simulation Correlation (a) (b) V iscosity Effect 1.5 2 2.5 3 3.5 50150250350450Re b max Simulation Correlation Surface Tension Effect 1.5 2 2.5 3 -50050100150200250Web max Simulation Correlation (c) (d) Figure 6.9: Comparison between the values of the maximum spreading coefficient obtained from the present simulation and from the energy correlation of Eq. (3.13) versus (a) contact angle , (b) impact velocity U, (c) Re, and (d) We. The following insight of the flow dynamics is obtained based on the pressure contours, and streamlines plots (Figure 6.10) at different instants during the drop impact, have given. As the liquid droplet impinges and spreads along the surface, it induces the gas, initially at rest, into motion, and a primary recirculation is shown in the streamline plot. As the droplet reaches its maximum spreadi ng position, a second reci rculation zone of opposite direction of the primary one grows. Th is behavior repeats itself until the droplet reaches its equilibrium position. The jump in pressure at the liquid-gas interface is
147 observed in the pressure contours. The nonuniform pressure distribution inside the droplet also explains the esta blished velocity field. The la rge curvature develops in the vicinity of the contact area causing high pres sure in the area, resu lting in the droplet recoil. As the system tends to equilibrium with time, the pressure gradient in the droplet diminishes. At equilibrium, the hydrostatic pressure jump across the interface is observed.
148 t=0.1 t=0.4 t=1.0 Figure 6.10: Pressure contours and streamlines at different time instants for Re=100, We=20, =60. Left side: pressure cont ours. Right side: streamlines in a fixed reference.
149 6.3 Results with a Dynamic Contact Angle Model 6.3.1 Static versus Dynamic Contact Angle In this section, computations are perfor med for the same case with Re=300, We=60, Fr= with both a static contac t angle model and a dynamic contact angle model and the results are compared. The computational fluid parameters are shown in Table 6.2 and the contact angle values are shown in Table 6.3. This case corresponds to the impact of a 20 m water droplet on a wax coated surface with an impact velocity of 14.6 m/s. The contact angle values have been measured by F ukai et al. (1995) for the case of water on a wax coated surface. For the dynamic contact angle model, the cont act angle value is a ssigned based on the contact line velocity. Fukai et al. (1995) and Bussmann et al. (1999) have employed such dynamic contact angle model based on measured values. In the present simulation, the contact line velocity is interpolated from th e neighboring fluid velocity field. The present dynamic model employed is illustrated in Figur e 6.11. If the instantaneous contact line velocity is within 5% of the impact velocity, then the contact angle varies linearly from the advancing/receding values. If it is greater than 5% of the impact velocity, then the contact angle is assigned the value of the advancing angle (A). Similarly, if the contact line velocity is smaller than -5% of the imp act velocity, the contac t angle takes the value of the receding angle (R). The instantaneous drop shapes, the spreading coefficient and the normalized maximum height are investigat ed as functions of time to assess the difference between the static and dynamic contact angle models. The dimensionless parameters are defined as: *U tt d (6.1)
150 based d (6.2) maxh d (6.3) where is the spreading coefficient, dbase is the length of the dr op on the solid surface, is the normalized maximum height correspond ing to the instantaneous maximum height of the drop, hmax. Table 6.2: Computational fluid properties. Properties LiquidVaporRatio Density () 300 0.3 1000 Viscosity () 0.2 0.01 Surface Tension () 1.0 Table 6.3: Contact angle values of water on a wax coated surface from Fukai et al. (1995). Static (S) Advancing (A) Receding (R) 76 92 60 -0.15-0.1-0.0500.050.10.15V<0V>0 A R Figure 6.11: Illustration of the dynamic contact angle model. Contact angle () in degree versus contact line speed (V). The value 0.05 is taken arbitrary as 5% of the impact velocity U.
151 First a grid refinement study is conducte d. Four grids are considered having respectively 20, 32, 40 and 50 cells per drop diameter (d/h). The effect of the grid size on the spreading coefficient and normalized ma ximum height are shown in Figure 6.12 for the static contact angle mode l and in Figure 6.13 for the dynamic contact angle model. In both cases, the results tend to be grid independe nt as the grid is refined. The results for the coarser grid of d/ h=20 with the static contact angl e model are in disagreement with the other grid sizes, clearl y showing that d/h=20 is too coarse to provide acceptable results. Results obtained with a grid ha ving a minimum of 40 cells per diameter are considered acceptable. The instantaneous shapes of the impinging drop at same instants are shown in parallel for both models in Figure 6.14. The spread ing coefficient and maximum height for both models are compared in Figure 6.15. The cont act angle versus time is plotted in Figure 6.16. Noticeable differences between the two contact angle models in recoiling behaviors are observed. As seen in Figure 6.15, the dynamic contact angle model in comparison to the static contact angle model dampens the osci llations appearing in the recoiling process. This is due to the values of the contact angl e. Smaller contact angle values mean that the fluid wets more the solid surface. This wetta bility property is indeed a resistance to the drop motion. Since the advancing contact angle is greater than the static contact value, the maximum extent of the drop is greater for the static contact angle.
152 Figure 6.12: Effect of grid refinement on sp reading coefficient and maximum height for Re=300, We=60 with the static contact angle model. Four grid sizes are considered with d/h=20, 32, 40, and 50. Figure 6.13: Effect of grid refinement on sp reading coefficient and maximum height for Re=300, We=60 with the dynamic contact angle model. Four grid sizes are considered with d/h=20, 32, 40, and 50. 0 10 20 30 40 50 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t* maximum height d/h=20 d/h=32 d/h=40 d/h=50 Re=300 We=60 dynamic contact angle model 0 10 20 30 40 50 0 0.5 1 1.5 2 2.5 t* spreading coefficient d/h=20 d/h=32 d/h=40 d/h=50 Re=300 We=60 dynamic contact angle model t* 0 10 20 30 40 50 0 0.2 0.4 0.6 0.8 1 1.2 1.4 maximum height d/h=20 d/h=32 d/h=40 d/h=50 Re=300 We=60 static contact angle model 0 10 20 30 40 50 0 0.5 1 1.5 2 2.5 3 t* spreading coefficient d/h=20 d/h=32 d/h=40 d/h=50 Re=300 We=60 static contact angle model
153 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 t*=0.5 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 t*=2.5 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 t*=5.0 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 t*=10.0 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 t*=15.0 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 t*=20.0 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0.2 t*=50.0 Figure 6.14: Instantaneous drop sh apes at same instant for the static contact angle model (left column) and dynamic contact angle mode l (right column) for Re=300, We=60 and a grid with d/h=40.
154 Figure 6.15: Spreading coefficient and normaliz ed maximum height ve rsus dimensionless time t* for Re=300, We=60 with a static contact angle model and a dynamic contact angle model and a grid with d/h=40. Figure 6.16: Contact angle versus dimensionl ess time t* for Re=300, We=60 and a grid with d/h=40. 0 10 20 30 40 5 0 56 60 64 68 72 76 80 84 88 92 96 t* static angle dynamic angle Re=300 We=60 0 10 20 30 40 5 0 0 0.2 0.4 0.6 0.8 1 1.2 1.4 maximum height static angle dynamic angle Re=300 We=60 t * t * 0 10 20 30 40 50 0 0.5 1 1.5 2 2.5 3 spreading coefficient static angle dynamic angle Re=300 We=60
155 6.3.2 Comparison with Experiment In this paragraph, two computational cas es are performed and are compared with existing experimental results from the literature in order to asses the present numerical method with a dynamic contact angle. The cons idered cases are presented in Table 6.4. The impact of a water drop on stainless steel and on wax coated surface is considered. For each of these cases a dynamic contact angle condition is applied. Table 6.4: Cases considered for the compar ison between the present numerical method and experimental results. Case 1 PasandidehFard et al. (1996, 2001) Re=2900 We=47 Fr=86 d0=2.0 mm U=1.3 m/s Water drop Stainless steel surface A=110 S=90 R=40 Case 2 Fukai et al. (1995) Re=6200 We=128 Fr=68 d0=3.72 mm U=1.58 m/s Water drop Wax coated plate A=92 S=76 R=60 A grid size with 80 cells per diameter is utilized because of the large Reynolds number and time steps smaller or equal 10-4 are required. The backgr ound fluid is taken to be 1000 times less dense and 20 times less viscous th an the liquid of the drop. These ratios are close to those of liquid water with steam. Whenever ava ilable, the measured data are used to facilitate the dynamic contact angle model. Note th at the measured contact angle values for Case 1 of the experiment are ta ken from Pasandideh-Fard et al. (1996) for the same liquid and surface studied in PasandidehFard et al. (2001). In Pasandideh-Fard et al. (2001), the advancing contact angle has also been measured using the same technique as in Pasandideh-Fard et al. (1996) and is 11010, however the receding and static contact angles are not reported. In the presen t simulation, an advancing contact angle of 100 is selected. Because of the present mesh resolution, the minimum value that the
156 contact angle can take in the simulation is limited to 60. Therefore, the receding contact angle for case 1 is taken as 60 instead of 40. For Case 2, the measured dynamic contact angle values and the ones from the simulation are the same. 220.127.116.11 Case 1: impact of a water drop on stainless steel A direct comparison of drop images between the present numerical results and the numerical and experimental re sults of Pasandideh-Fard (20 01) is shown in Figure 6.17. Initially, the liquid spreads radially, before it recoils. During the recoiling process, the accumulated mass at the droplet rim moves b ack to the center , which explains the formation of the elevated â€œdomeâ€ on its center at t=7 ms. From the numerical solutions of Pasandideh-Fard (2001), the drop does not se em to recoil after it has reached its maximum radius. In the present simulation, the drop recoils, as also observed in the experimental photographs. Howe ver, the drop retracts too mu ch in the pres ent simulation compared to the experimental results as clearl y shown in the spreadi ng coefficient plot of Figure 6.18. This is likely due to the fact that the receding contact angle in the present computation is taken as 60 instead as 40. A lower contact angle suggests a higher wettability. Pasandideh-Fard et al. (2001) have not reported the spr eading coefficient of their numerical simulation. The initial spreadin g and beginning of the recoil is in very good agreement. At later time, differences between computed and measured spreading coefficients appear again likely due to th e receding difference in the receding contact angle. To offer additional insight of the present numerical simulation, the contact angle is plotted versus the dimensionles s contact line velocity and ve rsus dimensionless time in Figure 6.19. Note that the contact line velo city exceeds the impact velocity during the initial spreading and that the contact line velocity never falls below -20% of the initial
157 impact velocity. The contact angle during the course of the computa tion varies depending on the contact line velocity. This variation is illustrated in Figure 6.19. Compared to the plot of Figure 6.16, which is for a case for Re=300, the higher Re in the present case induces contact angle oscillations. Figure 6.17: Comparison of the instantaneou s drop shapes of the present numerical results with the numerical and experimental re sults of Pasandideh-Fard et al. (2001) of a water drop impinging on stainless steel. First and second column: numerical and experimental results of Pasandideh-Fard et al. (2001). Reprinted with permission from Elsevier Science. Third column: present numerical simulation. t = 7.3 ms t = 5.0 ms t = 3.1 ms t = 1.2 ms t = 0.38 ms t = 0.7 s (mm)
158 0 5 10 15 20 0 0.5 1 1.5 2 2.5 3 t*=tU/d0spreading coefficient Exp. Pasandideh-Fard et al. (2001) Present numerical simulation Re=2900 We=47 Figure 6.18: Comparison of the spreading coe fficient of present numerical results with the experimental results of Pasandideh-Fard et al. (2001). Figure 6.19: Contact angle versus dimensionl ess contact line velocity and versus time during the course of the pr esent computation for Re=2900, We=47 case of PasandidehFard et al. (2001). t*=tU/d0 0 5 10 15 20 50 60 70 80 90 100 110 contact angle Re=2900 We=47 dimensionless contact line velocity -0.5 0 0.5 1 1.5 2 50 60 70 80 90 100 110 Re=2900 We=47 contact angle
159 18.104.22.168 Case 2: impact of a water drop on a wax coated plate The computed and experimental spreading co efficients and maximum heights for Case 2 are plotted versus time in Figure 6.20. This case has been studied experimentally and numerically by Fukai et al. (1995). From the plots of Figure 6.20, the present numerical solution is in very good agreement with the re sults of Fukai et al. (1995) up to t*=20. After that, the present numerical simulation pr edicts a thinner (more elongated) drop than the results of Fukai et al. (1995) eventually reaching a point where a drop is formed and is ready to break up. That is why, present simulation results are not reported after t*=34. The contact angle during the cour se of the computation variati on is illustrated in Figure 6.21.
160 (a) (b) Figure 6.20: Comparison of the present numer ical results with th e experimental and numerical results of Fukai et al. (1995). (a) Spreading coeffici ent and (b) maximum height. t*=tU/r0 0 10 20 30 40 50 60 0 1 2 3 4 5 6 maximum height Exp. Fukai et al. (1995) Num. Fukai et al. (1995) Present numerical simulation Re=6200 We=128 t*=tU/r0 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 3.5 4 spreading coefficient Exp. Fukai et al. (1995) Num. Fukai et al. (1995) Present numerical simulation Re=6200 We=128
161 Figure 6.21: Contact angle versus dimensionl ess contact line velocity and versus time during the course of the pres ent computation for Re=6200, We=128 case of Fukai et al. (1995). 6.4 Heat Transfer Analysis In this paragraph the heat transfer process occurring during the impingement of a drop on a constant wall temperature is investig ated. The boundary conditions for the energy equation are illustrated in Fi gure 6.1. The bottom wall of the computational domain, where the drop impacts, is kept at a constant temperature Tw=1.0, and the top and side walls are adiabatic. A static contact angle is imposed. Initially, both the liquid and gas phases are at T0=0. First, the heat transfer is an alyzed for the case with Re=300 and We=60. Then the effect of the Re and We on the wall heat flux and Nusselt number is investigated. 6.4.1 Heat Transfer for Re=300, We=60 and =60 Case The thermal parameters used in the computations are presented in Table 6.5 along with the corresponding Prandtl numbers. The flui d properties are listed in Table 6.2. The contact angle in the present computation is fixed to =60. A contact angle of 60 0 10 20 30 40 50 60 56 60 64 68 72 76 80 84 88 92 96 t*=tU/r0 contact angle Re=6200 We=128 -0.5 0 0.5 1 1.5 2 56 60 64 68 72 76 80 84 88 92 96 dimensionless contact line velocitycontact angle Re=6200 We=128
162 corresponds, for example, to the static equi librium contact angle of deionized water on silicon oxide (Kim and Chun, 2001). The case simu lated is representative of the impact of 20 m deionized water droplet on silicon oxide w ith an impact speed of 14.6 m/s. This case is representative of the type of drop im pact envisioned in the micro-scale cooling device. The gravity is neglected in the simu lation. The dimensionless parameters based on the liquid properties are Re=300, We=60, Fr=, and Pr=2. The grid contains 40 cells per drop diameter and the time step t is 10-3. The isotherms during the impingement of a droplet on a constant temperature flat surface are shown in Figure 6.22. They develop more quickly in the gas phase than in the liquid droplet. This can be expl ained by the fact that the thermal diffusivity of the gas is 100 times larger than the one in the liquid, a nd by the fact that the convection effect is dominating during the early stage of impact . Once the droplet has reached its maximum spreading position, the flow changes direc tion, and the liquid temperature starts to increase, as the droplet recoils in a slower motion than the initial spreading. With time, the temperature distribution tends to become linear as shown in Figure 6.22 at t=4.0. Table 6.5: Computational thermal properties. Properties Liquid Vapor Ratio Conductivity ( k ) 0.2 0.01 20 Specific heat (Cp) 2.0 1.0 2 Diffusivity () 3.33 x 10-43.33 x 10-210-2 Prandtl Number 2.0 1.0 2 In the study of drop dynamics for heat transf er applications, it is of great interest to know the wall heat flux. Here, the heat flux de fined by Fourierâ€™s law is calculated on the bottom wall. The wall heat flux is:
163 () T qtk x (6.4) The wall heat flux is plotted in Figure 6. 23 versus the radial position at different instants. From Figure 6.23, one notices that th e wall heat flux is concentrated around the liquid droplet, while that in the gas phase is negligible. High heat flux rates take place at the early stage of impingement, and decays with time, with â€œpeaksâ€ present next to the contact line area. The present wall heat fl ux results follow the same trend as the one obtained by Pasandideh-Fard et al. (2001) when considering the conjugate heat transfer problem in the solid and liquid phases. Theref ore, the present assumption of keeping the wall isothermal is not affecting much the wall heat flux. This is also because the thermal conductivity of a metallic solid is much larger than the one of a fluid. The effects of the Reynolds and Weber numbers, as well as the contact angle, on the wall heat flux distribution are discussed next.
164 t=2.0 t=3.0 t=4.0 t=1.4 t=1.0 t=0.6 t=0.1 t=0.2 t=0.4 Figure 6.22: Instantaneous droplet shapes and isothermals every T=0.1 during droplet impact on a constant temperature flat surface for Re=300, We=60, =60.
165 r q 0 0 . 1 0 . 2 0 . 3 0 . 4 0 . 5 0 1 0 2 0 3 0 4 0 5 0 6 0 7 0 8 0 t = 0 . 0 1 t = 0 . 1 t = 0 . 2 t = 0 . 3 t = 0 . 4 t = 0 . 5 t = 0 . 6 Figure 6.23: Wall heat flux (q) versus radial position ( r ) at different time instants during the impact of a droplet on a constant te mperature flat surface for Re=300, We=60, =60. 6.4.2 Parametric Variation Effects on the Heat Transfer The same parameters variations as in th e study of drop dynamics are considered here. The isotherms for the reference case Re=100 We=20 =60are plotted in Figure 6.24, and the wall heat flux distribution along the ra dial direction at different time and the corresponding drop shapes are shown in Figure 6.25. Similar observations can be made from Figure 6.24 and Figure 6.25 as the one made from the previous case.
166 t=0.1 t=0.2 t=0.3 t=0.4 t=0.7 t=1.0 Figure 6.24: Instantaneous droplet shapes and isothermals every T=0.1 during droplet impact on a constant temperature flat surface for Re=100, We=20, =60.
167 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0 0.1 0.2 0.3 0.4 0.5 0 0 0.1 0.2 0.3 0.4 0.5 0 0 0.1 0.2 0.3 0.4 0.5 0 0 0.1 0.2 0.3 0.4 0.5 0 0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.1 r q 0 0.1 0.2 0.3 0.4 0.5 0 10 20 30 40t=0.1 t=0.2 t=0.3 t=0.4 t=0.7 t=1.0 t=0.1 t=0.2 t=0.3 t=0.7 t=0.4 t=1.0 Figure 6.25: Drop shapes and corresponding wall heat flux (q) versus radial position ( r ) at different time instan ts for Re=100, We=20, =60. The wall heat flux distributions are plotte d in Figure 6.26, Figure 6.27, Figure 6.28, and Figure 6.29 for varying contact angle, imp act velocity, droplet viscosity, and surface tension, respectively. From these figures, one can deduce that the wall heat flux evolution
168 strongly depends on the droplet deformations . As shown in Figure 6.26, Figure 6.27, Figure 6.28, and Figure 6.29 the wall heat fl ux is distributed under a wider area for lower contact angle, larger We and larger Re, sin ce the drop spreads more in those cases. To help clarify the parametric variation effects, the Nusselt number is evaluated based on the following definition: () ()2() () 2Rt o Rt oNtrdr Nut rdr (6.5) where R(t) is the basal drop radius which is the distance between the center and the contact point, and N(t) is defined as: () () Tqt Nt x k (6.6) The Nusselt number versus time is plotted in Figure 6.30 for the different cases. Note that the scale used for the Nusselt number is logarithmic. From Figure 6.30, the Nusselt number at t=0.1 is about 10 for all the cas es, except for the case with larger impact velocity where Nu is about 20. From Figure 6. 30 (b), on notices that Nu decreases faster implying that the heat dissipates faster for larger impact velocity. Some oscillations are present in the plots of Figure 6.30 (a), (c) and (d), which is in accordance with the recoiling oscillations observed in the sp reading coefficient pl ots of Figure 6.8.
169 r q 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 30 35 40t=0.1 t=0.2 t=0.3 t=0.4Re=100We=20t=0.7 t=1.0 r q 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 30 35 40t=0.1 t=0.2 t=0.3 t=0.4Re=100We=20t=0.7 t=1.0 (a) =60.0 (b) =80.0 r q 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 30 35 40t=0.1 t=0.2 t=0.3 t=0.4Re=100We=20t=0.7 t=1.0 (c) =100.0 Figure 6.26: Effect of the cont act angle on wall heat flux. (a) =60.0, (b) =80.0, and (c) =100.0.
170 r q 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 30 35t=0.1 t=0.2 t=0.3 t=0.4 t=0.7 t=1.0Re=100We=20 r q 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 30 35t=0.1 t=0.2 t=0.3 t=0.4 t=1.0Re=200We=80t=0.7 (a) U=1.0 (b) U=2.0 r q 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 30 35t=0.1 t=0.2 t=0.3 t=0.4Re=300We=180t=0.7 (c) U=3.0 Figure 6.27: Effect of the impact velocity on wa ll heat flux. (a) U=1.0, (b) U=2.0, and (c) U=3.0.
171 r q 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 30 35t=0.1 t=0.2 t=0.3 t=0.4 t=0.7 t=1.0Re=100We=20 r q 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 30 35t=0.1 t=0.2 t=0.3 t=0.4Re=200We=20t=0.7 t=1.0 (a) =0.2 (b) =0.1 r q 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 30 35t=0.1 t=0.2 t=0.3 t=0.4Re=400We=20t=0.7 t=1.0 (c) =0.05 Figure 6.28: Effect of the droplet viscosity on wall heat flux. (a) =0.2, (b) =0.1, and (c) =0.05.
172 r q 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 30 35t=0.1 t=0.2 t=0.3 t=0.4 t=1.0Re=100We=2 r q 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 30 35t=0.1 t=0.2 t=0.3 t=0.4 t=0.7 t=1.0Re=100We=20 (a) =10.0 (b) =1.0 r q 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 30 35t=0.1 t=0.2 t=0.3 t=0.4 t=0.7 t=1.0Re=100We=200 (c) =0.1 Figure 6.29: Effect of the surface tension on wall heat flux. (a) =10.0, (b) =1.0 and (c) =0.1.
173 0 0.5 1 1.5 2 2.5 3 10-1 100 101 102 timeNu 60 80 100 Re=100 We=20 0 0.5 1 1.5 2 2.5 3 10-1 100 101 102 timeNu U=1 Re=100 We=20 U=2 Re=200 We=80 U=3 Re=300 We=180 (a) contact angle effect (b) impact velocity effect 0 0.5 1 1.5 2 2.5 3 10-1 100 101 102 timeNu Re=100 Re=200 Re=400 We=20 0 0.5 1 1.5 2 2.5 3 10-1 100 101 102 timeNu We=2 We=20 We=200 Re=100 (c) droplet viscosity effect (d) surface tension effect Figure 6.30: Parametric effects on the Nusselt number. Effect of (a) contact angle, (b) impact velocity, (c) droplet viscosity, and (d) surface tension.
174 CHAPTER 7 7 SUMMARY AND FUTURE WORK 7.1 Summary and Conclusions An original micro-scale thermal manageme nt device concept has been proposed. The concept also motivates the present effort in developing a multiphase simulation capability for drop dynamics with heat transfer. The micro-scale cooling device takes adva ntage of two efficient modes of heat transfer, namely, dropwise condensation and drop/spray cooling. Thermodynamic and simplified analyses have been performed to i nvestigate the feasibility and performance of the device. The proposed conceptual micro-de vice contributes to the development of an efficient cooling device utilizi ng dropwise phase change. Becau se of its small scale, the device can operate largely inde pendent of the gravity level, making it also attractive to thermal management of space vehicles. The present micro-scale cooling device has potential in contributing to a variety of thermal manageme nt applications related to computers, electronics, avionics and bio-medi cal applications. In addition, due to the deviceâ€™s modularity, multiple units can be assembled to provide cooling according to the individual mission requirement. It is expect ed that with the progress made in microdevice fabrications and materials sciences, such a system can be manufactured in the foreseeable future. The key thermo-fluid mechanisms involve d within the device are drop dynamics related to growth by condensation, ejecti on by actuated membrane, impingement with heat transfer and evaporati on on a solid surface. These physical mechanisms are reviewed
175 in a comprehensive manner with direct experi mental and numerical visualization, as well as quantitative information. The numerical method adopted in the present work is based on the immersed boundary technique (Peskin, 1977; Udaykumar et al. 1997). A single set of conservation equations is solved on fixed Cartesian grid. The interf ace is represented by marker points. The jump in properties is accommodated via a smoothed profile that follows the distribution of a Heaviside function over four computational cells. The inte rfacial boundary conditions are embedded in the flow equations via the c ontinuum surface force (CSF) model (Brackbill et al., 1992). The interface is advected in a Lagrangian way by the marker velocity calculated by interpolation of the flow veloc ity, via numerical approx imation of the delta function. The single se t of flow equations are solved us ing a second accurate projection method (Ye et al., 1999) along with the multig rid technique (Udaykumar et al., 2001) to accelerate the rate of convergence of the pressure Poisson e quation. It is demonstrated that the single-phase flow solver is second order accurate, consistent with the results of Ye et al. (1999). The static e quilibrium of a spherical drop is computed to evaluate the performance of the multiphase flow solver . A grid refinement study demonstrates the multiphase flow solver is first order accurate, which is consistent with the immersed boundary formulation. The effect of the densit y ratio on the error is shown to increase slightly as the density ratio is increased. To further assess the performance of th e present smeared in terface treatment for multiphase fluid flows with large property ju mps, the dynamics of isothermal drop and bubble in an enclosed cylindrical domain is investigated. Specificall y, results for rising bubble by buoyancy previously reported Ryskin and Leal (1989) and Ye et al. (2001),
176 using, respectively, moving grid and fixedgrid, sharp-interface t echniques are compared. The present solutions are in general agreemen t with those reported. Furthermore, results for a rising oil drop in water are also obtained, and found in agreement with the experimental results of Klee and Treyba l (1956) and the three-dimensional VOF computations of Lrstad et al . (2002). The motion of water a nd R-12 drops relevant to the micro-device condition are then simulated. It is found that the drops dynamics, with the assigned Weber and capillary numbers, remain largely spherical because of the dominance of the surfac e tension effect over the inertia effect. For the modeling of the impingement of a si ngle drop on a flat surface, the contact angle is prescribed based on th e values reported in the litera ture. Both static and dynamic contact angle models are adopted , and their impact contrasted. It is observed that there is small difference between the tw o models during the early st age of impact; however, the difference becomes substantial during th e recoiling process. The dynamic model dampens the recoiling characteristics of the drop. Computations of millimeter drop size are also performed and found to agree with pr evious experimental and numerical results (Fukai et al., 1995; Pasandideh-Fard et al., 2001). Computations for a 20 m water droplet, relevant to the microdevice, are then conducted. The effects of inertia, surface tension, liquid wettability and viscosity on the drop dynamics and heat transfer are investigated. The dominance of the inertia effect during the early stage is clearly demonstrated, which is also consistent with the previously published work (Chandra and Avedisian, 1991). The wall heat flux and Nusse lt number are evaluated and are found to agree qualitatively with the results of Pasandi deh-Fard et al. (2001). In particular, high heat transfer rate occurs during the early st age of impingement. In all stages, higher wall
177 heat fluxes are identified to be close to the contact line. The present findings will directly contribute to the design of the micro-scale cooling device. The present research has contributed to the understanding of the basic physics involved in drop dynamics with heat transfer, as well as further development of the immersed boundary method for problems invol ving drop impingement and large property variations. 7.2 Future Work Several aspects of the comput ational framework can be improved to further extend its performance: â€¢ The spatial resolution in the immersed boundary technique can be improved using adaptive mesh refinement techniques, es pecially around the interface and contact point. This aspect is discussed in chapter 4. â€¢ Parallelization algorithms based on domain decomposition will help improve the computational speed. In order to achieve satisfactory parallel efficiency, and to address the interfacial dynamics and ma rker movement, communications between domains need to be synchronized. The issue of large property jump and the resulting stiffness should also be carefully treated. â€¢ Algorithms to treat break-up and coalescen ce should be added in order to treat the topological change. At the continuum level, such topological changes result in singularities. Nevertheless, computationa lly, approximate treatments can be developed. â€¢ Three-dimensional modeling should be employed to study non-axisymmetric phenomena, such as splashing of the drop upon impact. To perform full simulations of the thermo-fluid dynamics involved in the present micro-scale device, the following physical mechanisms need to be considered: â€¢ Conjugate heat transfer between the solid, liquid and vapor phases. â€¢ Condensation and evaporation of the drop. â€¢ Drop ejection from an actuated membrane with break-up.
178 LIST OF REFERENCES Abu-Orabi M., Modeling of heat transfer in dropwise condensation, International Journal of Heat and Mass Transfer, 1998, 41: 81-87. Adams R.L., and Roy J., A one-dimensional numerical model of a drop on demand inkjet, Journal of Applied Mechanics, 1986, 53:193-197. Aftomis M.J., Berger M.J. and Adomavic ius G.A., Parallel multilevel method for adaptively refined Cartesian grids with embedded boundaries, AIAA 38th Aerospace Sciences Meeting and Exhibit, AIAA 2000-0808, 2000. Agresar G., Linderman J.J., Tryggvason G. an d Powell K.G., An adaptive, Cartesian, front tracking method for the motion, deform ation and adhesion of circulating cells, Journal of Computational Physics, 1998, 143:346-380. Almgren A.S., Bell J.B., Collela P., Howell L.H. and Welcome M.L., A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations, Journal of Computa tional Physics, 1998, 142:1-46. Ashraf N.S., Carter III H.C., Casey L.C., C how L.C., Corban S., Drost M.K.,Gumm A.J., Hao Z., Hasan A.Q., Kapat J.S., Kramer L., Newton M., Sundaram K.B., Vaidya J., Wong C.C. and Yerkes K., Design and an alysis of a meso-scale refrigerator, Department of Mechanical, Materials and Aerospace Engineering, University of Central Florida, FL 32816, USA, 1999. Aydin O. and Yang W.-J., On microtransport phenomena in minute droplets: a critical review, Applied Mechanic s Review, 2000, 53:101-116. Badea A., Carasso C. and Panasenko G., A model of a homogenized cavity corresponding to a multinozzle droplet gene rator for continuous ink-jet printers, Numerical Methods for Partial Differential Equations, 1998, 14:821-842. Badie R. and Frits de Lange D., Mechanis m of drop constriction in a drop-on-demand ink-jet system, Proceedings of the Royal Society of London, Series A-Mathematical, Physical and Engineering Sciences, 1997, 453: 2573-2581. Batchelor G.K., An introduction to fluid dynamics, Cambridge University Press, Cambridge, 1967. Berger M.J., Adaptive mesh refinement for hyperbolic partial-differential equations, Journal of Computational Physics, 1984, 53:484-512.
179 Bernardin J.D., Stebbins C.J. and Mudawar I., Mapping of impact and heat transfer regimes of water drops impinging on a polishe d surface, International Journal of Heat and Mass Transfer, 1997, 40(2):247-267. Bhaga D. and Weber M.E., Bubbles in visc ous liquids: shapes, wakes and velocities, Journal of Fluid Mechanics, 1981, 105:61-85. Bishop D., Gammel P. and Giles C.R., The littl e machines that are making it big, Physics Today, October 2001. Bolle L. and Moureau J.C., Spray cooling of hot surfaces, in Multiphase Science and Technology, edited by Hewitt G.F., Delhaye J.M. and Zuber N., Hemisphere, 1976, 1:1-97. Bonacina C., Del Giudice S. and Comini G ., Dropwise evaporation, Journal of Heat Transfer, 1979, 101:441-446. Brackbill J.U., Kothe D.B. and Zemach C ., A continuum method for modeling surface tension, Journal of Computat ional Physics, 1992, 100:335-354. Brackbill J.U., Juric D., Torres D. and Ka llman E., Dynamic modeling of microgravity flow, Proceedings of the fourth fluid physics and transport phenomena conference, NASA, August 12-14, 1998, Cleveland, Ohio. Burnside B.M. and Hadi H.A., Digital co mputer simulation of dropwise condensation from equilibrium droplet to detectable size, International Journal of Heat and Mass Transfer, 1999, 42:3137-3146. Burton T.M. and Eaton J.K., Analysis of a fractional-step method on overset grids, Journal of Computational Physics, 2002, 177:336-364. Bussmann M., Chandra S. and Mostaghimi J., Modeling the splash of a droplet impacting a solid surface, Physics of Fluids, 2000, 12:3121-3132. Bussmann M., Mostaghimi J. and Chandra S ., On a three-dimensional volume tracking model of droplet impact, Phys ics of Fluids, 1999, 11:1406-1417. Carey V.P., Liquid-vapor phase change phenom ena, an introduction of vaporization and condensation processes in heat transfer equipment, Taylor and Francis, Series in Chemical and Mechanical Engineering, 1992. Carey V.P., Oyumi S.M. and Ahmed S., Post -nucleation growth of water microdroplets in supersaturated gas mixtures: a molecula r simulation study, International Journal of Heat and Mass Transfer, 1997, 40:2393-2406. Cengel Y.A. and Boles M.A., Thermodynamics: An Engineering Approach, Third Edition, McGraw Hill, Boston, 1998.
180 Chandra S. and Avedisian C.T., On the coll ision of a droplet w ith a solid surface, Proceedings of the Royal Society, London, Serie A, 1991, 432:13-41. Chen L., Garimella S.V., Reizes J.A. a nd Leonardi E., The development of a bubble rising in a viscous liquid, Journal of Fluid Mechanics, 1999, 387:61-96. Chen P.H., Chen W.C. and Chang S.H., Bubbl e growth and ink jet ejection process of a thermal ink-jet printhead, In ternational Journal of Mechanical Science, 1997, 39:683695. Chilton T.H. and Colburn A.P., Mass transfer co efficients-prediction data on heat transfer fluid motion, Industrial Engin eering Chemistry, 1934, 26:1183-1187. Chorin A.J., Numerical solution of the Navier Stokes equations, Mathematics of Computation, 1968, 22:745-762. Clift R., Grace J.R. and Weber M.E., Bubbles , Drops, and Particles, Academic Press, New York, 1978. Collier J.G. and Thome J.R., Convective Bo iling and Condensation, Third Edition, Oxford University Press, New York, 1994. Crank J., Free and Moving Boundary Problems, Oxford University Press, New York, 1984. Dandy D. S. and Leal G. L., Buoyancy-dri ven motion of a deformable drop through a quiescent liquid at intermediate Reynolds numbers, Journal of Fluid Mechanics, 1989, 208: 161-192. Delhaye J.M., Jump conditions and entropy sour ces in two-phase systems: local instant formulation, International Journa l of Multiphase Flow, 1974, 1:395-409. De Vahl Davis G. and Jones I.P., Natural c onvection in a square cavity: a comparison exercise, International Journal for Nume rical Methods in Fluids, 1983, 3:227-248. Dijksman J.F., Hydro-acoustics of piezoelectr ically driven ink-jet print heads, Flow Turbulence and Combustion, 1999, 61:211-237. Dussan V. E.B., On the spreading of liquids on solid surfaces: static and dynamic contact lines, Annual Review of Fluid Mechanics, 1979, 11:371-400. Eucken A., Energie und Stoffaustausch an Grenzflchen, Naturwissenschaften, 1937, 25:209-218. Ferziger J.H. and Peric M., Computational methods for fluid dynamics, Second Edition, Springer, Berlin, 1996.
181 Floryan J.M. and Rasmussen H., Numerical methods for viscous flows with moving boundaries, Applied Mechanics Reviews, 1989, 42:323-340. Francois M., A study of the volume-of-f luid method for moving boundary problems, Master of Science Thesis, Embry-Ri ddle Aeronautical University, 1998. Fromm J.E., Numerical calcula tion of the fluid dynamics of drop-on-demand jets, IBM Journal of Research and De velopment, 1984, 28:96-110. Fukai J., Shiiba Y. and Miyatake O., Theoretic al study of droplet impingement on a solid surface below the Leidenfrost temperature, International Journa l of Heat and Mass Transfer, 1997, 40:2490-2492. Fukai J., Shiiba Y., Yamamoto T. and Miyata ke O., Wetting effects on the spreading of a liquid droplet colliding with a flat surf ace: Experiment and modeling, Physics of Fluids, 1995, 7:236-247. Fukai J., Zhao Z., Poulikakos D., Megaridi s C.M. and Miyatake O., Modeling of the deformation of a liquid droplet impinging upon a flat surface, Physics of Fluids, 1993, 5:2588-2599. Fujikawa S., Kotani M. and Sato H., Mol ecular study on phase transition phenomena of fluid, Thermal Science and Engineering, 1995, 3:45-50. Ghia U., Ghia K.N. and Shin C.T., High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid met hod, Journal of Computational Physics, 1982, 48:387-411. Ghodbane M. and Holman J.P., Experimental study of spray cooling with Freon-113, International Journal of Heat and Mass Transfer, 1991, 34:1163-1174. Gose E.E., Mucciardi A.N. and Baer E., Model for dropwise condensation on randomly distributed sites, Internat ional Journal of Heat and Mass Transfer, 1967, 10:15-22. Grassia, P., The design of droplet ejector, Journal of Fluids Engineering, 1999, 121:658664. Griffith P., Dropwise condensation, Handbook of Heat Transfer Fundamentals, edited by Rohsenow W.M., Hartnett J.P., Ganic E. N., McGraw-Hill, New-York, 1985, Chapter 11, 37-50. Gueyffier D., Nadim A., Scardovelli R. a nd Zaleski S., Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows, Journal of Computational Physics, 1999, 152:423-456. Guo Y. and Lee T.-Y.T., Experimental char acterizations of dyna mic properties of a silicon membrane structure, ASME Paper, 99-IMECE/EEP-21, 1999.
182 Han J. and Trygggvason G., Secondary brea kup of axisymmetric liquid drops. I. Acceleration by a constant force, Physics of Fluids, 1999, 11: 3650-3667. Harlow F.H. and Shannon J.P., The splash of a liquid drop, Journal of Applied Physics, 1967, 38:3855-3866. Harlow F.H. and Welch J.E., Numerical calculation of time dependent viscous incompressible flow of fluid with free surface, Physics of Fluids, 1965, 8:2182-2189. Harmathy T.Z., Velocity of large drops and bubbl es in media of infinite or restricted extent, American Institue of Chemi cal Engineers Journal, 6:281-288. Hatta N., Fujimoto H. and Takuda H., Deform ation process of a water droplet impinging on a solid surface, Journal of Heat Transfer, 1995, 117:394-401. Heffington S.N., Black W.Z. and Glezer A., Vi bration-Induced Droplet Atomization Heat Transfer Cell for High-Heat Flux Dissipa tion. Thermal Challenges in Next Generation Electronic Systems (THE RMES-2002), 2002, Santa Fe, NM. Hirt C.W. and Nichols B.D., Volume of Fluid (VOF) Method for the dynamics of free boundaries, Journal of Computati onal Physics, 1981, 39:201-225. Hnat J.G. and Buckmaster J.D., Spherical cap bubbles and skirt formation, The Physics of Fluids, 1976, 19:182-194. Ho C.-M. and Tai Y.-C., Micro-Electro-Mech anical Systems (MEMS) and fluid flows, Annual Review of Fluid Mechanics, 1998, 30:579-612. Holman J.P. and Kendall C.M., Extended studies of spray cooling with Freon-113, International Journal of Heat and Mass Transfer, 1993, 36:2239-2241. Ishii M., Thermo-Fluid Dynamic Theory of Two-Phase Flow, Eyrolles, Paris, 1975. Jakob M., Heat transfer in evaporation a nd condensation-II., Mech. Eng., 1936, 58:729. Jerman H., Electrically-activated, normally-closed diaphragm valves, Journal of Micromechanics and Microengineering, 1994, 4:210-216. Juric D. and Trygggvason G., Computations of boiling flows, International Journal of Multiphase Flow, 1998, 24:387-410. Kan H.-C., Udaykumar H.S., Shyy W., and Tran-Son-Tay R, Hydrodynamics of a compound drop with application to Leukocyt e modeling, Physics of Fluids, 1998, 10:760-774. Kaviany M., Principles of Convective Heat Transfer, Springer-Verlag, New-York, 1994. Kim H.-Y. and Chun J.-H., The recoiling of liquid droplets upon collision with solid surface, Physics of Fluids, 2001, 13:643-659.
183 Kim J. and Moin P., Application of a fract ional step method to incompressible NavierStokes equations, Journal of Com putational Physics, 1985, 59:308-323. Klee A.J. and Treybal R.E., Rate of rise or fa ll of liquids drops, A.I.Ch.E. Journal, 1956, 2:444-447. Kobayashi Y., Toyokawa S. and Araki T., Heat and mass transfer from thin liquid film in the vicinity of the interline of meniscus, 1994, Thermal Science and Engineering, 1994, 2 (1):45-51. Kothe D.B., Rider W. J., Mosso S. J., Brock J.S. and Hochstein J.I., Volume tracking of interfaces having surface te nsion in two and three dimensions, AIAA paper 96-0859, 1996. Kothe D.B. and Mjolsness R.C., RIPPLE: A new model for incompressible flows with free surfaces, AIAA Journal, 1992, 30:2694-2700. Lai C.L., Gibbs-Thomson effect on droplet conde nsation, Journal of Heat Transfer, 1999, 121:632-638. Lay J.H. and Dhir V.K., 1995, Shape of a va por stem during nucleat e boiling of saturated liquids, Journal of Heat Transfer, 1995, 117:394-401 Lee J., Kim J. and Kiger K.T., Timeand space -resolved heat transfer characteristics of single droplet cooling using microscale heater arrays, Internationa l Journal of Heat and Fluid Flow, 2001, 22:188-200, 2001. Le Fevre E.J. and Rose J.W., A theory of heat transfer by dropwise condensation, Proceedings of the Third International H eat Transfer Conference, Chicago, 1966, 2:362-375. Le Quere P., Accurate solutions to the squa re thermally driven cavity at high Rayleigh number, Computers Fluids, 1991, 20:29-41. Liu G.W., Morsi Y.S., and van der Walt J.P., An alysis of spray cooling heat flux, Journal of Heat Transfer, 1999, 121:742-745. Lfdahl L. and Gad-El-Hak M., MEMS applic ations in turbulence and flow control, Progress in Aerospace Sciences, 1999, 35:101-203. Lrstad D., Francois M., Shyy W. and Fuch s L., Immersed boundary and volume of fluid investigations of single rising droplets, 2002, submitted to AIAA. Majumdar A. and Mezic I., Instability of ultra-thin water films and the formation on hydrophilic surfaces, Journal of He at Transfer, 1999, 121:964-971. Marzo (di) M. and Evans D.D., Evaporation of a water droplet deposited on a hot high thermal conductivity surface, Journal of Heat Transfer, 1989, 111:210-213.
184 Marzo (di) M., Tartarini P. and Liao Y., Ev aporative cooling due to a gently deposited droplet, International Journal of He at and Mass Transfer, 1993, 36:4133-4139. Matsumoto M. and Fujikawa S., Nonequili brium vapor condensation: molecular simulation and shock-tube experiment, Microscale Thermophysical Engineering, 1997, 1:119-126. Matsumoto M., Yasuoka K. and Kataoka Y., Molecular mechanism of evaporation and condensation, Thermal Science and Engineering, 1995, 3:27-31. McCormick J.L. and Westwater J.W., Nucl eation sites for dropwise condensation, Chemical Engineering Science, 1965, 20:1021. Melton J.E., Automated three-dimensional Ca rtesian grid generation and Euler flow solutions for arbitrary geometries, Ph.D. Thes is, University of California, Department of Mechanical Engineering, Davis CA, 1996. Middleman S., Modeling Axisymmetric Flows, Dynamics of Films, Jets, and Drops, Academic Press, San Diego, 1995. Nagai N. and Carey V.P., Assessment of surf ace wettability and its relation to boiling phenomena, Proceedings of the 2001 ASME International Mechanical Engineering Congress and Exposition, HTD-24132, N ovember 11-16, 2001, New York, NY, USA. Osher S. and Sethian J.A., Fronts propa gating with curvature dependent speed: Algorithms based in Hamilton-Jacobi fo rmulations, Journal of Computational Physics, 1988, 79:12-49. Pasandideh-Fard M., Aziz S.D., Chandra S. an d Mostaghimi J., Cooling Effectiveness of a water drop impinging on a hot surface, International Journal of Heat and Fluid Flow, 2001, 222:201-210. Pasandideh-Fard M., Qiao Y.M., Chandra S. and Mostaghimi J., Capillary effects during droplet impact on a solid surface, Physics of Fluids, 1996, 8:650-659. Patankar S.V., Numerical heat tr ansfer and fluid flow, Series in computational methods in mechanics and thermal sciences, Taylor and Francis, Washington DC, 1980. Peeters E. and Verdonckt-Vandebroek S., Th ermal ink-jet technology, IEEE Circuits and Devices, 1997, 13: 19-23. Percin G., Levin, L. and Khuri-Yakub, B.T ., Piezoelectrically actuated droplet ejector, Review of Scientific Instruments, 1997, 68:4561-4563. Peregrine D.H., Shoker G. and Simon, A., The bifurcation of liquid bridges, Journal of Fluid Mechanics, 1990, 212:25.
185 Peskin C.S., Numerical analysis of blood flow in the heart, Jour nal of Computational Physics, 1977, 25:220-252. Peterson G.P., Swanson L.W. and Gerner F.M., Micro heat pipes, in Microscale Energy Transport, Series in Chemical and Mechanic al Engineering, edited by Tien C.-L., Majumdar A., Gerner F.M., Ta ylor and Francis, 1997, 295-337. Peterson P.F., Bai R.Y., Schrock V.E. and H ijikata K., Droplet condensation in rapidly decaying pressure fields, Journal of Heat Transfer, 1992, 114:194-200. Press W.H., Teukolsky S.A., Vetterling W.T. , Flannery B.P., Numerical Recipes in FORTRAN the art of scientif ic computing, Second Edition, Cambridge University press, Cambridge, 1992. Rembe C., aus der Wiesche S. and Hofer E.P., Thermal ink jet dynamics: modeling, simulation and testing, Microelect ronics Reliability, 2000, 40:525-532. Rieber M. and Frohn A., A nume rical study on the mechanism of splashing, International Journal of Heat and Fluid Flow, 1999, 20:455-461. Rizza J.J., A numerical solution to dropwise evap oration, Journal of Heat Transfer, 1981, 103:501-507. Rose J.W., Utaka Y. and Tanasawa I., Dropwise condensation, in Handbook of Phase Change: Boiling and Condensation, edited by Kandlikar S., Shoji M., Dhir V.K., Taylor and Francis, Philadelphia, 581-594, 1999. Ryskin G. and Leal L.G., Numerical so lution of free-boundary problems in fluid mechanics. Part 1, Part 2, Part 3, Journal of Fluid Mechanics, 1984, 148:1-43. Sadhal S.S., Ayyaswami P.S. and Chung J. N., Transport phenomena with drops and bubbles, Springer, New York, 1997. Sadhal S.S. and Plesset M.S., Effect of so lid properties and contact angle in dropwise condensation and evaporation, Journa l of Heat Transfer, 1979, 101:48-54. Saif, M.T.A., Alaca, B. E., Sehitoglu, H., An alytical modeling of Electrostatic Membrane Actuator for Micro Pumps, IEEE Journal of Microelectromechan ical systems, 1999, 8:335-345. Sarno C. and Moulin G., Thermal management of highly integrated electronic packages in avionics applications, Electronics Cooling Magazine, November 2001, 7(4), http://www.electronics-cooling.com/html/2001_nov_a1.html , accessed on 8/2/02. Scardovelli R. and Zaleski S., Direct numerical simulation of free-surface and interfacial flow, Annual Review of Fluid Mechanics, 1999, 31:567-603.
186 Semiconductor Industry Association (SIA), International Technology Roadmap for Semiconductors: 1999 edition, Austin, Texas: International SEMATECH. Senturia S.D., Microsyste m Design, Kluwer Academic Publishers, Boston, 2001. Sethian J.A., Level set methods: evolving in terfaces in geometry, fluid mechanics, computer vision, and materials science, Cambridge University Press, Cambridge, 1996. Shannon M.A., Philpott M.L., Miller N.R., Bullard C.W., Beebe D.J., Jacobi A.M., Hrnjak P.S., Saif T., Aluru N., Sehitoglu H., Rockett A. and Economy J., Integrated mesoscopic cooler circuits (IMCCS), AES-Vol. 39, Proceedings of the ASME Advanced Energy Systems Division, 1999. Shield T.W., Bogy D.B. and Talke F.E., A numerical comparison of one-dimensional fluid jet models applied to drop-on-dema nd printing, Journal of Computational Physics, 1986, 67:327-347. Shield T.W., Bogy D.B. and Talke F.E., Drop formation by DOD ink-jet nozzles: a comparison of experiment and numerical si mulation, IBM Journal of Research and Development, 1987, 31:96-110. Shyy W., Computational modeling for fluid fl ow and interfacial Transport, Elsevier, Amsterdam, 1994. Shyy W., Francois M., Udaykumar H.S., Nâ€™Dri N. and Tran-Son-Tay R., Moving boundaries in micro-scale bi ofluid dynamics, Applied Mechanics Reviews, 2001, 54:405-453. Shyy W., Udaykumar H.S., Rao M.M. and Sm ith R.W., Computational Fluid Dynamics with Moving Boundaries, Taylor and Francis, Washington DC, 1996. Smith M.K., Vukasinovic B. and Glezer A., Vibration-induced droplet atomization, Proceedings of the fourth fluid physics a nd transport phenomena conference, NASA, August 12-14, 1998, Cleveland, Ohio. Son G.H., A numerical method for bubble moti on with phase change, Numerical Heat Transfer Part B, 2001, 39(5):509-523. Son G., Dhir V.K. and Ramanujapu N., Dynami cs and heat transfer associated with a single bubble during nucleate boiling on a horizontal surface, Journal of Heat Transfer, 1999, 121:623-632. Sussman M. and Smereka P., Axisymmetric free boundary problems, Journal of Fluid Mechanics, 1997, 341: 269-294.
187 Sussman M., Smereka P. and Osher S., A leve l-set approach for computing solutions to incompressible two-phase flow, Journal of Computational Physics, 1994, 114: 146159. Szita N., Sutter R., Dual J. and Buser R., A micropipettor with integrated sensors, Sensors and Actuators A, 2001, 89:112-118. Tanaka H., A theoretical study of dropwise c ondensation, Journal of Heat Transfer, 1975, 97:72-78. Tanasawa I., Advances in condensation heat transfer, Advances in Heat Transfer, Academic Press, 1991, 21:55-140. Tien N.C., Silicon micromachined thermal sensors and actuators, in Microscale Energy Transport, Series in Chemical and Mechanic al Engineering, edited by Tien C.-L., Majumdar A., Gerner F.M., Ta ylor and Francis, 1997, 369-386. Tryggvason G., Bunner B., Esmaeli A., Juric D., Al-Rawahi N., Tauber W., Han J., Nas S. and Jan Y.-J., A front-tr acking method for the computations of multiphase flow, Journal of Computational Physics, 2001, 169:708-759. Tseng F.-G., Kim C.-J. and Ho C.-M., A Mi croinjector Free of Satellite Drops and Characterization of the Ejected Droplet s, ASME DSC-Vol.66, Micro-EelectroMechanical Sytems (MEMS), 1998. Udaykumar H.S., Kan H.-C., Shyy W., and Tran-Son-Tay R., Multiphase dynamics in arbitrary geometries on fixed Cartesian gr ids, Journal of Computational Physics, 1997, 137:366-405. Udaykumar H.S., Mittal R. and Shyy W., Computation of solid-liquid phase fronts in the sharp interface limit on fixed grids, Journal of Computa tional Physics, 1999, 153:535574. Umur A. and Griffith P., Mechanism of dropw ise condensation, Journal of Heat Transfer, 1965, 87:275-282. Unverdi S.O. and Trygggvason G., A front tr acking method for viscous incompressible flows, Journal of Computational Physics, 1992, 100:25-37. Versteeg H.K. and Malalasekera W., An in troduction to computa tional fluid dynamics. The Finite Volume Method, Longman, 1995. Vincenti W.G. and Krueger C.H., Introduc tion to physical gas dynamics, Krieger Publishing Company, Malabar, Florida, 1975. Vukasinovic B., Glezer A. and Smith M.K., Vibration induced droplet atomization, Physics of Fluids, 2000, 12, S12.
188 Vukasinovic B., Heffington S.N., Smith M.K. and Glezer A., Vibration-induced Droplet Atomization (VIDA) for Two-phase Thermal Management, Proceedings of the ASME International Mechanical Engi neering Congress and Exposition (IMECE) 2001, Electronic and Photonic Packaging Division, November 11-16, 2001, New York, NY. Wan D.C., Patnaik B.S.V. and Wei G.W., A new benchmark quality solution for the boundary-driven cavity by discrete singular co nvolution, Numerical Heat Transfer, Part B, 2001, 40:199-228. Wayner Jr. P.C., Spreading of a liquid film with a finite contact angle by the evaporation/condensation proce ss, Langmuir, 1993, 9:294-299. Wayner Jr. P.C., Thermal and mechanical effects in the spreading of a liquid film due to a change in the apparent fin ite contact angle, Journal of Heat Transfer, 1994, 116:938945. Watchers L.H. and Westerling N.A.J., The he at transfer from a hot wall to impinging water drops in the spheroidal state, Chemical Engineering Science, 1966, 21:10471056. Weiss D.A. and Yarin A.L., 1999, Single drop im pact onto liquid films: neck distortion, jetting, tiny bubble entrainment, and crown formation, Journal of Fluid Mechanics, 385:229-254. Welch J.E., Harlow F.H., Shannon J.P. and Daly B.J., The MAC Method: a computational technique for solving viscous , incompressible, transient fluid-flow problems involving free surfaces, Los Alamos Scientific Laboratory, Report No. LA3425, 1966. Welch J.F. and Westwater J.W, Micros copic study of dropwise condensation, International Developments in Heat Transf er, Proceedings of the International Heat Transfer Conference, ASME, Part II, 1961, 302-309. Wilkes E.D., Phillips S.D. and Basaran A., Co mputational and experimental analysis of dynamics of drop formation, Physics of Fluids, 1999, 11:3577-3598. Wilkes E.D. and Basaran O.A., Drop ejection from an oscillating rod, Journal of Colloid and Interface Science, 2001, 242:180-201. Williams M.W., Kothe D.B. and Puckett E.G, Accuracy and convergence of continuum surface-tension models, Fluid Dynamics at Interfaces, Edited by W. Shyy and R. Narayanan, Cambridge University Press, 1999. Yang J., Chow L.C. and Pais M.R., Nucleate boiling heat transfer in spray cooling, Journal of Heat Transfer, 1996, 118:668-671.
189 Ye T., Mittal R., Udaykumar H.S. and Shyy W., An accurate Cartesian grid method for viscous incompressible flows with co mplex immersed boundaries, Journal of Computational Physics, 1999, 156:209-240. Ye T., Shyy W. and Chung J.N., A fixe d-grid, sharp-interface method for bubble dynamics and phase change, Journal of Computational Physics, 2001, 174:781-815. Zang Y., Street R.L. and Koseff J.R., A nonstaggered grid, fractional step method for time dependent incompressible Navier-Stoke s equations in Curvilinear coordinates, Journal of Computational Physics, 1994, 114:18-33. Zapalowicz Z., Analysis of the process of wetting a horizontal h eated surface area by a single water droplet, International Communi cations in Heat and Mass Transfer, 1995, 22:713-720. Zhao Z., Poulikakos D. and Fukai J., H eat transfer and flui d dynamics during the collision of a liquid droplet on a substrate I. Modeling, Internati onal Journal of Heat and Mass Transfer, 1996a, 39:2771-2789. Zhao Z., Poulikakos D. and Fukai J., H eat transfer and flui d dynamics during the collision of a liquid droplet on a substrate II. Experiments, International Journal of Heat and Mass Transfer, 1996b, 39:2791-2802. Zheng L.L. and Zhang, An adaptive leve l set method for moving-boundary problems: application to droplet spread ing and solidification, Numeri cal Heat Transfer, Part B, 2000, 37:437-454. Zhu X. and Kim E.S., Microfluidic motion ge neration with acoustic waves, Sensors and Actuators A, 1998, 66:355-360.
190 BIOGRAPHICAL SKETCH Marianne Franois was born in 1975 in Normandy, France. She received her Baccalaurat in mathematics and physics in 1992. Then, she attended Ecole Polytechnique Fminine in Sceaux, near Paris. In 1996, for her last year of study, she was selected to be an exchange student at Embr y-Riddle Aeronautical University in Daytona Beach, Florida, and started her Master of Sc ience studies in Aerospace Engineering. She received her Diplme dâ€™Ingnieur in 1997 a nd her Master of Science in Aerospace Engineering with distinction in 1998. In fa ll 1998, she entered the Ph.D. program at the University of Florida and joined the computational thermo-fluid group of the Aerospace Engineering, Mechanics and Engineering Scie nce Department. There, she has performed research in the area of computational flui d dynamics with moving boundaries applied to the fluid-structure interaction pr oblem of a soft contact lens with the tear film and bubble and drop dynamics for heat transfer enhancem ent. She is a 2000 Amelia Earhart Fellow.