UFDC Home  myUFDC Home  Help 



Full Text  
PAGE 1 NANOSCALE TRANSPORT MECHANISMS ASSOCIATED WITH THREEPHASE CONTACT LINES By SHALABH CHANDRA MAROO A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLOR IDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2009 1 PAGE 2 2009 Shalabh Chandra Maroo 2 PAGE 3 To my parents for their endless love and support 3 PAGE 4 ACKNOWLEDGMENTS I am greatly thankful to my advisor, Dr. J acob N. Chung, for all the time and efforts he has put in towards mentoring my dissertation wo rk. I am indebted to him for his invaluable guidance, support and encouragement. His dedicati on and hard work towards research will be an everlasting source of inspiration for me. I express my sincere gratitude to Dr. S. Balachandar, Dr. Youping Chen, Dr. HaiPing Cheng and Dr. James F. Klausner for serving as my supervisory committee members. Their precious advice and recommendations have considerably helped improve this work. My special acknowledgement goes to all my colleagues in my research group, especially, TaeSeok Lee, Dr. Yun Whan Na and Dr. Chengfeng Tai, with whom discussions have been very insightful and helpful. My utmost appreciation goes to my pare nts and family for their endless love and support which I have been blessed with in my life. 4 PAGE 5 TABLE OF CONTENTS page ACKNOWLEDGMENTS...............................................................................................................4 LIST OF TABLES................................................................................................................. ..........8 LIST OF FIGURES.........................................................................................................................9 ABSTRACT...................................................................................................................................14 CHAPTER 1 INTRODUCTION AND LITERATURE REVIEW..............................................................16 Introduction................................................................................................................... ..........16 Need for Molecular Dynamics Simulation.............................................................................18 Literature Review.............................................................................................................. .....20 Nucleate Boiling..............................................................................................................2 0 Critical Heat Flux............................................................................................................27 Droplet Impact and Spreading.........................................................................................28 Wall Models....................................................................................................................34 Ultrathin Films, Disjoining Pr essure and Hamaker Constant.........................................37 Capillary Pressures at Nanoscale.....................................................................................41 Research Objectives............................................................................................................ ....46 2 THEORY AND APPROACH................................................................................................50 Molecular Dynamics Simulation............................................................................................50 Force Calculation Algorithm..................................................................................................51 Brute Force Method.........................................................................................................52 Linked List Method.........................................................................................................52 Integrator Method...................................................................................................................53 NonDimensional Units..........................................................................................................53 Boundary Conditions..............................................................................................................55 Periodic Boundary Condition (PBC)...............................................................................55 Mirror Boundary Condition.............................................................................................55 Thermal Wall Boundary Condition.................................................................................56 Initial Structure.............................................................................................................. .........56 Simulation Stages.............................................................................................................. .....57 Statistical Sampling........................................................................................................... .....58 Density........................................................................................................................ .....58 Energy..............................................................................................................................58 Temperature.................................................................................................................... .59 Heat Flux.........................................................................................................................59 Thermal Conductivity......................................................................................................60 Pressure....................................................................................................................... .....60 5 PAGE 6 Pressure Verification..............................................................................................................61 Distinguishing Between Liquid At oms and Vapor Atoms for Argon....................................65 FluidWall Thermal Eq uilibrium Model................................................................................66 Model Validation.............................................................................................................70 Internal energy and thermal conductivity.................................................................70 Heat equation............................................................................................................73 Maxwell speed distribution......................................................................................82 Adiabatic Walls...............................................................................................................8 5 3 TRANSIENT FILM EVAPORATION..................................................................................88 Transient Film Evaporation....................................................................................................8 8 Heat Flux and Evaporation Rates....................................................................................96 Evaluation of Hamaker Constant from Molecular Dynamics Simulation.......................96 Parameters Governing the Forma tion of NonEvaporating Film.........................................100 Heat Flux, Evaporation Flux, Hamaker Constant for Nanochannel..............................105 Varying Nanochannel Height h .....................................................................................106 Varying Film Thickness tfilm..........................................................................................108 Vapor PressureDensity Plots........................................................................................111 Case 7 Simultaneous Evaporation and Condensation................................................113 4 NANODROPLET IMPACT ON A HOMOGENEOUS SURFACE..................................116 Evaluation of Drop Temperature..........................................................................................117 Evaluation of Interface Markers, In terfacial Fit and Contact Angle....................................117 Droplet Impingement Case 1.............................................................................................118 Droplet Impingement Case 2.............................................................................................121 Droplet Impingement Case 3 and Case 4..........................................................................123 5 MENISCUS EVAPORAT ION AT NANOSCALE.............................................................125 Simulation Model............................................................................................................... ..125 Results and Discussion Case I...........................................................................................133 Results and Discussion Case II..........................................................................................142 Meniscus Boundary Condition.............................................................................................148 Analytical Meniscus Boundary Condition....................................................................148 CurveFitted Meniscus Boundary Condition.................................................................150 6 CONCLUSION AND FUTURE WORK.............................................................................153 Conclusion............................................................................................................................153 Proposed Thermal Wall Model.....................................................................................153 Transient Film Evaporation...........................................................................................154 Simulation of Seven Cases in Nanochannel Film Evaporation.....................................154 NanoDroplet Impact on a Homogeneous Surface........................................................155 Meniscus Evaporation at Nanoscale..............................................................................156 Future Work..........................................................................................................................157 Molecular Dynamics Simulation (Nano Region)..........................................................157 6 PAGE 7 CoarseGrained Method (Micro Region)......................................................................159 Finite Volume Method (Macro Continuum Region).....................................................159 MultiScale Model and Matching Procedure................................................................160 MD simulation and MDCG matching (nanomicro matching).............................160 CG and FVM matching (micromacro matching)..................................................161 LIST OF REFERENCES.............................................................................................................163 BIOGRAPHICAL SKETCH.......................................................................................................171 7 PAGE 8 LIST OF TABLES Table page 21 Dimensionless Parameters.................................................................................................54 22 Simulation parameters fo r pressure verification................................................................62 23 Change in internal energy for validat ion of proposed wall heat transfer model................73 24 Thermal conductivity for validation of proposed wall heat transfer model.......................73 31 Parameters for seven cases simulated in nanochannel.....................................................101 32 Case 13 results for film evaporation in nanochannel.....................................................108 33 Case 2, 46 results for film evaporation in nanochannel.................................................109 8 PAGE 9 LIST OF FIGURES Figure page 11 Flowchart showing the basic struct ure of molecular dynamics simulation.......................17 12 Schematic of bubble growth, (a) overall picture at macroscale (shaded region depicts the area covered by nonevaporating thin film), and (b) zoomed in nanoand microscale regions at the threephase contact line...........................................................21 13 Single Bubble Ebullition Cycle.........................................................................................23 14 Expected variation of film thickness with film thickness..................................................37 15 Gradual and complete film evapor ation simulated by Yi et al. (2002)..............................41 16 Phase diagram of water (Angell, 1988).............................................................................45 17 Schematic of water phase diagram show ing a metastable cap illary bridge. The capillary bridge under tension acts like a wate r tube sealed with a piston, which leads to negative pressure. (Yang et al., 2008)...........................................................................45 18 Optical micrograph of a water plug in a 100 nm wide nanochannel (Tas et al., 2003......46 19 (a) Optical image (b) Schematic of wate r meniscus between a cantilever and a flat surface (Nosonovsky and Bhushan, 2008).........................................................................46 21 Arrangement of Pt wall atoms in fcc(111) structure; (a) 3D view of the first three layers, (b) XY projection (t op view) of the wall atoms....................................................57 22 Sketch for pressure contribution of atoms i j ....................................................................61 23 Simulation domain at an intermediate timestep for pressure verification.........................63 24 Temperaturedensity plot for li quid slab surrounded by vapor atoms...............................63 25 Pressuredensity plot for liq uid slab surrounded by vapor atoms......................................64 26 Surface tensionpressure plot for liquid slab surrounded by vapor atoms.........................64 27 Plot to obtain the threshold number dens ity and threshold radius for saturated liquid argon at 90 K.................................................................................................................. ....66 28 Fluidwall thermal equilibrium model...............................................................................68 29 Computational domain at an intermediate timestep for validation of proposed wall heat transfer model............................................................................................................ .71 9 PAGE 10 210 Variation of temperature and energy with time for validation of proposed wall heat transfer model (a) lower wall temperature is increased to 110 K while the upper wall temperature is decreased to 90 K; (b) lo wer wall temperature is increased to 130 K and the upper wall decreased to 90 K................................................................................72 211 1D heat problem...............................................................................................................74 212 Average temperature and energy vari ation with time of liquid argon time for validation of proposed wall heat transfer model................................................................76 213 Temperature gradient along zdir ection at different timesteps for Ti=100 K, Tup=110 K and Tlow=90 K (smooth lines represen t analytical solutions).........................................79 214 Temperature gradient along zdirectio n at different timesteps for h=23.657 nm, Ti=100 K, Tupp=120 K and Tlow=120 K (smooth lines represent analytical solutions)......80 215 Temperature gradient along zdirectio n at different timesteps for h=44.357 nm, Ti=100 K, Tup=120 K and Tlow=120 K (smooth lines represen t analytical solutions).......81 216 Average domain temperatures from MD simulation and analytical equation with heating from both platinum walls for different nanochannel heights................................82 217 Variation of number of liquid atoms with time in the simulation for verification with Maxwell Speed Dist ribution function................................................................................84 218 Velocity distribution of vapor argon atoms at fi nal equilibrium with a nonevaporating film at the lower wall.....................................................................................85 219 Variation of average temperature and en ergy with time for adiabatic wall with higher wettability J.............................................................................................86 210.89410220 Variation of average temperature and en ergy with time for adiabatic wall with lower wettability J..............................................................................................87 210.3581031 Snapshot of transient film evaporation domain.................................................................89 32 Variation of average system temperature and energy with time for transient film evaporation.........................................................................................................................90 33 Variation of liquid film argon atoms in the domain with time for transient film evaporation.........................................................................................................................91 34 Projection of the YZ plane of the simula tion domain at different time intervals for transient film evaporation..................................................................................................93 35 Number density profiles at different time steps for transient film evaporation.................94 36 Averaged temperature distribution pr ofile for transient film evaporation.........................94 10 PAGE 11 37 Averaged number density profile for transient film evaporation.......................................95 38 Evaporation flux and heat flux rate s for transient film evaporation..................................97 39 Sketch showing the coordinates for calc ulating the interactio n energy between (a) a Pt atom and Ar block (b) semii nfinite Pt block and Ar block..........................................98 310 Snapshot of simulation domain of nanochannel..............................................................101 311 Variation of average system temper ature and energy with time for case 4.....................102 312 Variation of vapor argon atoms with time for case 4.......................................................103 313 Snapshots of domain for case 4.......................................................................................104 314 Average density of domain at different times for case 4.................................................104 315 Densityheight plot showing nonevaporat ing liquid film thickness and equimolar surface position for case 6................................................................................................105 316 Heat flux rate with time for va rying nanochannel hei ght (cases 13)..............................107 317 Evaporation rate with time for vary ing nanochannel height (cases 13).........................107 318 Heat flux rate with time for varying film thickness (cases 2, 46)..................................110 319 Evaporation rate with time for va rying film thickness (cases 2, 46)..............................110 320 Pressuredensity plot of initial and final states for varying nanochannel height (cases 13)...................................................................................................................................112 321 Pressuredensity plot of initial and fina l states for varying film thickness (cases 2, 46)......................................................................................................................................113 322 Vapor pressure variation with time for case 7.................................................................114 323 Snapshots of domain at different time intervals for case 7..............................................115 41 Variation of drop temperature and ener gy with time for drop impact (case 1)...............119 42 Variation of dropcenter he ight and drop velocity with time for drop impact (case 1)...119 43 Snapshots of xz plane for drop impact (case 1)..............................................................120 44 Variation of contact angle with time for drop impact (case 1)........................................121 45 Snapshots of xz plane for drop impact (case 2)..............................................................122 46 Variation of drop center he ight and drop velocity with time for drop impact (case 2)...122 11 PAGE 12 47 Variation of contact angle with time for drop impact (case 2)........................................123 48 Snapshots of yz plane for drop impact (case 3)..............................................................124 49 Snapshots of yz plane for drop impact (case 4)..............................................................124 51 Schematic showing configuration of the computational domain.....................................127 52 3D schematic of the computational domain.....................................................................128 53 Schematic depicting division of nanochanne ls and liquid meniscus into 11 regions......130 54 Schematic showing two regions of an evaporating meniscus at times t1 and t2..............131 55 Snapshots of xz plane at different time steps for Case I.................................................134 56 Variation of liquid atoms with time for regions 510 for Case I.....................................135 57 Variation of vapor pressu re with time for Case I.............................................................137 58 Average temperature of regions at initi al and final time periods for Case I....................137 59 Average heat and evaporation flux rate s for regions 110 during the total heating period for Case I...............................................................................................................139 510 Heat flux distribution al ong nanoand microlayers at the base of the bubble...............139 511 Variation of disjoining pressure along meniscus length for Case I.................................141 512 Variation of capillary pressure along meniscus length for Case I...................................141 513 Snapshots of xz plane at different time steps for Case II...............................................143 514 Variation of liquid atoms with time for regions 510 for Case II....................................144 515 Average temperature of regions at initi al and final time periods for Case II..................146 516 Average heat and evaporation flux rate s for regions 110 during the total heating period for Case II............................................................................................................. 146 517 Variation of disjoining pressure along meniscus length for Case II................................147 518 Variation of capillary pressure along meniscus length for Case II..................................147 519 (a) Schematic of a meniscus to com pute forces acting by the extended meniscus region on liquid atoms in the bounded menisc us region, and (b) configuration for force interaction between an atom and a rectangular region...........................................150 12 PAGE 13 520 Schematic showing meniscus and nanoc hannels used in the present work to determine the force distribution along th e meniscus length by the nanochannels...........151 521 Force distribution on the meniscus by the nanochannel as a function of the distance from the nanochannel (a) Fx and Fy, (b) Fz......................................................................152 61 Schematic showing multiscale method of solution for proposed work..........................158 62 Domain Decomposition method......................................................................................162 13 PAGE 14 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 NANOSCALE TRANSPORT MECHANISMS ASSOCIATED WITH THREEPHASE CONTACT LINES By Shalabh Chandra Maroo December 2009 Chair: Jacob N. Chung Major: Mechanical Engineering Nucleate boiling is one of the most efficien t heat transfer mechanisms and has been widely featured in many high heatflux engine ering applications. But, the knowledge on the underlying basic transport mechanisms is quite limite d and rather incomplete, due to the lack of fundamental understanding of the physics occurring at the threepha se contact line. The contact line region comprises of noneva porating thin film region, evapor ating thin film region and intrinsic meniscus region. The nonevaporating thin film is the main focus in this work. Its thickness is expected to be in th e size of nanometers, and thus is studied here using molecular dynamics simulations. A novel fluidwall thermal equilibrium model, to numerically simulate heat transfer between wall and fluid atoms, is proposed and verified against available experimental data and analytical methods. Nanoscale phase change process is studied by simulating evaporation of a thin liquid argon film on a platinum wall using the proposed model. A nonevaporating film is obtained on the pla tinum surface, and its Hamaker constant is evaluated. Additional simulations of thin film evaporation are done in a nanochannel to study the effect of varying nanochannel height and film thickness. Hamaker constant, vapor pressure, film thickness, and net evaporation and heat fluxes are evaluated. It is conclude d that the creation of nonevaporating film is dependent on vapor pressure, substrate temperature and solidliquid 14 PAGE 15 molecular interaction strength. The vapor pres sure and substrate temperature are dominant factors. Next, a nanoscale evaporating meni scus is simulated. Heat and mass transfer characteristics and pressure variation in the nonevaporating and inter line regions are studied. High disjoining pressure a nd negative capillary pressures are ob tained. It is shown that cavitation cannot occur as the capillary size is smaller than the critical cavitation radius, and the meniscus can exist in metastable state. A curvefitted menisc us boundary condition is developed; a force function of the form can be applied at the boundaries of a liquid film to create curvature and form a meniscus. Also, impact of argon nanodroplet on platinum wall with varying wettability is simulated to study the variation of contact angle with time. 3 nFAnCn2 15 PAGE 16 CHAPTER 1 INTRODUCTION AND LITERATURE REVIEW Introduction The advent of the computer sparked the intere st in simulations on a molecular scale. Not many years after the digital computers appeared th at the very first molecular simulations were performed. Since then, the spectacu lar growth in the computationa l power of the computer has made possible the usage of molecular simulations in solving problems of increasing complexity. As every substance is made from elementary particles (atoms/molecules) and if the basic dynamic parameters of these atoms (i.e. positi on, velocity and interaction force) can be determined, the physical properties of the substa nce, like temperature, pressure, etc can be obtained from those dynamic parameters using st atistical methods. This forms the basis of Molecular Dynamics simulations. Using molecula r dynamics, it is possible to solve a dynamic problem without making any approximations, with in the limits exposed by machine accuracy. All molecular dynamics simulations follow Newton s second law in its simple form for simple atomic systems or in a more generalized form (Lagrangian/Eulerian formulation) for complex geometry systems. Thus, molecular dynamics can also be understood as a virtual experiment since the movement of atoms occu rs naturally and the system e volves with time as in a real experiment. Molecular dynamics simula tion consists of three basic parts: 1) A model for interactions between various system components like atoms, molecules, surfaces, etc which is commonly known as molecular interaction potential function. 2) An integrator method to obtain the positions a nd velocities of the mol ecules as the system evolves with time. 3) A statistical analysis of the dyna mics parameters to obtain prop erties such as temperature, energy, etc 16 PAGE 17 Figure 11 shows a basic molecular dynamics algorithm flowchart. The force on each molecule by all other molecules in the system can be determined by using the molecular interaction potential function. From this force, th e equation of motion can be integrated over time to obtain the new positions and velocities of each molecule. This information can be statistically analyzed by space/time averaging to provide macr oscopic physical properties, and thus the evolution of the system can be studied. It should be noted here that molecular dynami cs approximates electron ic distributions in a rather coarsegrained fashion by putting either fixed partial charges on inte raction sites or by adding an approximate model for polarization effect s. In both cases, the time scale of the system is not dominated by the motion of electrons, but by the time of intermolecular collision events, rotational motions or intramolecula r vibrations which are orders of magnitude slower than those of electron motions (Sutmann 2002). Most of the mo lecular dynamics simulations till date have been restricted to nano length and time scales as they strongly depend on computing power available. Molecular dynamics is nowadays applied to a huge class of problems, e.g. phase change, properties of liquids, defects in solids, fracture, surface properties, friction, molecular clusters, polyelectrolytes and biomolecules. F U Figure 11. Flowchart showing the basic st ructure of molecular dynamics simulation 17 PAGE 18 Need for Molecular Dynamics Simulation Nucleate boiling is one of the most efficien t heat transfer mechanisms and has been widely featured in many high heatflux engineering applications. In many engineering applications, nucleate boiling heat transfer is the mode of choice. Boiling heat transfer has the potential advantage of being able to transfer a large amount of energy over a relatively narrow temperature range with a sma ll weight to power ratio. The critical heat flux, the maximum nucleate boi ling heat flux, is closely related to the design and safe operation of boiling system s. But the knowledge on the underlying basic transport mechanisms is quite limited and rather incomplete, that has hindered development of a universal theory and corresponding model govern ing the boiling process from bubble dynamics and heat transfer during bubble gr owth and departure to critical heat flux phenomenon. This lack of fundamental understanding is the sole reason behind this work that is intended to unlock the core of the transport physics in boiling a moving contact line with simultaneous heat and mass transfer. The moving contact line is a threephase phenomenon wher e immiscible gasliquid interface meets the solid surface. Moving contact lin e is involved, usually playing the dominant role, in many natural and engineering processes such as rivul et, coating flows, and droplet impingement. Spray droplet cooling and hete rogeneous nucleate boiling, both are highly effective heat transfer mechanisms, are dominated by the moving contact line with simultaneous heat and mass transfer around the contact line region. During nucleate boiling from a heated surface, vapor bubbles are nucleated and they st art to grow owing to infusion of vapor mass generated mainly from the moving contact line region. The moving contact line moves outward during bubble growth period while it retracts back during the bubbl e departure stage. At higher heat fluxes, bubblebubble interact ion and coalescence are responsible mechanisms for achieving 18 PAGE 19 high heat transfer rates in heterogeneous nucle ate boiling where moving c ontact line also plays an important role. Bubblebubble co alescence creates strong disturba nces to the fluid mechanics and heat transfer of the microlayer and the moving contact line beneath the bubbles. Because of the complicated nature, the detailed physics a nd its effects of moving contact line on the bubblebubble interaction process have not yet been unveiled. The critical heat flux on eart h, which is the upper heat transfer limit in the nucleate boiling regime, represents a st ate of momentum and mass balance. Because the buoyancy force strength is relatively constant on earth, for heat fluxes lower than th e critical heat flux, this force is more than that required for a complete rem oval of vapor bubbles formed on the heater surface. At the critical heat flux, the buoyancy force is exactly equal to the force required for a total removal of the vapor bubbles. For heat fluxes gr eater than the critical heat flux, the buoyancy force is unable to remove all the bubbles, t hus resulting in the accumu lation and merging of bubbles on the heater surface, which simultaneously leads to a total blan keting of the heater surface by a layer of superheated vapor. It is believed that bub ble departure and vapor removal from the heater surface is dominated by vapor generation and bubble coalescence, both of which are hinged on the moving contact line dynamics, as a result, the moving c ontact line is closely related to the critical heat flux. The noslip boundary condition between a flui d flow and a solid surface is generally regarded as a cornerstone in the continuum hydrodynamics, owing to its proven applicability in diverse fluid flows. However decades ago, it was discovered that in immiscible gasliquid twophase flows over a solid surface, the moving contact line, defined as the intersection between the moving gasliquid interface and th e solid surface, is incompatib le with the noslip boundary condition (Moffatt 1964, Huh and Scriven 1971, Dussan and Davis 1974, Dussan 1976 and 19 PAGE 20 1979). These reports showed that there is veloci ty discontinuity at the MCL which is now wellrecognized as the contactline singularity (Shi khmurzaev 2006 and Blake 2006). In the past two decades, it has been confirmed through molecular dynamic simulations that a nearcomplete slip indeed takes place at the MCL (Koplik et al. 1988, Thompson and Robbins 1989, and Thompson et al. 1993). This finding pres ented a challenge for the fluid m echanics community, due to a lack of viable solutions apart from ad hoc fixes that include analytical models (Phan et al. 2006, Qian et al. 2006, and Fuentes and Cerro 2007) and nume rical models (Ding and Spelt 2007a,b). In the absence of a viable boundary condition which can reproduce the molecular dynamics results, an accurate continuum description of the moving cont act line flows at the microor nanoscale remains an elusive goal which may never be solv ed explicitly by contin uum fluid mechanics. Literature Review Nucleate Boiling It has been widely accepted in literature that the contact lin e at the base of the bubble in nucleate boiling can be divided into macroand microregions. Th e micro region is the ultrathin liquid film between the solid su rface and the evolving liquidvapor interface. The macro region is the region occupied by vapor and liquid, except the microlayer. The microregion can be subdivided into two additional re gions. Thus the three regions, as shown in figure 12, are: a) Nonevaporating thinfilm region liquid is adsorbed on the heater surface and forms a nonevaporating layer with molecular forces having controlling influence b) Evaporating thinfilm region maximum evapor ation and heat transfer occurs in this region and the liquid is fed from the bulk liq uid through the intrinsic meniscus region c) Intrinsic meniscus region fluid mechan ics in this region is governed by the conventional equation of capillarity 20 PAGE 21 There is a region between the nonevaporating thin film region a nd the evaporation thin film region over which the film varies in thickness and curvature to accommodate the transition between the two regions. This is called the inte rline region and is the thinnest portion of the meniscus over which vaporization can occur. Since it is the thinnest, it is also the location where the evaporation rate is the highest. Figure 12. Schematic of bubble growth, (a) overall picture at macroscale (shaded region depicts the area covered by nonevapora ting thin film), and (b) zoomed in nanoand microscale regions at the th reephase contact line Boiling has been studied extensively during the last few decades which include experimental work, and theoretical models to predict the exact nature of boiling in various conditions. Since the boiling process is very complex, its complete understanding still poses significant challenges and researchers have not converged on a defin ite precise model. The past efforts can be categorized into three broad groups (Dhir 2006): 1) Empirical correlations experi mental data was obtained for boiling heat transfer as a function of several independent variable like heater size, heating method, heater thickness, system pressure, liquid temperature, flow conditions, etc. This data was correlated by many researchers. Although the correlations have been helpful in application to practical situations, they serv e their useful purpose only when applied in the range of the database used in developing the correlations. 21 PAGE 22 2) Mechanismbased correlations This involved studying some of the individual subprocesses to develop mechanistic models for predicting heat fluxes as a function of wall superheat or other independe nt variables. An expression for partial nucleate boiling heat flux was written based on the contribu tions of three mechanisms by: transient conduction into the liquid during the waiting period over the region influenced by vapor bubbles, natural convection over th e remaining regions of the h eater and evaporative heat transfer from the microlayer underneath the bubble. The expr ession required the knowledge of three key parameters: number de nsity of active sites, bubble diameter at departure and bubble release frequency. A numb er of attempts have been made to develop correlations or models for th ese parameters, but with limited success. 3) Numerical simulation In order to have a credible predictive m odel of nucleate boiling, one must include all subprocesses which occur and not guess how one parameter influences the overall process. Thus, a comple te numerical simulation of the process is asked for a tool which has been develope d only recently. The region of interest is divided into microand macroregions as me ntioned before. Lubrication theory has been widely used to solve for the radial variati on of microlayer thickness, while complete conservation of mass, momentum and en ergy are solved for both phases in the macroregion. As reported in Chen and Chung (2002, 2003a), a t ypical ebullition cycle for a single bubble is composed, as shown in Figure 13, of two basic stages: (1) bubble nuc leation (Step A) which takes place spontaneously with the departure of the proceeding bubble (Step F), and (2) bubble growth (Steps B. C, D and E). There is only one heat flux spike for a single bubble ebullition cycle that is the result of the combined bubble departure and nucle ation process. The boiling heat 22 PAGE 23 transfer is closely associated with the bubble's contact line moveme nt and dry area on the heater surface. Chen and Chung (2003b) and Chen et al. (2004) further examined the effects of heater size on the boiling characteristics. They conclu ded that the boiling cu rve obtained from the microheater is composed of two regimes which ar e separated by a peak heat flux. It is suggested that in the lower superheat regime, the boiling is dominated by liquid re wetting and microlayer evaporation, while in the highe r superheat regime, conduction th rough the vapor film and microconvection plays the key heat transfer role as the heater is covere d by vapor all the time. Time(seconds) Heatflux(w/cm2) 0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 20 30 40 50 60 70 80 90 100[A] [B] [C] [D] [E] [F] [F],[A]:departureandnucleation [B],[C],[D],[E]:growthHeater#1at110 C bullition Cycle Liao et al. (2004) presented a physical model to predict the early st age bubble growth in cleate boiling. The microlay er was taken to be in the form of a wedge. mal boundary layer is revealed by the asymptotic and num A B C D E F Figure 13. Single Bubble E saturated heterogeneous nu The structure of the thin unsteady liquid ther erical solutions. The author s find that in some cases, the heat transfer from the thin 23 PAGE 24 unsteady thermal boundary layer near the rapidly growing bubble can be larger than that from the microlayer. Son (2001) proposed a numerical method based on a level se t technique for capturing th phase interface w e hich included the effect of phase change at the interface as well as to achieve mass co d. The microregion was wedge a t. The coupled algorith e equatio o ve to nservation. The modified level set met hod was applied for numerical simulation of bubble rise and growth in a stati onary liquid. Wu et al. (2007) c oupled the level set method with the moving mesh method to simulate subcooled nucleate boiling. Genske and Stephan (2006) modeled the macroregion by solving the NavierStokes equation for both vapor and liquid phase using finite element metho shaped and moving mesh was used for free surface simulation. Flow inside and around growing bubble and its influence on the ove rall heat transferr was considered. Tomar et al. (2005) used the coupled levelset and volumeoffluid method to simulate bubble growth in water at near cr itical pressure for different de grees of superhea m was said to conserve ma ss and capture the complicated in terfaces very accurately. Th effect of superheat and saturation pressure on frequency of bubble formation was analyzed. Mukherjee and Dhir (2004) numerically studied bubble dynamics and associated wall heat transfer during lateral bubble merger. Th ey solved the incomp ressible NavierStokes n and captured the liquidva por interface using the level se t formulation. Mergers of tw bubbles, three bubbles in a line and three bubb les in a plane were simulated assuming axisymmetry. The contact angle between the wa ll and the bubble was taken to be constant throughout the simulation. The authors acknowledge th e fact that the resu lts are sensiti contact angle. 24 PAGE 25 Mukherjee and Kandlikar (2007) used a sta tic contact angle and dynamic contact ang models based o le n the contact line velocity and th e sign of contact line ve locity to numerically study th e be ular l cannot be achieve nd on a ics in a na nochannel. Liquid argon was placed between parallel n ined in a nanochannel with emphasis on nucleation phenomenon. A etasta ty. The pic e growth and departure of single vapor bubbles from a heated wall. The liquidvapor interface is simulated using the levelset technique. The effect of microlayer evaporation at th bubble base was neglected. The authors concluded th at correct motion of the contact line can simulated when the dynamic contact angle is re presented as a function of interface velocity. However, the numerical model was not va lidated with experimental results. Recently, a few researchers attempted to obs erve nanobubble nucleation using molec dynamics simulations. Vapor nucleation in a li quid confined in a nanochanne d by increasing the temperature of the wall on ly. Hence, the studies were associated with varying system pressure by using metastable liquid, increasing the nanochannel volume a thermally contract the liquid by lowering the te mperature. The nucleation occurs when the pressure reached the spinodal line. Maruyama and Kimura (2000) simulated hete rogeneous nucleation of vapor bubble solid surface using molecular dynam solid platinum surfaces and graduall y expanded until a vapor bubble was nucleated. The wettability of the lower surface was varied by changing the potential parameter between argo and solid molecule. Nagayama et al. (2006) carried out molecular dynamics simulations to examine the bubble behavior conf mble liquid confined in a nanochannel w ith an inlet driving force was simulated. Nano sized bubbles were observed wider different conditi ons of solidliquid inte rfacial wettabili authors also found that the vapor pressure inside the nanobubble disagrees with the macrosco 25 PAGE 26 YoungLaplace equation. Park et al. (2001) formed a similar bubble using molecular dynamics with the liquid initially in metastable condition. No solid surfaces were included in this simulation. The surface tension of the bubble, calcul ated from the density profile and from the normal and transverse components of pressure te nsor, was found to increase very slightl decrease of bubble radius. Schoen et al. (2005) investigated the phase transition of a simple liquid bounded betwee two parallel walls a few nan y with n ometers apart with molecular dynamics simulation. The authors used l eating with a constant mean negative pressure using molecular dynamics. Velocityscaling nd o species of flui d particles between two solid walls and then shearing the walls with a constant speed in opposite directions. This simple geometry was chosen the fact that since some fluids (argon in their case) thermally contract on cooling, the lower wall temperature was reduced to the extent that the fl uid started to crystalli ze starting from the bottom surface creating a frozen bubble in crystallized state. It was found that the confining walls primarily affect the crystallization process. Th e authors also noted that vapor nucleation in a liquid confined in a nanochannel cannot be achieved by increasing the temperature of the wal only. Wu and Pan (2003) studied bubble nucleati on in a homogeneous LennardJones liquid under h temperature control is used to raise th e system temperature while the system volume is expaned to keep the mean pressure nearly cons tant. No solid surfaces were included in the simulation. Bubble nucleation is observed as th e system pressuretemperature reaches the spinodal line. The authors find very good agreemen t of nucleation rate between simulation a prediction from classical theory. Ren and E (2007) studied the physical pr ocess near a moving contact line using molecular dynamics by placing tw 26 PAGE 27 based o ion at d el tact of the major focus of boiling heat transfer research. cal correlations for the critical heat flux are now available in the literature, althoug of a liquid layer formed on a heat ing surface. They also suggested that a liquid macrolayer is n the assumption that the outer flow does not change the physics of the inner reg this length scale. They evaluated the friction force and stresses around the contact line region an suggested boundary conditions, based on force balan ce, at the contact line and at the wall away from the contact line. The authors implemented these boundary c onditions in a continuum mod by solving the incompressible NavierStokes equation for a similar setup by having two immiscible fluids confined between two parallel solid walls with external force acting on the fluid. The two fluids interact identically with the wall, and values for the different parameters in the boundary conditions were assumed. The depe ndence of the apparent and dynamic con angles on the contact line velo city and friction coefficient wa s studied. The authors concluded that: (a) shear rates are well de scribed by linear models for all practical purposes, (b) Youngs stress dominates over the viscous stress inside the contact line regi on, and (c) the Navier boundary condition describes the friction dynamics at the fluidsolid in terface away from the contact line region accurately. The authors point out the fact that further work needs to be done in order to study more complex processes. Critical Heat Flux The boiling curve was first identified by N ukiyama (1934) more than sixty years ago. Since then the critical heat flux has been one A plethora of empiri h each is applicable to somewhat narrow ranges of experimental conditions and fluids. There have been three mechanistic models suggested in the past : 1) vapo r escape path instability model (Zuber, 1959, Lienhard and Dhir, 1973), 2) macrolayer dryout model (Haramura and Katto, 1983) and 3) vapor stem merg ing model (Dhir and Liaw, 1989). Sakashita and Kumada (1993) proposed that the critical heat flux is caused by the dryout 27 PAGE 28 formed due to the coalescence of bubbles for most boiling systems, and that the dryout of the macrolayer is controlled by the hydrodynamic beha vior of coalesced bubb les on the macrolayer. Based o o he s satisfactory overall mechanistic descrip scientific significance. The act of g can be seen in a wide range of applications like spray cooling ic n these considerations, a new critical heat flux model is proposed for saturated pool boiling at higher pressures. In the model, they su ggest that a liquid macrolayer is formed due t coalescence of the secondary bubbles formed from the primary bubbles. The detachment of the tertiary bubbles formed from the secondary bubbles determines the frequency of the liquid macrolayer formation. The critical heat flux occurs when the macrolayer is dried out before t departure of the tertiary bubbl es from the heating surface. Sadasivan et al. (1995) suggest ed that identifying the physical characteristics such as active nucleation site distributi on, static versus dynamic contact angles, and advancing versu receding angles would help understand the heater surface rewetting behavior, thus again showing the importance of microlayer region. The consensus is that a tion for the critical heat flux still remains elusive. Droplet Impact and Spreading Droplets are encountered in daytoday life and the phenomena related to them are so common that they are usually overlooked. The ba sic physics of droplet dynamics is extremely intriguing and understanding them in detail is of immense droplet impingement and spreadin of surfaces (semiconductor chips, elect ronic devices, etc), spray combustion, inkjet printing, lubrication, crop spraying, etc. Thus, droplet impact dynamics can be accompanied by mass, momentum and energy transfer. The sp reading phenomenon is governed by the dynam contact angle and the moving contact line. Full characterization of the effect of the surface wettability and roughness on droplet growth is still to be achie ved, and there is no complete 28 PAGE 29 agreement on the boundary conditions to be impl emented at the moving contact line. Oblique impacts of droplets are also insufficiently studied and understood. Droplet impact and spreading on horizontal surfaces have been studied on the continu level by many researchers. Fukai et al. (1993) studied the deform ation of a spherical droplet impingement on a surface using deforming finite elements and grid um generation. The singularity at the con t thor around d sing apes compared well with the previous experimental results. tact line is removed by imposing a net in terfacial force given by the equilibrium surfacetension coefficients of the join ing phases during flow with the capillary force neglected at the contact line. They used two fl uids: water and tin in their simulations and showcased the occurrence of droplet rebounding and mass accumu lation around its periphery on spreading. Davidson (2000) analyzed the axisymmetric sp reading of an inviscid drop impingemen on a solid surface using a boundary integral method for different Weber numbers. The au concluded that surface tension forces are responsible for the accumulation of liquid in a rim the periphery of the spreading drop for low Weber numbers, whereas viscous forces govern the accumulation for high Weber numbe rs. Reznik and Yarin (2002) studied the spreading of an axisymmetrical drop on a horizontal surface under the influence of gravity an surface tension for different Bond numbers by assuming inertial forces to be negligible and u the creeping flow regime approximation. Francois and Shyy (2003) applied the imme rsed boundary method to simulate droplet impact on a flat surface with heat transfer. Bo th the static and dynamic contact angle models were considered. The computed 3D drop sh The wall heat flux distribution under the droplet was governed by its dynamics with a larger spread resulting in a wider area of heat exchange. 29 PAGE 30 Fedorchenko et al. (2005) studied the lates tage of drop spreading by modifying the she flow equations to account for th e viscous drag and obtaine et d an analytical solution. They also propose spherical cap wh ) (2007) simulate d the droplet impact and spread ing using the volumeoffluid erical (1999) ct d a simplified kinematic model which take s into account viscous and capillary effects and shows the influence of Reynolds and Weber numbers on drop spreading dynamics. Park et al. (2003) developed a theoretical model based on energy ba lance to predict the maximum spreading ratio at low impact velocity assuming the shape of the drop to be a ile spreading on the surface, and validated their model with experimental results showing good agreement between the two. Fukai et al. (1 998) studied the maximum spreading diameter and spreading time for droplet impact on a surface by solving the NavierStokes equations for an axisymmetric system and neglecting the gravitati onal effect for a wide range of Reynolds and Weber number. PasandidehFard et al. (1996), Bussmann et al. ( 1999), Sikalo et al. (2005) Gunjal et al. (2005 and Lunkad et al (VOF) method. PasandidehFard et al. (1996) solved the Na vierStokes equation using a modified SOLAvolume of fluid method to m odel drop impact on surfaces. Experimentally measured values of dynamic contact angles were used as a boundary condition for the num model. The computer generated images matc hed well with photographs. Assumption of a constant contact angle made the model prediction s less accurate. A simple analytical expression was developed to estimate the maximum spread of a droplet on a surface. Bussmann et al. developed a 3D model of droplet impact onto asymmetric surfaces. Two cases are studied: oblique impact of the droplet onto an inclined surface and onto a sharp edge. Contact angles, measured from photographs, are applied as boundar y conditions at the contact line. A simpler model is proposed which requires only the values of the rapidly advanci ng and receding conta 30 PAGE 31 angles. The simulation results matched well with the photographic data. Sikalo et al. (2005) assumed that the liquid and gas pha ses share the same velocity fi eld. They noted that numerical simulations which assume a fixed advancing/re ceding contact angle do not predict the receding phase of the droplet spreading well. The singularity near the cont act line was replaced by a local force with some defined dependence on the in stantaneous advancing/receding contact line velocity. Gunjal et al. (2005) not ed that the agreement between the experimental and numerical results improved when the average contact angle variation with time was implemented inste using an equilibrium contact angle value. Lunka d et al. (2007) used both dynamic contact angle (DCA) and static contact angle (SCA) models. For horizontal surfaces, they concluded that results from both the models matched well with experiments for less wettable surface whereas for highly wettable surface the DCA model is neede d. They also noted that SCA models fail predict the observed partial or full rebound regi mes of droplet impact on inclined surfaces. As can be seen from the different con tinuum simulations of droplet spreading, no complete consensus has been reached so far on the dynamic of the problem. This can be ow ad of to ed due to t ve not es in mbersome he fact that it is governed by a num ber of factors and comprises of many different phenomenon occurring at diverse length and ti me scales. Among the various parameters, the variation of the contact angle a nd the motion of the contact line ar e most intriguing and ha been clearly understood. The divergence in the solution caused by the noslip model near the contact line has lead to researchers proposing different boundary conditions to remove the singularity but their physical origin and valid ity is not clear. Some researchers have experimentally measured the advancing and recedi ng contact angle, and used these two valu their models to account for the variation in the contact angle. This pro cess is both cu and still an approximation as the dynamic cont act angle changes continuously with time. The 31 PAGE 32 microscopic region near the contact line is experimentally inacce ssible. It should also be noted that at the actual molecular level, the three phas es meet in a small region with finite dimension and not at a line as assumed in continuum analysis. Thus this mi croscopic contact angle may be different from the experimentally measured cont act angle (Sikalo et al., 2005). All these facts lead to the necessity of molecular modeling in th e region near the contact line. With the advent of highspeed computing, molecular dynamics has become an increasingly important tool in gaining a better understanding of fundamental processes by simulation each atom and performing virtual experiments. It becomes po ssible to vary any governing parameter and to study its dominance on the overall process. Franco is and Shyy (2003) acknowledged th multiscale modeling that couples the continuum simulation with molecular dynamics in regio surrounding the contact line as a more fundamental approach to understand droplet dynamics. Maruyama et al. (1998) used molecular dynami cs to simulate equilibration of an argon liquid droplet on a platinum wall. The wall was represented by three ty pes of models: three s e need of ns layers o l. The eled with the SPC/E potential and the platinum wall was represented by three layers of harmonic molecules modeled as a f harmonic molecules, one layer of fixed molecules and a simple onedimensional function. The fluidfluid and fluidsolid interaction was governed by the LennardJones potential. The contact angle was measured by fitt ing a circle to the twodimensional density profile. The contact angle was co rrelated with the integrated de pth of the surface potentia effect of temperature on the contact angle coul d not be obtained. The authors also studied the simultaneous evaporation and condensation of dr oplets, and concluded that their measured contact angles were the same as that at the equilibrium condition. Kimura and Maruyama (2003) studied the sp reading of a water droplet on a platinum surface using molecular dynamics. The water molecu les were mod 32 PAGE 33 constan t early g rating the potential energy between an LJ particle ium which interacts strongl rs, h f t temperature heat bat h. The waterplatinum interaction was simulated by two differen potential functions. The spreadi ng rate of the water droplet was determined by measuring the area of the waterplatinum interf ace, and was proportional to the onethird power of time at stage and onefifth power of time at the later stage. Ingebrigtsen and Toxvaerd (2007) determined th e contact angles of droplets constitutin of LennardJones particles on a solid surface usin g molecular dynamics simulations. The surface, called the 3 LJ potential, was modeled by in teg s and a semiinfinite co ntinuum of uniformly distribute d LJ particles. The domain volume was divided into parallel sheets and the local density variation was obtained using an empirical tanh function in order to determine the equilibrium contact angle. The angles determined matched fairly well with those obt ained from Youngs equa tion for medium and weak fluidwall interactions but over predicted them for strong interactions. Adao et al. (1999) used molecular dynamics simulations to show that the equilibr contact angle of a sessile drop placed on a heterogeneous substr ate follows Cassies law. The homogeneity is modeled by having the substrate composed of two species, one y with the liquid and other which interact s weakly. The LennardJones potential is used for all interactions. The contact angle is determin ed by dividing the droplet into horizontal laye computing the density of the particles as a f unction of the distance from the centre and then locating the extremity of the laye r as the distance where the density falls below a cutoff value of 0.5 times the liquid density. The authors concluded that the equilibrium contact angle does not change if the sizes of the heterogeneities are of molecular scale, and that the dynamics of suc MD drops can be described by the molecular kine tic theory. Gentner et al. (2004) extended this work to include drop impact for cases having low inertia with the droplet molecule consisting o 33 PAGE 34 16atom chains. The shape of the drop was fitted by a sphere for zero impact velocity, and approximated by a spheroid when the initial spee d is not zero. They measured the associated base radii and contact angle versus time and concluded that the spreadi ng velocity is a unique function of the driving capillary force. Thus, molecular dynamics has been used in the recent years to model drop impact and spreading providing some insight into the moving contact line and the over all process. Droplet spreading coupled with evapor ation, ob lique drop impacts, droplet condensation are a few problem ng the fluidsolid wall boundary co ndition across which the heat transfer takes place. N rd of first two planes of an fcc s s which need to be captured accurately using molecular dynamics. Extensive effort is still required to understand the physics and fundamentals of the phenomenon. Wall Models As molecular dynamics is becoming an increas ingly important tool in gaining a better understanding about numerous heat transfer research problems many problems in this field require modeli umerous models have been pr oposed and developed in literature. Abraham (1978), in a pioneering study, used four wall models to study the interfacial fluid density profile of solid/liquid interface sy stems using Monte Carlo simulations: a) Lenna Jones (LJ) atoms were placed and constrained to rema in at the lattice sites olid termed LJ (100) wa ll b) discrete centers of intera ction in each solid plane were replaced by a continuous distribut ion of LennarJones atoms with uni form planar density leading to the LJ 104 wall model, c) proposed B oltzmann weighted wall whose potential is proportional to the probability that a single mole cule is at a distance z from the (100) wall irrespective of its x, y co ordinates, d) the Hard wall whose potential is infinite if an atom penetrates the wall structur e, and zero otherwise. 34 PAGE 35 Xu and Li (2007), Thompson and Troian (199 7), Priezjev (2007), Voronov et al. (20 arranged the solid wall atoms in a crystalline struc 06) ture, and remained fixed during the whole comput ws heavy c st ructure and remain fixed during the simulation. Accord et n by random g them ely linder. The ation process. The fluid atoms were coupl ed with a thermostat to assign a specific temperature to the fluid. Koplik et al. (1989), to study microscopic aspects of slow viscous flo past a solid wall, set the wall atoms into an initi al fcc lattice structure and assigned them a mass mH=1010m allowing them to move according to the equations of motion. The walls would retain their integrity during the c ourse of the simulation as well as conserve energy in collisions, but eventually disintegrate theoretically. Xu and Zhou (2004) studied Poiseuille flow in a nanoscale channel using a thermal wall model. The wall atoms are arranged in a fc ing to the model, when a liquid particle st rikes the wall, all three components of velocities are reset to a biased Maxwellian distribution and th e direction of the zcom ponent velocity is s away from the wall; the criteria for liquid at om striking the wall was not mentioned. Markvoort and Hilbers (2005) looked into the dependence of wetting on solidgas interface for different density gases. The walls were constructed in a different simulatio ly placing the wall atoms, governed by LJ potential, in a simula tion box and givin a high temperature. The domain was then cooled down and the system crystallized which was placed in a wider box forming a wall. None of the wall atoms were fixed or restricted. However, it was seen that as the mass of the wall was so la rge compared to a single gas particles that a single collision hardly affected the wall and the wall atoms keeps their positions. Drazer et al. (2005) used molecular dynamics simulation to study the behavior of clos fitting spherical and ellipsoidal particles moving through a fluidfilled nanoscale cy 35 PAGE 36 cylindr all and allows d sites by a harm onic restoring force. The wall atoms interacted by the LJ potenti a ing t. In nd solid wall as (111) fcc structure made up of three layers o ical wall was compos ed of atoms of mass mw=100 tethered to fixed lattice sites by a stiff linear spring. Heat transfer between wall and fluid was not considered in the model. Spijker et al. (2008) develope d a new particle wall boundary condition to replace solid walls. The new condition is based on averaging the contributions of an explicit solid w for different crystal lattices to be in cluded. Wemhoff and Carey (2005) replaced the soli wall by deriving a wall potential to bind a liquid argon film to the imaginary wall, in order to reduce computational cost. The wall potential depends on the approximated values of Hamaker constant for metalmetal and fluidfluid interac tion. The system pressure was adjusted to an equilibrium value using a flux boundary. The model did not include heati ng/cooling of fluid by the wall atoms. Yang (2006) simulated a plane wall where th e wall atoms were anchored to their respective lattice al, and a Gaussian isokinetic thermostat was used to adjust the wall temperature to constant value at each time step. Ziarani and Mohamad (2008) used a similar method of spr force but used velocity rescaling to achieve cons tant wall temperature instead of a thermosta addition, the first five molecules at the two end of each layer were completely frozen to avoid walls being washed away by fluid molecules. Maruyama and Kimura (1999), Ohara and Suz uki (2000), Yi et al. (2002), Nagayama a Ceng (2004), Schoen et al. (2005) modeled the f solid atoms governed by a harmonic pot ential. Outside the three layers, phantom molecules were used to maintain the solid atom s at a specified mean temperature, resembling a Langevin thermostat. 36 PAGE 37 The reason that such various models have b een introduced in the past is that no one particular model has been accepted by all resear chers; the biggest concern being the actual physics n the liquid molecules and the molecules to the layer as if the pressure in the l the rating anoof the any heat transfer model may not be properly depicted. Wall models need to be thoroughly validated to show their physical soundness. Ultrathin Films, Disjoining Pressure and Hamaker Constant For ultrathin films, attractive forces exist betwee of the solid surface. These attractive forces act to pull the liquid in ayer were reduced by an amount Pd, known as the disjoining pressure. Variation of disjoining pressure with film thickness is s hown in figure 14. A th in film can exist in equilibrium on a solid surface even when its vapor pressure in the surrounding gas is below normal saturation pressure for the system temperature (Carey, 1992). Such a nonevapo film can be seen in the microlay er region at the contact line in nuc leate boiling at th e base of the bubble as discusses before and shown in figure. The thickness of a nonevaporating film is a function of the surrounding vapor pressure a nd substrate temperature (Wayner, 1999). The thickness of these nonevaporating th in films is expected to be in the size of nanometers, and thus can be studied using mol ecular dynamics. Only a handful of researchers have studied n thin film evaporation on a solid surf ace using molecular dynamics simulation. Figure 14. Expected variation of film thickness with film thickness 37 PAGE 38 Maruyama and Kimura (1999) applied molecu lar dynamics simulation to study thin film evaporation and condensation in a nanochannel. An argon vapor region was sandwiched between two argon liquid layers which were in contact with two platinum solid walls. The platinum wall was represented by three layers of harmonic molecules, and the temperature of the wall was controlled by a layer of phantom molecules wh ich are excited by a random force of Gaussian distribution. The aim of the study was to demonstr ate that thermal resist ance cannot be neglected of t temperature using the Langev otential over a solidliquid interface when th e system size is very small. Wu and Pan (2006) simulated a thin argon film evaporation on Pt surface, without the temperature effect, in vacuum. They obtained th e evaporation/condensation coefficient; noticed the mean temperature for the interface region to be 23 K lower than that inside the liquid and the interface thickness to increase with temperature. They concluded that the net evaporation rate thin liquid film into vacuum in a closed sy stem may be modeled based on Schrage approach. Yi et al. (2002) studied the evaporation of a thin liquid argon layer on a Pt surface in an argon vapor atmosphere. The Pt surface was mo deled as constan in thermostat, similar to Maruyama and Kimura (1999). The wall was heated to two temperatures of 150 K and 300 K respectively from an initial equilibrium temperature of 110 K. In the former case a gradual evaporation of th e liquid film occurred, while in the latter the Leidenfrost phenomenon was observed. In both the cases, a complete evaporation of the liquid film was obtained. The gradual evap oration is shown in figure 15. Wemhoff and Carey (2005) replaced the solid wall by deriving a wall potential to bind a liquid argon film to the imaginary wall, in order to re duce computational cost The wall p depends on the approximated values of Hamake r constant for metalmetal and fluidfluid interaction. The system pressure was adjusted to an equilibrium value using a flux boundary. 38 PAGE 39 Equilibrium properties were evaluated at different temperatures; film evaporation was not studied. In the literature available ( on molecular dynamics), there has been no mention of any observed nonevaporating thin film in simulations. Van der Waals forces play a central role in all phenomena involving intermolecular forces. Van der Waals force is constituted by three d types of forces: induction fo rce, orientation force and dispersion force. The dispersion contribution almost completely dominates the va n der Waals energy in the interaction betw two dissim istinct een ilar atoms of which one is nonpolar (Israelachvilli, 1994). Dispersion forces act betwee the the n part n all atoms and molecules, are always present and contribute to events like adhesion, surface tension, physical adsorption, wetting, liquid and thin films, etc. These forces may be repulsive or attractive, and may also tend to mutually align or orient the atoms between which they are acting. Hamaker constant is introduced to relate the interfacial free energies to dielectric properties of the materials and to use the surface excess convention for thin films. Hamaker (1937) performed an analysis where he at tributed the adhesion between particles to dispersion forces (also called London forces). The London potential is th e interaction betwee two atoms and by assuming additivity, he obtained the Van der Waals interaction energy between bodies consisting of many molecule s by integrating Londons pair potential of dispersion interaction. All of the Hamakers equati ons can be split into a purely geometrical and a constant A (widely known as Hamaker constant) which depends on the nature of interaction between the two bodies. The Hamaker cons tant is a very helpful parameter and can be used to calculate the van der Waals force of in teraction, dispersion energy, adhesive pressure, etc. It can also determine the magnitude of spreading of a film on different surfaces (Israelachvilli, 1994). A larger va lue of the Hamaker constant will result in a wetting film 39 PAGE 40 whereas a comparatively smaller value leads to a nonwetting film. The Hamaker constant A is defined as (Wayner, 1999): 2 12AC (1.1) where and 12 are the number of atoms per unit volume in the two volumes and C is a coefficient describing the interaction between two atoms. Nonevaporating thin film important role in the microlayer region of nucle ate boiling and meniscus evaporation, and thus the knowledge of Hamaker constant is essential. For example, numerical simulation of nucleate boiling requires the valu e of Ha s play an maker constant to determin e the nonevaporati ng film thickness (Dhir, 2001). Not ma hitz 9). ny Hamaker constant values for different fluidsolid combinations are available in literature. Hamaker constant values can be evaluated using Derjaguin and Lifs theories, but serve as approximations at best for nonretarded van der Waals forces. A more exact procedure would be to obtain it experiment ally using an atomic force microscope (AFM) and then determining its value from measured forcedistance curves (Butt, 2005). The Hamaker constants obtained will thus be limited to those particular se t of pressure, temperature and materials for which the AFM was used restric ting the usage of the va lues. Thus, additional development in the field of computational mo lecular dynamics to emphasize interfacial mass transfer is required to have an impact on changeofphase heat transfer systems (Wayner, 199 40 PAGE 41 Figure 15. Gradual and complete film ev aporation simulated by Yi et al. (2002) Capillary Pressures at Nanoscale The properties of materials at nanoscale are di fferent from those at macroscale due to the lengthscale effects like high surfacetovolume ratios and different dominating physical mechanisms acting at nanoscale. Capillary forces are often present in small mechanical systems and can become very significant at nanoscale. Fo r example, capillary bridges form when water condenses (even under low relative humidity) at the contact spots during the contact of two solid 41 PAGE 42 bodies in air. At nanoscale thes e capillary bridges can usually be under high negative pressures. The phase diagram of water, shown in Fig. 16, shows the stable, metastable and unstable regions (Angell, 1988). In the top righthand corner water is stable as a liquid. To the left of Tb, the locus of boiling points, water can exist in a metastable state as a superheated liquid in the absence of nucleation centers. To the left of Th, ordinary density fluctu ations are sufficient to cause prompt spontaneous cavitation resulting in homogeneous boiling. Similar arguments hold below Tm, the locus of melting points. Ts is the spinodal of speculated stability limits of water. In the past, the region between Th and Ts was experimentally accessibl e as water would cavitate to the left of Th at macroscale. However, it has been recently possible to have water exist at highly metastable states at nanoscale. When pressure is below the liquidvapor transition line of water phase diagram, vapor is the most stable stat e although the transition involves energy barriers. Random fluctuations are large e nough to overcome these barriers at macroscale, but at nanoscale the barriers are large compared to the scale of th e system and thus metastable state can exist for long intervals of time. An interest ing example of metastable state is stretched water, i.e., water under tensile stress of ne gative pressure. Pressure and surfac e tension affect the formation of bubbles inside a liquid which act towards expa nding and collapsing the bubble respectively. Formation of a cavity of volume V requires work equal to PV where P is the pressure of the liquid. Formation of liquidvapor interface boun ding the bubble requires work equal to LVA where LV is the interfacial tension and A is the area of the interface. The work required to fill the bubble reversibly with vapor of pressure Pv is negative and equal to PvV Value of Pv is negligible to P and is neglected (Fisher, 1948). Hence, th e net work associated with reversible formation of a spherical vapor bubble of radius r is: 42 PAGE 43 23443LV LVWAPV rr P Thus, for maximum work there is a cr itical radius of the cavitation Rc defined as: 2 L V cR P If the radius of a bubble is greater than Rc, it can grow unrestricted. No cavitation can occur if the size of the capillary is smaller than Rc. Thus, the liquid can be in a metastable state to the left of Th line (in Fig. 16) and the pressu re inside it may be negative. Briggs (1955) heated water in a thinwalled capillary tube open to atmosphere. The author heated the tubes upto 267oC for about 5 sec before explos ion occurred, and concluded that during the time before explosion occurs, the wate r must be under an internal negative pressure. When the temperature is raised to a point where the forces of cohesion no longer equal or exceed the increasing negative pressure, the system expl odes. Thus, water was able to withstand an internal negative pressure of over 51 atm. Yang et al. (2008) measured the pulloff force between a Si AFM tip and Si samples, as shown in Fig. 17, in UHV and in air under differe nt relative humidity values. As there is no water condensation and there no cap illary bridges in UHV, subtra cting the pulloff force in UHV from that in air results in the capillary force values. Very high negative pressures, down to 160 MPa, were obtained in the capillary bridge. The authors mention that no vapor bubble nucleation occurred in the experiments since the size of the capillary bridges was smaller than the critical cavitation radius. Tas et al. (2003) created water plugs by f illing water in a hydrophilic silicon oxide nanochannel of approximate height of 100 nm. Th e authors found out that the tensile capillary forces were so strong on this sc ale that the water plugs were at negative pressure. This negative 43 PAGE 44 pressure induced bending of the thin nanocha nnel capping layer which results in a visible curvature of the liquid meniscus as shown in Fig. 18 The amount of bending of the channel capping was calculated and the negative pressure was found out to be 17 10 bar. The absence of cavitation was explained by the fact that the cr itical radius for cavitation is comparable to the nanochannel height. Tyree (2003) discussed the transport system that drives sap ascent from soil to leaves in plants and trees. For very tall trees like redwoods ( Sequoia sempervirens ) which can grow over 100 m high, H. H. Dixon (1986) prop osed that a pulling force was generated at the evaporative surface of leaves and that this force was tr ansmitted downwards through water columns under tension to lift water from the roots. This theo ry, called the cohesiontens ion theory, assumes both adhesion of water to the conduit walls and cohesi on of water molecules to each other. Plants transport nearly pure water in the xylem conduits at negative fluid pressures of the order of 1 to 10 MPa. Thus, the cohesiontension theory propos es that water is transported in metastable state. The previous mentioned experimental tech niques have experimenta lly shown the existence of water in metastable state, which can be a major step towards confirming the principles of cohesiontension theory. Another important consequence of the water pres sure reduction in capi llary bridges is the possibility of boiling at temperatures lower than saturation temperatures at corresponding vapor pressure. If the capillary is larg e enough for cavitation to occur, the pressure inside the capillary bridge will be at a reduced pressure compared with the ambient and the boiling point will be depressed. Figure 19 shows the reported obser vation of boiling in a meniscus between the cantilever and a flat Si sample under laser irradiation th at provided a temperature rise of only up to 45oC (Nosonovsky and Bhushan, 2008). 44 PAGE 45 Figure 16. Phase diagram of water (Angell, 1988) Figure 17. Schematic of water phase diagram showing a metastable capillary bridge. The capillary bridge under tension acts like a wate r tube sealed with a piston, which leads to negative pressure. (Yang et al., 2008) 45 PAGE 46 Figure 18. Optical micrograph of a water plug in a 100 nm wide nanochannel (Tas et al., 2003) Figure 19. (a) Optical image (b) Schematic of water meniscus between a cantilever and a flat surface (Nosonovsky and Bhushan, 2008) Research Objectives In spite of the numerous numerical and experime ntal advances and decades of research in the field of multiphase physics, the moving contact line problem observed in a) bubble nucleation and growth and b) droplet spreading and impact has still no t been fundamentally understood due to the complex subprocesses invol ved which work at various length and time scales simultaneously. In numeri cal simulations, the classical noslip boundary condition leads to singularity at the contact line and the validity of the alternative ad hoc models proposed is far 46 PAGE 47 from being established. Exceptionally highreso lution experimental hardware, that is not currently available, is required to study the motion of the contact line at the molecular level. Researchers have agreed on the fact that molecu lar dynamics simulations of the moving contact line can lead to gaining d eeper insight into the funda mentals of the problem. The most accurate model can be developed if the problem is solved using molecular dynamics for the entire domain of interest. But cu rrently the roadblock to this method of solution is that simulating a bubble on the millimeter scale in molecular dynamics would require more than billions of atoms which is not currently possible with the available computing hardware. The contact line region is of the or der of a few molecular diameters ( nanoscale ), whereas the bubble grows on the millimeter scale. Thus, the best possible solution would be to couple the molecular dynamics with the continuum mechanic s. Two feasible approaches, with respect to nucleate boiling, are: 1. The problem is formulated by applying th e molecular dynamics simulation to the nanoscale region near the vicin ity of the contact line and then coupling that with the continuum mechanics for the remainder of the domain. The coupled system is then solved simultaneously. 2. Simulate the nanoscale vicinity around the moving contact line with molecular dynamics in conditions that mimic nucleate boiling and extract useful inform ation/models to be used at larger lengthtime scales. The first approach, also called the hybrid at omisticcontinuum model, has been implemented by a few researchers [OConnell and Thompson ( 1995), Hadjiconstantinou (1999), Flekkoy et al. (2000), Wagner et al. (2002), Buscalioni a nd Coveney (2004), Ren and E (2005), Nie et al. (2006), Liu et al. (2007)]. The difficulty with this approach is that it would be still very 47 PAGE 48 computationally demanding to simulate a bubble on the millimeter scale and would definitely require simulating billions of atoms on the molecula r level. Decoupling the timescales is also an unresolved issue since time steps in molecular dynami cs is of the order of picoseconds, while in continuum mechanics it is in micro to milliseconds. The work done so far by researchers using the hybrid method has been based on nanolength scales only for the entire domain, with the aim of showing the validity of the hybrid approach by comparing the results obtained with fully molecular ones. Thus, solving problems involving coupling nano and micro/millimeter length and time scales using the first approach has yet to be becomes a success. The second approach, which is also known as the information passing scheme (Chen and Fish, 2006), is adopted in this propo sal. The objective of this appr oach is to simulate the local nanoscale vicinity around the moving contact line by molecular dynamics, then study and analyze the molecular dynamics results to extrac t useful information/models which can be used in larger lengthtime scale simulations. Thus, the basic idea here is to extract useful information from molecular dynamics simulation to obtain a deeper understanding of the physics occurring which cannot be captured in experiment s due to the length scale involved. The nonevaporating thin film is the main focus in this work. Its thickness is expected to be in the size of nanometers, and thus is stud ied here using molecular dynamics simulations. The specific goals of this proposal are to: 1) simulate thin film evaporation using mol ecular dynamics, study the physics obtained and evaluate useful properties using statistical mechanics 2) simulate droplet impingement and spreadi ng on surfaces of varying wettability and evaluate moving contact a ngle variation with time 48 PAGE 49 3) use molecular dynamics to model an evaporati ng meniscus in order to gain a fundamental understanding of the physics occu rring at the threephase interline region coupled with the nonevaporating region Chapter 2 outlines the molecular dynamics theory and the algorithms used along with statistical sampling approach. A novel thermal wall model is introduced and validated to simulate heat transfer between wall and fluid. Chapter 3 shows the results obtained from simulation of thin film evaporation and the occu rrence of nonevaporating thin film. A method to calculate the Hamaker constant from molecular dynamics simulations is outlined, which has not been attempted in literature befo re. The transient film evaporati on is extended to seven cases in order to study the effect of na nochannel height, initial film thickness and wall temperature on the formation of nonevaporating film. Chapter 4 studie s the impact of a nanodroplet on a wall with varying surface wettability. Variation of contact angle with time is evaluated. The wall temperature is also changed to observe the Leidenfrost phenomenon. Chapter 5 simulates the evaporation of a meniscus at nanoscale. Two cas es, with partial and complete wall heating, are studied. Heat and mass transfer characteristics an d pressure variation in the nonevaporating and interline regions are studied. A curvefitted me niscus boundary condition is also developed. The dissertation ends with conclusions and future work in Chapter 6. 49 PAGE 50 CHAPTER 2 THEORY AND APPROACH Molecular Dynamics Simulation The most rudimentary model for a substance capab le of existing in any of the three states of matter solid, liquid and gas is based on atoms interacting with one another, without considering their quantum orig ins. The interactions are re sponsible for providing the two principal features of an interatomic force: a) resistance to compression (re pulsion at close range), and b) binding atoms together is liquid and solid states (attrac tion over a range of separation). Potential functions exhibiting these characteri stics can have many forms, and when chosen properly, provides accurate models for real s ubstances. Usually this function comes from experimental data or quantum mechanics calcul ations. The interaction force between the atoms can be derived from the potential function U by the following relation: FU Summing up all the forces acting one atom by all other atoms yields the total force acting on that atom. The best known potential function is the LennardJones (LJ) potential. This func tion was originally proposed for liquid argon. For a pair of atoms i and j located at and irjr the potential energy as defined by the LJ potential is: 1264 () 0 ijcut ij ij LJij ijcutrr rr Ur rr (2.1) where ijijrrr The parameter governs the strength of the interaction, and defines a length scale. The interacti on repels at close range, then attracts and is eventually cut off at some separation rcut (known as the cutoff radius) to reduce computational cost while serving as a reasonable approximation and still providing good accuracy to the simulations. The LJ potential 50 PAGE 51 is the most widely used potential for molecula r simulations, and provides accurate representation for the properties of argon (Sadus, 1999). In order to smoothly truncate the LJ potential and remove the discontinuity at the cutoff distance, a few modifications have been suggested in literature. The modified LJ model proposed by St oddard and Ford (1973) has been used in this work for both argonargon and argonplatinum interactions: 12 6 2 12 6 126()4 63 74MLJ cut cutcut cut cutr Ur rrrrrrr (2.2) The above potential form is employed for bot h argonargon and argonplati num atom interaction with the following values [6]: 103.410,ArArm 211.6710,ArAr J The mass of an argon atom is 103.08510,ArPtm 210.89410ArPtJ 266.6910 kg, 120.95BkK and time step is 5 fs. The and parameters in Lennard Jones (LJ) potential are obtained by analysis of the experimental da ta for second virial coefficients, transport coefficients, constants of the critical, melting or boiling points, or other properties (Hirschfelder et al. 1964). The author presents the equations (in the sections as mentioned above), which when fitted with the experimental data, give the LJ parameters. The experimental values of viscosity of argon are listed by Johnston and Grilly (1948), an d isotherm data (for virial coefficient evaluation) is listed by Michels et al. (1949). Recently, other met hods have been introduced to obtain the parameters, like using abinitio qua ntum mechanical and molecular mechanical potentials (Marti n et al. 2002). Force Calculation Algorithm In a typical molecular dynamics simulation, as much as 95% of the computing time is spend in examining the complete set of N(N1)/2 pairs, identifying thos e pairs separated by less than the cutoff distance and computing the forces for this subset (Allen and Tildesley, 1989). 51 PAGE 52 Hence, the algorithm used for force evaluation becomes very important. For simulations regarding liquidvapor phase chan ge (where density continuously keeps changing), the two most widely used algorithms are: Brute Force Method In this method, to compute the total force on an atom i the force acting on i due to all other atoms 1 to N has to be calcu lated. Thus, all pairs of atoms are examined. This is simplest to implement, but very inefficient when the intera ction cutoff radius is small compared to the linear size of the domain. The amount of computation needed grows as order of N2 (or N(N1)/2 by taking advantage of Newtons third law). Linked List Method This algorithm is a cellbased method and invol ves data organization, and is particularly efficient for simulations involving a large number of molecules (N 1000) (Allen and Tildesley, 1989). The computational domain is divided into cells with each side of the cell being greater than or equal to rcut. At any time step, the atoms are binne d into their appropriate cells. Two arrays, headofchain array (HEAD) and linkedlist array (LIST), are created during the binning process. The HEAD array contains the identificati on number of one of the atoms binned into that cell, and this number acts as a pointer to the LI ST array which contains the number of all the atoms in that cell in a sequentia l order. Thus, this method redu ces the task of finding neighbors of a given atom to checking in 27 cells the cell the atom is in and 26 neighboring cells. Particular care has to be taken in this method when dealing with vapor or liquid regions, as some of the cells might not contain any atom. For a three dimensional system, there are on average Nc = N/M2 atoms in each cell. The computing power is of the order of 27NNc (or 13.5 NNc). This algorithm results in faster computation speeds and is implemented in this work. 52 PAGE 53 Integrator Method The equation of motion needs to be integrated in order to obtain the posi tions and velocities of the atoms at every time step. Thus, after calculating the force acting on each atom i from the potential function, Newtons second law is used to determine the acceleration: 2 2 ii iii idrF Fma dtm (2.3) where is the force acting on atom i, and and F ma are its mass and acceleration respectively. Using the acceleration, the position and velocities are obtained via the integrator method for time t + t The integrator method used in this work is the VelocityVerlet method: 2() ()()() 2 tat rttrttvt (2.4a) () (2)() 2 tat vttvt (2.4b) ( ()(2) 2 tatt vttvtt ) (2.4c) The advantage of this method is that the calculati on of velocities is in phase with that of the positions. Also, this method calculates the velocities more accurately than the other common Verlet and leapfrog methods (Sadus, 1999). Thus after the new positions and velocities are obtained at time t + t the whole process is repeated agai n i.e. force computation from linked cell algorithm, using Newtons second law and then integrator method to advance the system in time. NonDimensional Units The use of dimensionless parameters in molecular dynamics is very helpful due to several reasons (Rapaport, 2004). The values of various parameters are extremely small as they are associated to the atomic scale. Implementation of dimensionless parameters makes it possible to 53 PAGE 54 work with values which are not far from unity. Also, these units remove any risk of facing computer hardware related problem of not be ing able to represent such small values. Dimensionless units simplify the equations of mo tion as the parameters defining the model are reduced. Also, it helps in scaling of certain computed properties and a single model can be used to describe a whole class of problems. The dime nsionless parameters, denoted with an asterisk (*), used here are defined in Table 21. Table 21. Dimensionless Parameters Parameter Dimensionless form Energy ArU U Mass refm m m Distance Arr r Velocity 1/2 ref Arm vv Time Ar Arreft t m Temperature **1 ;. (dim)B ii ArkT TT N vv Density *3 Pressure 3 *P P The force acting on an atom i due to interactions for argon argon and argonplatinum can be deduced (from derivation of potential function) to the following equations: ** 1 4* 8 1 4* ()480.5 0.5ii ji jc u tc u t ji8 i j A rgonArgoninteraction rrrrr r (2.5) 54 PAGE 55 12 6 12 6 ** 1 4 8 1 4 ()48 0.5 0.5ArPt ArPt ArPt ArPt ArPt ii j i j c u t ji Ar Ar Ar Ar ArArgonPlatinuminteraction rrr r 8 c u t i jr r (2.6) Boundary Conditions In the present work, the computational domains are in the form of a cuboid. This imposes a total of six boundaries, two each in the x, y a nd z directions. Thus, the same number of boundary conditions is required. Three boundary condi tions are used in this work depending on the problem simulated. Periodic Boundary Condition (PBC) A system that is bounded but free of physical walls can be constructed by using this boundary condition. This is equivalent to consideri ng an infinite, spacefilling array of identical copies of the simulation region. An atom that leaves the simulation re gion through a particular boundary face immediately reente rs the region through the opposite face. Also, atoms lying within the cutoff distance of a boundary intera ct with atoms near th e opposite boundary. Thus, periodic boundaries have to be in multiples of two. Mirror Boundary Condition This boundary condition simulates an imaginar y adiabatic wall. The collision of atom with this wall is perfectly elastic; momentum and energy are conserved. If an atom, which is close to the wall at time t crosses this boundary at atom t + t it will be placed inside the boundary as if it has been reflected by a mirror. The atom will be assigned the same position as it was at time t with the velocity directed away from the wall. 55 PAGE 56 Thermal Wall Boundary Condition A novel fluidwall thermal equ ilibrium model to numerically simulate heat transfer between wall and fluid atoms in molecular dynamics is introduced in this work. The details of this proposed model, and its validat ion is discussed in later section. Initial Structure The domain in the simulations consists of argon fluid and platinum wall(s). Platinum atoms, arranged in fcc (111) structure as shown in figure 21, constitute the platinum wall with the interatomic distance as defined in literature (Yi et al., 2002). As will be discussed later, the platinum atoms remain stationary throughout the simulation. For argon fluid, the distance between the argon atoms is calculated as follows: knowing the volume occupied by liquid (or vapor) argon, the saturation dens ity at the equilibration temper ature (from argon thermodynamic tables) is used to obtain the total mass of li quid (or vapor) argon, which divided by the mass of an argon atom gives the total number of liquid (o r vapor) argon atoms. The distance between the atoms is then obtained assuming th ey are arranged in a simple cubic structure. The argon atoms are arranged by this interatomic distance thus de fining the initial structure of the domain. The argon atoms are assigned initial random velocities such that the average temperature of the domain is the desired temperature. No initial ac celeration is given to th e argon atoms. The results of a simulation of ample duration are not sensitive to the initial state. 56 PAGE 57 Figure 21. Arrangement of Pt wall atoms in fcc(111) structure; (a) 3D view of the first three layers, (b) XY projection (t op view) of the wall atoms Simulation Stages All the simulations in this work are given an initial temperature. A simulation procedure in implemented to make the argon fluid reach equilibrium at the desired temperature from the initial structure outlined in the previous section. The simu lation process is divided into two stages to achieve the above. In the first stage, a velocityscaling period is used where the velocity of each argon atom is scaled at every time step to maintain the system at a constant temperature (Sadus, 1999). The velocity of each component of an atom i is scaled from it actual velocity to a new velocity in the following way: new d ixix aT vv T ; new d iyiy aT vv T ; new d iziz aT vv T (2.7) where Td and Ta are the desired and actual system temperatures respec tively. This is followed by a second stage, called the equilibration period, wher e the velocity scaling is turned off and the atoms are allowed to move naturally. This is done to remove the artificial motion which the atoms were subjected to in the first stage. The system temperature fluctuates about the desired temperature in this period. A third stage can foll ow in the simulation pro cedure if an external 57 PAGE 58 wall is present in the computationa l domain. In this stage, named the heating period, heat transfer from the wall occurs and the actual physic s of the problem is observed and studied. Statistical Sampling Statistical mechanics provid es the link between therm odynamics and random atomic behavior. Thermodynamic properties such as pressure, temperature, energy, etc can be calculated by averaging the atoms positions, velocities, inte ratomic potential and forces. In this work, several properties have been evaluated by using statistical mechanics equations. Density With phase change being the main focus in this work, calculati on of density becomes quite important. The domain is divided into equa l slices in zdirection. The number of atoms in each slice is counted to compute the dens ity. Hence, if the slice thickness is z and Ni is the number of atoms in slice i then the instantaneous density of that slice can be found out as: ()i i xymN z llz (2.8) where lx and ly are the x and y dimensions of the doma in respectively, and m is the mass of the atom. The average density can be computed by timeaveraging the instantaneous density as follows: ()ii z The average density can be plotted al ong the height (zd irection) of the domain to see the density gradient at different times. Number density calculates the average number of atoms in a slice, and can be found out using the same equation as above but without the m term. Energy The total energy E of the domain can be obtained by adding up the kinetic Ekin and potential Epot energies. The kinetic energy can be com puted from the individual momenta of the atoms, while the potential energy is the sum of interatomic interactions. 58 PAGE 59 2 11 2N kin i i E mv (2.9) 12 6 2 12 6 12646 37pot cut cutcut cut cutr E rrrrrrr 4 (2.10) kinpot E EE (2.11) Temperature The equipartition principle is used to calculate the temperature. For N atoms, each with three degrees of freedom, the instanta neous temperature of the domain is: 22 1 33kin i i BBE T NkNk m v (2.12) Temperature gradient along the height of the domain can be computed by dividing the height into equal slices, and co mputing the average temperature of e ach slice, similar to the calculation of density as mentioned before. Heat Flux The microscopic heat flux Q for a onecomponent system of N particles interacting via pair potentials in a control volume V can be written as follows (de Andrade and Stassen, 2004): 1kPwQQQQ V (2.13) where, 2 11 2N k iQmi ivv is the kinetic energy transported by mass flow 11 () 2NN P iji ijiQU rv is the flux of potential energy due to mass flow 59 PAGE 60 1() 1 2NN ij wi j iji ijdUr Qrr dr i j i i jvr is directed along the intermol ecular separation by the force displacing the pair of particles i and j by rij, and accounts for the coupling of this work term to the mass flow Thermal Conductivity From the definition of temperature, the temperature gradient dTdzcan be computed when the system reaches equilibrium. The heat flux can be computed as described above. Thus, the thermal conductivity can be calculated as follows: simulkQdTdz (2.14) Pressure Weng et al. (2000) suggested an improved method to estimate local pressure distributions by combining the work by Nijmeijer et al. (1988) and Daiguji and Hihara (1999) with Kirkwood and Buffs theory (Rowlinson and Widom, 1982). This method takes into account that an intermolecular force contributes to local pressure in each plane between the two molecules. The domain is divided into equal slices in the desi red direction. The local normal pressure component pN and tangential pressure component pT are expressed as: 2 ,1 ()() ()ij k NB ij slab ijz pslabknkkT Urf Vr i j k i j (2.15) 22 ,1 1 2 ()() ()ijij k TB i ij slab ijxy pslabknkkT Urf Vr j k i j (2.16) where n(k) is the number density in slab k Vsl is the volume of slab k and fk,ij is kijij L z. The first terms in both the equations are the contribu tion from kinetic motion of molecules while the 60 PAGE 61 second terms are the contribution from the interatomic force. The length Lk,ij is defined as the size of the interval in which the in termolecular force between molecules i and j is effective in the slab k For example, consider the force between atoms i and j as shown in figure 22. If neither atom is in a slab, such as slab 3, L3,ij = Lslab; if one atom is in a slab, such as slab 1, L1,ij = L1. If both atoms are in the same slab, such as atoms p and q, L7,ij = ijz. Surface tension, defined as the difference between the local normal and tange ntial pressure components, can be expressed as: 01 2zL NT p pdz (2.17) Slab 1 Slab 2 Slab 3 Slab 4 Slab 5 Slab 6 Slab 7 zijrijL1Lslabi pq j Figure 22. Sketch for pressu re contribution of atoms i j Pressure Verification The method described in above section is used to evaluate pressure distributions. For verification, pressure and surface tension are calcu lated in a computationa l domain consisting of a liquid layer surrounded by vapor, as shown in figure 23, for si x different temperature cases. The domain parameters for all the cases are lis ted in table 22. Periodic boundary condition is 61 PAGE 62 applied to all six boundaries. The domain is divi ded into slabs of thickness 0.34 nm in the zdirection. The timestep is taken as 5 fs. Velocity scaling is done for the initial 500 ps and then the system is allowed to equilibrate for an addi tional 1000 ps. The Statistica l values are averaged over the final 100,000 timesteps. Tempratureden sity, pressuredensity and surface tensionpressure plots are shown in fi gures 24, 25 and 26 respective ly, in which the results from simulation are compared to experimental values from thermodynamic tables for argon (NIST webbook). It can be seen that the values match very well. Table 22. Simulation parameters for pressure verification Temperature (K) lx (nm) ly (nm) lz (nm) Liquid layer thickness (nm) Number of argon atoms 90 6.20 6.20 14.63 6.20 4949 100 6.20 6.20 14.63 6.20 4949 110 6.20 6.20 14.63 6.20 4949 120 6.20 6.20 14.63 6.20 4949 130 4.79 4.79 13.33 4.79 2213 140 6.20 6.20 12.51 4.08 3215 62 PAGE 63 Figure 23. Simulation domain at an intermed iate timestep for pressure verification Figure 24. Temperaturedensity plot for liquid slab surrounded by vapor atoms 63 PAGE 64 Figure 25. Pressuredensity plot for liquid slab surrounded by vapor atoms Figure 26. Surface tensionpre ssure plot for liquid slab surrounded by vapor atoms 64 PAGE 65 Distinguishing Between Liquid At oms and Vapor Atoms for Argon An important part of the simulati on is to define a criterion to la bel a fluid atom as liquid or vapor. The argon atoms in the liquid phase are differentiated from the argon atoms in the vapor phase based on the minimum number of neighboring atoms within a certain radius. Thus in order to obtain this threshold number density and thres hold radius for a saturation temperature of 90 K, a computational domain of size 4.72 nm 4.72 nm 4.42 nm is constructed with periodic boundary conditions in x, y and z directions. The domain consists of 2028 argon atoms in saturated liquid state at 90 K. Velocityscali ng is done for 0100 ps and then the system is allowed to equilibrate for an additional 1000 ps. The minimum threshold number of all the atoms is determined within a radius (ranging from 0.41 nm to 0.59 nm at steps of 0.01 nm at every 0.25 ps (i.e. total of 4000 times) during the equilibration period. The averag e of these values is plotted against the radius in figure 27. The plot also shows the number of times (in %) when the minimum threshold number at a certain time step falls below the average. Threshold radius of 0.53 nm is chosen as the criterion and the averag e threshold number is rounded off to the lower integer. Thus, if an argon atom has 7 or more neighboring atoms within a radius of 0.53 nm it is said to be a liquid atom or else a vapor atom. 65 PAGE 66 Figure 27. Plot to obtain the th reshold number density and thres hold radius for saturated liquid argon at 90 K FluidWall Thermal Equilibrium Model The heat transfer from the heated wall to the fluid is simulated with the following wallfluid interaction mechanism. The heated wallfluid thermal energy transfer interaction is constituted of two parts: 1) interaction potential model, and 2) proposed fluidwall thermal equilibrium model (to numerically simulate h eat transfer between wall and fluid atoms). The interaction potential model, using the LJ potenti al, is a standard method as it has been used by many researchers for argonplati num interaction in molecular dyna mics approach. The fluidwall thermal equilibrium model is presented in details below. Figure 28 shows the schematic of the fluidheated wall thermal equilibrium model. The Pt atoms constituting the heated wall are assumed to be fixed in space, and do not move at any time during the simulation as for a solid, their range of movement is relatively much restricted than those of the fluid atoms. In addition, it is al so aided by the fact that the mass of a Pt atom is 66 PAGE 67 roughly five times that of an argon atom. The flui dwall thermal energy exchange is simulated as follows: the Pt wall is at a specified constant temperature, and an argon atom reaches this temperature (i.e. the atom is assigned a corres ponding velocity by instantaneous velocity rescaling; direction of the atom is not changed) by coming in thermal equilibrium with the temperature of the wall The force due to the wall atoms on th e fluid atoms is either attractive or repulsive depending on the distance between the two due to the nature of the LennardJones potential. A distance rcr less than rcut, is defined as the neutra l position of wall influence. Accordingly, any fluid atom on the z = rcr, as shown in figure 28, would experience a zero force from the first layer wall atom that is located at the shortest distance from this fluid atom. It is noted that this fluid atom would experience onl y attractive forces from other wall atoms as the distance between this fluid atom and other wall at oms are all larger than that shortest distance Therefore any fluid atom will have to get below z = rcr (i.e. in region 3) to have any chance of experiencing a repulsive force from the wall atoms. Based on the above scenarios, the fluid atoms can be divided into three groups based on th eir locations in the respective regions near the wall as shown in figure 28. Region 1. z > rcut : The fluid atoms which are above th e cutoff radius do not experience any force from the wall atoms. Region 2. z < rcut & z > rcr : In this region, the fluid atoms only experience an attractive force by the wall atoms due to the nature of the LJ potential. Region 3. z < rcr : In this region, the fluid atoms can experience either an attr active force or a repulsive force due to the wall atoms. 67 PAGE 68 z = 0 z = rcrz = rcut Fixed wall atoms Fluid atoms experience attractive or repulsive force by wall atoms Fluid atoms experience only attractive force by wall atoms Fluid atoms experience no force from wall atoms Figure 28. Fluidwall th ermal equilibrium model The physical mechanism of a fluid receiving heat from a solid wall is composed of a twostage force interaction that result s in an energy transfer from th e hot wall to the fluid atoms. During the first stage, the fluid atom is accelerat ed towards the wall atoms by the attractive force. The fluid atom would then reach a zero force point before being bounced back by a repulsive force. If the fluid atom is heated up, the re bounded fluid atom would possess more kinetic energy than its original kinetic energy be fore the interaction. Th erefore, a fluid atom is considered to be heated up to the wall temperature only when: (a) the atom is in region 3 as it can only experience a repulsive force from wall atoms in this regi on and the repulsive force is the indispensible second stage of the heat excha nge interaction mentioned above (b) the magnitude of the zcomponent force due to the wall atoms zfluidwallris greater than the zcomponent of the attractive/repulsive force due to the neighboring fluid atoms, zfluidfluidr. This condition is based on the concept that the fluid atom must be prim arily dominated by the wall force in order to satisfy the requirement that it is receiving energy from the wall atoms. These proposed criteria allow us to numerically simulate the energy exch ange due to collisions between the fluid atoms and the wall. The criteria are summarized as: 68 PAGE 69 fluid atom is in Region 3 Thus, if AND fluid atom is heated up at the wall temperature.zfluidwallzfluidfluidrr At each time step, the fluid atoms which meet the abov e criteria are identified and randomly assigned temperatures within K of the wall temperature such that the mean value of the temperatures for these atoms equals the wall temperature. The stepchange in wall temperature incorporated in this model theoretically corresponds to a wall with an infinite thermal conductivity, and the heating time peri od of the wall is neglected. This was again taken as a reasona ble assumption for highly conductive metallic materials since the heating time period for threemolecular layers of platinum would be in the order of 110 pico seconds or less. Thus, it is reasonable to neglect the wall heating time and not to affect the physical conclusions. Maruyama a nd Kimura (1999) used Langevin thermostat to increase the temperature of Pt wa ll atoms by 10 K; they conclude that the solid wall temperatures quickly respond to the temperature change, which strengthens our assumption of neglecting the heating time period for platinum. Regarding th e role of electrons in the platinum heater, Yamamoto et al. (2000) studied the structure of a Pt(111)/dipolar liquid interface by fully selfconsistent combination of firstprinciples calculation based on quantum mechanics for the metal and reference hypernettedchain (R HNC) theory for the liquid. The metal region is represented by a supercell with repeated slab s each of which consists of elev en Pt(111) layers. The liquid is modeled as hard spheres of diameter 0.28 nm embedded with point dipoles. For the liquid molecules, the metal acts as a hard wall. It is assumed that the two phases interact only through electrostatic fields and there is no charge or mass transfer between the two phases. The electrostatic density profile for the metal, density and orientational struct ure of liquid molecules, and electrostatic potential acros s the interface and evaluated and discussed. The authors compare 69 PAGE 70 the density and orientatio n profiles in liquid for metal wall and inert wall (where the liquid feels no electric field). They see that a layer of liqui d molecules, which is orientationally ordered, is formed near the metal surface, but extends to about only three molecular diameters from the surface. The density and orientational structure de cays very rapidly with distance from the metal surface. The authors mention that their result is in good accord with the recent experimental observations. Model Validation The proposed model is validated by evaluati ng both statedependent and timedependent properties. Molecular dynamics simulations are conducted by placing liquid argon between two platinum walls. The validation is divided into two sections. Internal energy and thermal conductivity The computational domain is a cubo id with dimensions of 4.794 nm 4.794 nm 8.238 nm. The domain comprises of 3718 argon atoms at liquid saturation density of 100 K. The Pt walls are placed at z = 0 and z = 8.238 nm and is made up of a total of 2101 Pt atoms. The fluidwall boundary condition suggested in this paper is appl ied here. The initial system temperature is 100 K. The domain, at an intermediate time step, is shown in figure 29. The time step is 5 fs. The velocityscaling period is from 0500 ps while the equilibration period is from 5001000 ps. The domain reaches equilibrium during the equili bration period, which is then followed by the heating period. At the start of the heating pe riod (10002000 ps), the wall temperature is stepchanged to a different temperature and the he atingcooling process of the liquid argon is observed. Two different cases are simulated: (CASE a) lower wall temperature is increased to 110 K while the upper wall temperature is decr eased to 90 K, and (CASE b) lower wall temperature is increased to 130 K a nd the upper wall decreased to 90 K. 70 PAGE 71 Figure 29. Computational domain at an intermed iate timestep for validation of proposed wall heat transfer model The average energy and temperature of the syst em is evaluated at the initial and final states, and the plots are shown in figure 210 (scaling period is not shown). The system temperature for both the cases reaches a value clos e to the average temperature of the two walls. The values are listed in Table 23. The difference in internal energy calculated in the simulation is compared to that obtained from the experi mental thermodynamic properties of argon using Engineering Equation Solver which implements a high accuracy equation of state (Klein, 2007). The thermal conductivity of argon is also calculated in the simulations, as outlined in section 2.8.5. The temperature gradient dTdzis evaluated along the height of the channel and thermal conductivity is determined. These values are co mpared to the experimental values from thermodynamic properties as listed in Table 24. The experimental k values listed are for saturated liquid argon at 100 K and 110 K. The th ermal conductivity and change in internal energy values obtained match well with the experimental values for the two cases simulated, which can only happen if the proposed wallflui d boundary condition is nu merically correct. It should be noted here that the internal energy of the fluid is evaluated over the complete fluid domain (including the atoms in the vicinity of the wall). 71 PAGE 72 (a) (b) Figure 210. Variation of temperature and energy with time for validation of proposed wall heat transfer model (a) lower wall temperature is increased to 110 K while the upper wall temperature is decreased to 90 K; (b) lo wer wall temperature is increased to 130 K and the upper wall decreased to 90 K 72 PAGE 73 Table 23. Change in internal energy for valid ation of proposed wall heat transfer model Equilibration Period (avg. over 5001000 ps) Heating Period (avg. over 30005000 ps) Temp. (K) Internal Energy (kJ/mol) Temp. (K) Internal Energy (kJ/mol) usimul (kJ/mol) uEES (kJ/mol) % Error Case a 99.40 3.8546 99.54 3.8546 0 0 Case b 99.95 3.8430 107.83 3.6925 0.1505 0.155 2.90 Table 24. Thermal conductivity for validati on of proposed wall heat transfer model Heating Period (avg. over 30005000 ps) Temp. (K) Heat flux (MW/m2) Temp. gradient (K/nm) ksimul (W/m.K) kexp (W/m.K) % Error Case a 99.54 265.96 2.331 0.114097 0.10887 [at 100 K] 4.80 Case b 107.83 471.20 4.316 0.109175 0.09589 [at 110 K] 13.85 Heat equation Although the previous validation effort provides us with phys ical soundness of the proposed model, both the therma l conductivity and the change in internal energy are state dependent but time independent. Thus to furthe r validate our model, a similar simulation is performed where the temperature gradient in liquid argon, at different time intervals, is compared with that of the 1D analytical solution obtained from the heat equation. A 1D problem, as shown in figure 211, can be defined as: 73 PAGE 74 2 2, 1 TztTzt tz (2.18) Initial condition: (,0)iTztT Boundary conditions: (0,)lowTztT (,)uppTzhtT z = 0 z = h TlowTup T(t=0) = Ti Figure 211. 1D heat problem The solution to the above equations can be shown to be: 212 (,) 1 sinn t n H lowuplow upi lowi nzn TztTTT TTTTze Hn H (2.19) where H is the nanochannel height and is the thermal diffusivity of fluid. Equation 2.19 will be used to find the analytical solution of temperature gradient in liquid argon at different time intervals. Any of the previous works in literature have not attempted to validate their wallfluid heating/cooling models by compari ng temperature gradient in liquid to analytical solution of the heat equation. H eat conduction is conventionally tr eated as a diffusion process on the concept of local thermal equilibrium. The lo calequilibrium condition breaks down when the characteristic length L is smaller than a mechanistic lengt h, such as the mean free path. For 74 PAGE 75 example, heat conduction across a di electric thin film when the thickness is much smaller than the phonon mean free path or energy transport across a thin film whose thickness is much smaller than the mean free path, conventional Fouriers law cannot be directly applied. The hyperbolic heat equation, with no heat generation, is defined as: 2 2 21qT T tt T (2.20) where q is a kind of relaxation time similar to the average time between atomic collisions. The solution of the hyperbolic heat equation is in the form of a propagating wave whose amplitude decays exponentially as it travels. Zhang (2007) shows a sample calculation to compare the heat transfer between Fouriers heat equation and th e hyperbolic heat equati on for very small and large q For small q the parabolic heat equation incorrec tly predicts a continuous temperature distribution without any wavefront But as time passes on, the first order time derivative (i.e. the diffusion term) in the hyperbolic heat equation above dominates. After a sufficiently long time, usually 5 to 10 times q a local equilibrium is reestablished and the thermal field can be described by the parabolic heat equation. Also, at steady state, the hype rbolic and parabolic equations predict the same results. In the present work, the temperature increase of the wall occurs at 1000 ps. The time step of the MD simulation is taken as15510 s, which is around the same order as the average atomic collision time in liquid argon (1014 to 1015 s). Based on the above theory, the parabolic heat equation should be applicable after ~10 to 20 timesteps. The MD simulation is run till 2500 ps, i.e. for an additional 1500 ps (=300,000 timesteps) after the temperature of wall was increased. Thus, the pa rabolic heat equation results hold during the majority of the time period over which the temper ature gradients are evaluated and averaged, as it is much greater than q 75 PAGE 76 The velocityscaling period is from 0500 ps the equilibration period is from 5001000 ps followed by the heating period (10002500 ps) where the wall temperatures are stepchanged to desired temperatures. The computational domain is a cuboid with dimensions of 5.939 nm5.939 nm23.657 nm. The domain comprises of 16384 argon atoms at liquid saturation density of 100 K. The Pt walls are placed at z = 0 and z = 23.657 nm and are made up of a total of 3224 Pt atoms. The initial system temperature is 100 K. The time step is 5 fs. At the end of equilibration period (i.e. t=1000 ps), the upper wall temperature is increased to 110 K and the lower wall temperature is decreased to 90 K. Figure 212 shows the variation of average syst em temperature and energy with time. The system temperature fluctuates about a mean va lue of 100.43 K and the system mean energy is 0.041164 eV/atom, during the equilibration period. The average system temperature and energy evaluated in the final 500 ps is 99.60 K and 0.041374 eV/atom. Figure 212. Average temperature and energy variation with time of liquid argon time for validation of proposed wall heat transfer model 76 PAGE 77 In order to evaluate the temp erature gradient, the domain wa s divided into 70 equal slices along the height of the channel (i.e. in z directio n) and the average temper ature of each slice was evaluated at every timestep. The temperature gradient was then averaged over 100 ps (20000 timesteps) to obtain good statistical averages. Hen ce, if the temperature gradient is averaged from time t ps to t+100 ps, the statistical averag e can be said to hold at t+50 ps and the analytical solution is calculated at that timestep. The simu lation values are compared to the heat equations analytical solution (equation 2.19) in figure 213. As it can be see n, the results match quite well thus validating the proposed fluidwall therma l equilibrium model on a temporal basis. Liquid argon heating from both platinum walls An additional case was simulated by h eating liquid argon from both walls. The computational domain is similar to previous simulation: 5.939 nm 5.939 nm23.657 nm, 16384 argon atoms, and 3224 Pt atoms constituting the walls The initial system temperature is 100 K. The upper and lower wall temperatures are increa sed to 120 K at the end of equilibration period (i.e. t=1000 ps). The temperature gr adient is averaged over 100 ps and the simulation results are compared to analytical solution in figure 214. It can be seen that heat transf er occurs at a faster rate in the nanochannel than evaluated from heat equations analytical result. This occurrence can be attributed to the length sc ale of the current case as the heat equation is usually applied at a continuum scale. Hence, faster heating/cooling rate at nano length scale is not unexpected. The results compared well in the previous case (w ith heating/cooling at the walls) since higher heating rate at the upper wall is matched with higher cooling rate at the lower wall. It is fairly easy to imagine that heat transf er across a nanochannel of an even smaller height (say about 510 molecular diameters) would happen at a much highe r rate. Also from this conclusion it should be expected that with an increasing nanochanne l height, the e rror between molecular dynamics 77 PAGE 78 simulation result and analytical re sult would keep decreasing up to a particular channel height, beyond which the results would perfectly match. To verify this point of view, another simulation in conducted with increased nanocha nnel height as explained below. The dimensions of computational domain are 5.939 nm 5.939 nm44.357 nm. The height is about ~ 1.9 times the previous case. The domain now consists of 30720 argon atoms, and 3224 Pt atoms. The initial system temperature is 100 K. Similar to previous run, the upper and lower wall temperatures are increased to 120 K at the end of equilibration period (i.e. t=1000 ps). Figure 215 shows the comparison between simula tion and analytical results. It can be seen that the results match better than previous si mulation of shorter nanocha nnel (figure 213). In order to find the error, the average temperatur e of the domain is obtained from the temperature gradients of simulation and analyt ical results and plotted in fi gure 216 at intervals of 100 ps. The simulations were only run till 2500 ps due to computational limitations and final steady state is not obtained in either case. It can be noticed from the initial error (at 50 ps) that for a shorter channel height, liquid argon gets heated up faster as pointed out earlier. The maximum error for the shorter nanochannel occurs at 1350 ps (analy tical temperature = 109.14 K) and the value is 5.01 K. For the larger nanochannel, the error rema ins nearly the same ~ 3.4 K during the final 1000 ps of simulation. This justif ies that the error between MD simulation and heat equations analytical result depends on length scale of the problem. 78 PAGE 79 Figure 213. Temperature gradient along zdirection at different timesteps for Ti=100 K, Tup=110 K and Tlow=90 K (smooth lines repres ent analytical solutions) 79 PAGE 80 Figure 214. Temperature gradient along zdir ection at different timesteps for h=23.657 nm, Ti=100 K, Tupp=120 K and Tlow=120 K (smooth lines represent analytical solutions) 80 PAGE 81 Figure 215. Temperature gradient along zdir ection at different timesteps for h=44.357 nm, Ti=100 K, Tup=120 K and Tlow=120 K (smooth lines repres ent analytical solutions) 81 PAGE 82 Figure 216. Average domain temperatures from MD simulation and anal ytical equation with heating from both platinum walls for different nanochannel heights Maxwell speed distribution The Maxwell speed distribution (MSD) is a prob ability distribution describing the spread of molecules which are moving around at a given spee d. It is derived under the assumption of an ideal gas. Using elementary statistical mechanic s, it is easy to find out that the MSD must be proportional to the probability that a particle is moving at a given speed. The MSD can be derived and shown to be equa l to the following function: 23/2 22()4 2mvkTm Dv ve kT (2.21) where m is the mass of the molecule, k is Boltzmanns constant ( ) and T is the temperature. As this formula is a normalised probability distribution, it gi ves the probability of a molecule having a speed between v and v+dv. 2311.3810JK 82 PAGE 83 The MSD will be used to provide another verification of the proposed fluidwall thermal equilibrium model. At high temperatures, near the critical temperat ure, argon vapor would behave like an ideal gas. Since the proposed h eat transfer model is responsible for imparting kinetic energy to the argon atoms and evaporation, the MSD will serve as a good check as to if the molecular speeds of the vaporized atoms follow the expected normalized probability distribution. In order to do s o, a computational domain is set up in the form of a cuboid expressed in Cartesian c oordinates with the zdir ection denoting the height of the domain. The x, y and z dimension are 11.061 nm, 11.061 nm a nd 29.954 nm respectively. The four boundaries in x and y directions are periodic. A Pt wall constitutes the lower boundary in the zdirection. A mirror boundary condition is set up at the upper boundary in the zdire ction; this implies that the upper boundary is elastic and adiaba tic in nature. A thin liquid ar gon film of thickness 2 nm is placed on the Pt wall. The remaining volume is occupied by argon vapor. The simulation domain contains 5638 argon atoms and 5560 Pt atoms. The time step is 5 fs. Both liquid and vapor are initially at their respective saturation states at 100 K, and equilibrium is attained in the first 1000 ps. After this time period, the liquid atoms at the wall are heated up us ing the proposed fluidwall thermal equilibrium model with the wall te mperature set to 150 K. The simulation is run till 3000 ps. Figure 217 shows the va riation of number of liquid at oms with time. As it can be seen, evaporation occurs and a nonevaporating film forms on the lower wall similar to as obtained in the previous simulations. The numbe r of liquid atoms form ing the nonevaporating film is ~ 513, and the final average temperatur e is ~ 135 K. Increasing the domain height (i.e. increasing the vapor volume) will results in comp lete evaporation. However, since the aim here was to compare the velocity distribution of the argon vapor atoms in the presence of the nonevaporating film, the domain height was take n accordingly. Figure 218 shows the velocity 83 PAGE 84 distribution of the vapor atoms at t=3000 ps along with the MSD function at the final system temperature. In Fig. 218, N is the number of argon vapor atoms whose velocity values lie between V and V+dV, N is the to tal number of vapor atoms and V is the increment of velocity values. It can be seen that the simulation result s match well with the MSD function; this verifies that the velocities imparted by the proposed fl uidwall heat transfer model follow expected physics and at high temperatures argon vapor behaves like an ideal gas. Figure 217. Variation of number of liquid atoms with time in the simulation for verification with Maxwell Speed Dist ribution function 84 PAGE 85 Figure 218. Velocity distribu tion of vapor argon atoms at final equilibrium with a nonevaporating film at the lower wall Adiabatic Walls Adiabatic boundaries are very common in heat transfer problems. Thus, the necessity of mimicking adiabatic walls in molecular dynamics simulations is important. As per definition, an adiabatic process is one where no heat transfer occurs across the system boundaries. We simulate a nanochannel formed of platinum walls with liquid argon in between. The platinum walls are present in the zdirection, while the four boundaries in x and y di rections are peri odic in nature. The dimensions of the argon space in x, y and z directions are 5.78 nm, 5.78 nm and 8.585 nm. The initial equilibrium temperature is 100 K, and the density of liquid argon corresponds to the liquid saturation density at 100 K. A total of 5632 argon atoms and 3066 platinum atoms are present. The wallfluid interac tion is governed only by the Le nnardJones potential, and the fluidwall thermal equilibrium model is not applied in addition. This forms the basis of adiabatic wall in molecular dynamics simulations. After the velocityscalin g step for the initial 85 PAGE 86 500 ps, the argon atoms are allowed to interact fr eely for an additional 900 ps. Figures 219 and 220 show the variation of average temperatur e and energy with time for two different wall wettability values. The wettabil ity of a surface depends on the value in the LennardJones equations. As it can be seen, the temperatur e and energy remain constant throughout the simulation thus showing that, by simulating the fluidwall interaction only by the LennardJones potential, an adiabatic wall can be simulated. Figure 219. Variation of average temperature and energy with time for adiabatic wall with higher wettability 210.89410 J 86 PAGE 87 Figure 220. Variation of average temperature an d energy with time for adiabatic wall with lower wettability 210.35810 J 87 PAGE 88 CHAPTER 3 TRANSIENT FILM EVAPORATION Transient Film Evaporation Evaporation of a thin liquid argon film in ar gon vapor surrounding is simulated to study the phase change process at nanoscale. As theori zed in the literature, nonevaporating thin films on the order of nanometers occur in the micro region during bubble nucleation, as well as during meniscus evaporation. This part of the work aims to study nanoscale film evaporation and observe the physics as the film evaporates. Various parameters like temperature, density, etc are evaluated during the process which is not possible to do in experime nts due to the length scale of the problem. The simulation domain is a cuboid of dimensions 8.011 nm 8.011 nm18.961 nm. An ultrathin liquid argon film of 2 nm thickness is placed on the Pt wall. The boundaries in x and y directions are periodic in nature. The upper boundary in the zdirection is an imaginary adiabatic wall (mirror boundary condition), and the Pt wall forms the lower boundary. The remaining volume is occupied by argon vapor. Both liquid a nd vapor are initially at their respective saturation states at 90 K. The simulation domai n initially comprises of 2645 liquid film argon atoms, 112 argon vapor atoms and 2929 Pt atoms. Th e initial structure of the domain is shown in figure 31 along with the domain at an intermed iate timestep. The velocityscaling, equilibration and the evaporation (heating) periods are from 01000 ps, 10002000 ps and 20006000 ps respectively. The time step is 5 fs. The Pt wall temperature is stepjumped from 90 K to 130 K at the starting of the evaporation period. 88 PAGE 89 Figure 31. Snapshot of transi ent film evaporation domain Figure 32 shows the variation of average sy stem temperature and energy with time. The system temperature fluctuates about a mean va lue of 90.43 K and the system mean energy is 0.03145 eV/atom, during the equilibra tion period. When the Pt wall temperature is increased, the liquid film starts evaporating and the system temperature and energy increase. The average system temperature evaluated in the final 1000 ps is 120.62 K. The curve fit for energy shows that the final equilibrium energy is 0.001545 eV/atom. The stepchange in wall temperature incorporat ed in this work theoretically corresponds to an infinitely thermal conductive wall, and th e heating time period of the wall is neglected. This was again taken as a reasonable assumpti on for highly conductive metallic materials since 89 PAGE 90 the heating time period for threemo lecular layers of platinum would be in the order of 110 pico seconds or less. The heating period for the thin argon film is roughly 600 pico seconds (as seen in Figure 32). Due to this large difference, it is reasonable to neglect the wall heating time and not affect the physical conclusions of the cu rrent work. Maruyama and Kimura (1999) used Langevin thermostat to increase the temperature of Pt wall atoms by 10 K; they conclude that the solid wall temperatures quick ly respond to the temperatur e change, which validates our assumption of neglecting the heating time period for platinum. Figure 32. Variation of average system temperature and energy with time for transient film evaporation The number of liquid film argon atoms and va por argon atoms were determined at each time step and the plot is shown in figure 33. During the equilibration period, the number of liquid argon atoms varies to some extent due to the continuous transfer of atoms from the liquid phase to the vapor phase and vice versa at the in terface. The lower film starts to evaporate at 2000 ps and the evaporation rate exponentially decreases with time until no further evaporation takes place. The film is said to be nonevapor ating on the Pt surface. The curve fit shows this 90 PAGE 91 nonevaporating liquid film to be made up of an average of 868 atoms. This behavior will be explained later on. Figure 33. Variation of liquid film argon atoms in the domain with time for transient film evaporation The snapshots of YZ plane of the comput ational domain are shown at different time steps in figure 34. After the end of equilibrati on period, the liquid film can be distinguishably seen clearly defining an interface (Fig. 34a). Surface tension assist s in stabilizing the interface. The Pt wall atoms are denoted by the triangular symbols in red, the argon atoms forming the liquid film are shown in dark grey and the argon atoms in the vapor phase in light grey. Figures 34b through 34e demonstrate vaporization of the liquid film. The interface is distorted during the process, and the film become unstable as it evaporates. Figure 34f through 34h show the nonevaporating liquid film. A gap between the liquid film and the Pt wall can also be noticed which is due to the repulsive na ture of LennardJones potential. The number density profile in the zdirection of the domain is shown in figure 35 at different time steps. It is dete rmined by dividing the height of the domain into 55 equal slices, 91 PAGE 92 and counting the number of argon atoms in each slic e. The plot also shows the position and the thickness of the interface. The liquid film thickness can be seen to be decreasing initially and the system attains equilibrium towards the end of the simulation as there is negligible change in the number density profiles. Figure 36 shows the averaged temperature di stribution along the height of the domain for different time intervals. For the time in terval t=20003000 ps, during which maximum net evaporation takes place, a positive temperature gradient is established along the height of the domain. When the evaporation process begins, atoms from the liq uid film move to the vapor region; as the height of the domain is comparable to the mean free path of the atoms, the region near the upper imaginary adiabatic wall gets heat ed up more compared to the rest of the domain. This also explains the initial temperature jump which can be seen in figure 32. For the later times when no net evaporation occurs, the vapor atoms mix among themselves to obtain a nearlyconstant temperature distribution. Thus it can be concluded that the dominating mode of phase change heat transfer in nanochannels is molecular evaporation accompanied with high evaporation rates. Figure 36 also shows a large temperature grad ient in the nonevapora ting film. There is a temperature jump at the solidliquid interface, similar to Mar uyama and Kimura (1999) which was regarded as the thermal resistance over the interface. The temperature jump obtained in the present study is 2.8 K, and the temperature gradient in the nonevapora ting film is 9.64 K/nm. 92 PAGE 93 Figure 34. Projection of the YZ plane of the s imulation domain at different time intervals for transient film evaporation 93 PAGE 94 Figure 35. Number density profil es at different time steps for transient film evaporation Figure 36. Averaged temperature distributi on profile for transient film evaporation 94 PAGE 95 Figure 37 shows the averaged number density profiles along the height of the domain for different time intervals. The av erage interface thickness is defined as the distance over which the largest gradient in density occurs. The interface po sition is taken to be midway over the interface thickness. The interface position an d the thickness during the equilibration period are 2.241 nm and 0.689 nm respectively, whereas their valu es during t=50006000 ps are 1.362 nm and 1.0343 nm respectively. Figure 37. Averaged number density pr ofile for transient film evaporation The reason for the nonevaporating film can be attributed to the following factors: pressure in the vapor region, temperature of th e wall and the wallfluid interaction force. The pressure increase in the vapor region causes hi ndrance to evaporation, and thus increases the condensation rate. The Van der Wa als force of attraction (govern ed by the wallfluid interaction model) also plays a role in determining the th ickness of the noneva porating film making the value of energy scale parameter in the LJ potential crucial. The temperature of the wall determines the kinetic energy provided to the liqui d atoms. The nonevaporating film acts as an insulation due to which no further heat transfer occurs and the heat transfer coefficient value 95 PAGE 96 become zero, similar to the results obtained before the interline regi on for a nonevaporating liquid film. Heat Flux and Evaporation Rates The heat transfer rate and the net evaporati on rate per unit area are shown in figure 38. They are obtained by differentiating the curve fits of energy (Figure 32) and number of vapor atoms (Figure 33). Figure 38 also shows the time when 90% and 95% of total evaporation and heat transfer is completed. The maximum heat flux and evaporation rates are 603.39 MW/m2 and 4208.5 kg/s.m2 respectively. It should be noted that th e net evaporation rate is the difference between the evaporation and conde nsation rates. Thus when the net evaporation rate becomes zero it means that the evaporation rate equals the condensation rate. As this is an ultra thin film, the majority of atoms in the film keep eva porating and are replaced by the condensing vapor atoms. The film can be considered to be adso rbed at a given instant in time between ~ 40006000 ps. Evaluation of Hamaker Constant fr om Molecular Dynamics Simulation A Hamakertype analysis is done here to obtain the value of A for PtAr interaction which is not available in the literature. To the best of our knowledge, this kind of analysis involving the calculation of Hamaker constant using molecula r dynamics has not been attempted before. In addition, we would like to stress the point that th is method outlined below is applicable to other solidwall combination as long as their interaction potential is appropriate and depicts real situations. The LennardJones potential has been widely used and accepted for ArPt interaction in literature. The interaction between the Pt and Ar atoms is represented by the Lennard Jones potential: 126()4LJUr rr (3.1) 96 PAGE 97 where ArPt and ArPt The interactive energy per unit area between an argon plate of thickness t at a distance d from a semiinfinite Pt wall is determined analogous to the way shown by Stokes (2001). Figure 38. Evaporation flux and heat flux rates for tr ansient film evaporation 97 PAGE 98 The interaction between a Pt atom and an Ar bloc k is first considered as shown in figure 39a. The volume of the ring shaped element equals 2 ydyd Figure 39. Sketch showing the coordinates for ca lculating the interaction energy between (a) a Pt atom and Ar block (b) semii nfinite Pt block and Ar block The interaction potential between the single Pt atom and all the atoms in the small volume element becomes: 126()4 2av ArN dUr ydyd rrM (3.2) where av A rNM is the number of atoms per unit volume in the Ar block. Representing r in terms of z and and rewriting equation 3.2. 12 6 22 22()4 2av ArN dUz ydyd M zyzy (3.3) Integrating equation 3.3 for all values of y from zero to infinity and of from zero to t: 98 PAGE 99 12 6 22 22 00()4 2t av ArN dUz ydyd M zyzy (3.4) Carrying out this integration gives: 66 339151522 ()4 90()()av ArN Uz6 9 M zztzzt (3.5) Now to determine the interaction energy between th e semiinfinite Pt block and the Ar block, the atom shown in Figure 39b is considered to be one of many atoms in the Pt block. Thus all the atoms in the Pt block located at a distance x from the Ar block have an energy given by equation 3.5. A volume element dz in the Pt block will have av P tNMdz atoms per unit area, and hence the interaction potential per unit area can be given by: 66 3399151522 ()4 90()()av av Ar PtNN dUz dz MMzztzzt 6 (3.6) Integrating equation 12 for z = d to gives the interaction energy per unit surface area of a semiinfinite Pt block and an Ar block of thickness t as: 66 22811 () 12()3030() A Ud dtddtd 8 (3.7) The last two terms in equation 3.7 appear due to the repulsive terms of the LJ potential. Here A is a Hamakerconstant parameter which is derived to be: 624av av Ar PtNN A MM (3.8) This expression is similar to the general expressi on in literature as seen from equation 1.1. The modified LJ potential (equation 3.1 ), used in this simulation, is used to determine the energy of interaction per unit surface area between the liquid argon atoms and Pt wall atoms for the time 99 PAGE 100 period 50006000 ps at time intervals of 50 ps. Th e nonevaporating liquid f ilm is treated as the plate of thickness t and the thickness of the gap being d is the distance between the wall and the closest liquid film atom to it. The average value of A is obtained, from equation 13 as with a standard deviation of 7.68%. The typical value of Hamaker constants are of the order of 1020 and 1021 J (Butt et al., 2005). For further validation, the value of A is also obtained from equation 3.8 which does not us e force evaluation. This value equals which shows good agreement with the one obtained above. This analysis attempts to outline a general procedure for evaluating the Hamaker constant by molecular dynamics simulation, thus providing a link for us ing molecular dynamics simulation results in continuum modeling. 205.275910 J208.556310 J Parameters Governing the Format ion of NonEvaporating Film The previous section is extended further to gain better understanding on the occurrence of the nonevaporating film, and the parameters which govern its formation. The mirror boundary condition is replaced by an addi tional platinum wall with a thin film, thus making the problem symmetric in nature. Vapor pressure is also evaluated in the simulations. Seven cases are simulated to study the effect of nanochannel he ight, initial film thickness and temperature on heat flux rate, evaporation fl ux rate and Hamaker constant. The simulation domain is a nanochannel in th e form of a cuboid as shown in figure 310. The x and y dimensions are 6.117 nm 6.117 nm respectively; zdirect ion represents the channel height. The four boundaries in x and y directions are periodic; two Pt walls constitute the boundaries in the zdirection. A th in liquid argon film of thickness tfilm is placed on each Pt wall. The remaining volume is occupied by argon vapor. Both liquid and vapor are initially at their respective saturation states at 90 K. The simula tion domain contains 3465 Pt atoms. The time 100 PAGE 101 step is 5 fs. The velocityscali ng period is 0500 ps, the equilibration period is from 5001000 ps and the heating period is from 10004000 ps. Nanochannel height h (case 13) and initial film thickness tfilm (case 46) are varied, and simultaneous evaporationcondensation (case 7) is simulated. The details of seven cases studied are listed in Table 31. Figure 310. Snapshot of simu lation domain of nanochannel Table 31. Parameters for seven cases simulated in nanochannel h (nm) tfilm (nm) Tup wall (K) Tlow wall (K) No. of argon atoms Case 1 16.32 3 130 130 4669 Case 2 25.5 3 130 130 4714 Case 3 35.7 3 130 130 4768 Case 4 25.5 2 130 130 2989 Case 5 25.5 4 130 130 6439 Case 6 25.5 5 130 130 8164 Case 7 25.5 3 130 90 4714 Figure 311 shows the variation of average sy stem temperature and energy with time for case 4. The system temperature fluctuates abou t a mean value of 90.99 K and the system mean 101 PAGE 102 energy is 0.0342397 eV/atom, duri ng the equilibration period. When the Pt wall temperature is increased, the liquid film starts evaporating and the system temperature and energy increase. The average system temperature evaluated in the fi nal 1000 ps is 123.99 K. The curve fit for energy shows that the final equilibrium energy is 0.0088331 eV/atom. Figure 311. Variation of average system temperature and energy with time for case 4 The number of liquid film argon atoms and va por argon atoms were determined at each time step and the plot for vapor atoms is shown in figure 312. The films start to evaporate at 1000 ps and the evaporation rate exponentially decreases with time until no further evaporation takes place. The snapshots of xz plane of the co mputational domain are shown at different time steps in figure 313, similar to figure 34. Afte r the end of equilibration period, the liquid film can be distinguishably seen clearly defining an interface (Fig. 313a). Surface tension assists in stabilizing the interface. The Pt wall atoms are denoted by the triangular symbols in red, the 102 PAGE 103 argon atoms forming the liquid films are shown in dark grey and the argon atoms in the vapor phase in light grey. Figures 313d and 313e show the formation of nonevaporating liquid film. Figure 314 shows the averaged number densit y profiles along the height of the domain for different time intervals. It is determined by dividing the height of th e domain into 90 equal slices, counting the number of argon atoms in each slice and averaging it over time. The liquid film thickness can be seen to be decreasing initially. The system attains equilibrium towards the end of the simulation as there is negligible ch ange in the number density profiles confirming the formation of a nonevaporating film. Figure 312. Variation of vapor argon atoms with time for case 4 103 PAGE 104 Figure 313. Snapshots of domain for case 4 Figure 314. Average density of domain at different times for case 4 104 PAGE 105 Heat Flux, Evaporation Flux, Ha maker Constant for Nanochannel The heat transfer rate and the evapor ation rate per unit area are obtained by differentiating the curve fits of energy (Figur e 311) and number of vapor atoms (Figure 312) for each case. The net heat flux and evaporati on flux are evaluated by ca lculating the area under the respective curve fits. Pressure in vapor region is co mputed at time intervals of 50 ps using the method verified before. The initial and final vapor pressure is calculated by averaging values from 7001000 ps and 30004000 ps respectively. The liquid film thic kness is determined from the density profile at the position where the liquid de nsity begins to decrease. The lo cation of the equi molar surface, obtained from the derivative of de nsity profile, coincides with the extremum of the gradient of the density profile (Mareschal et al., 1997), as is shown in figure 315 for case 6. Hamaker constant is evaluated as described in section 3.1.2. Figure 315. Densityheight plot showing noneva porating liquid film thickness and equimolar surface position for case 6 105 PAGE 106 Varying Nanochannel Height h Figures 316 and 317 show heat flux rate vs. time and evaporation rate vs. time respectively for varying nanochannel height with constant initial f ilm thickness (cases 13). From the figures, it can be seen that nanochannels ha ve the capacity to produce very high heat transfer rates compared to macroscale heat transfer. The highest starting heat transfer rate occurs in the shortest nanochannel (case 1) wh ich is due to the occurrence of high initial evaporation rate since the channel height is small, the high en ergy vapor atoms which evaporate from, say, the lower liquid film have greater chances of reach ing and striking the upper liquid film in turn causing additional evaporation. Heat flux rate and evaporation rate decrease exponentially with time. The decrease is more prominent in shor ter height nanochannels as the vapor space available for evaporation is less. Thus heat tr ansfer and evaporation continues on for a longer time with increasing nanochannel height. Table 32 shows final equilibrium properties, net heat flux and evaporation flux for cases 13. The thickness of obtained nonevaporating f ilm decreases with increase in nanochannel height. The interface becomes less spread out with increasing nanochannel height. The final vapor pressures for the three cases are nearly the same, which justifies the calculated Hamaker constants to be similar as it depends on pressure temperature and solidfluid interaction. Net heat and evaporation fluxes (area under the curves) increase with increasing nanochannel height, which is as expected. 106 PAGE 107 Figure 316. Heat flux rate with time for varying nanochannel he ight (cases 13) Figure 317. Evaporation rate with time for varying nanochannel height (cases 13) 107 PAGE 108 Table 32. Case 13 results for film evaporation in nanochannel Nonevaporating film thickness (nm) Equimolar liquidvapor position (nm) Final vapor pressure (MPa) Net heat flux (J/m 2 ) Net evaporation flux (kg/m 2 ) x 10 6 Hamaker constant A (J) x10 20 Case 1 1.93 3.412 1.7065 0.1811 0.640086 6.756 Case 2 1.61 2.975 1.7031 0.2154 1.00314 6.614 Case 3 1.57 2.713 1.7574 0.2654 1.46444 6.623 Varying Film Thickness tfilm Figures 318 and 319, respectively, show the heat flux rate vs. time and evaporation rate vs. time for varying initial film thickness with a c onstant nanochannel height (cases 2, 46). The initial evaporation rate, from figure 319, is maximum for tfilm = 2 nm and minimum for tfilm = 5 nm. This shows that initial evaporation rate is in versely proportional to th e thickness of the liquid film as heat has to be transfe rred across larger number of atom s. Also, compared to figure 317, it can be said that thickness of the liquid film plays a more dominant role to higher initial evaporation rate than the height of nanochannel. The inset of figure 318 shows a somewhat odd tr end for the initial heat transfer rates. For cases 13, discussed above, ini tial evaporation rate was directly proportional to initial heat transfer rate, but it should also be noted that the number of liqui d film atoms were almost same. In the present discussion, the numbe r of liquid atoms vary to a grea t extent. Thus the thickness of the liquid film also determines how much heat ca n be initially transferred from the wall to the liquid as more liquid atoms store mo re energy. This explains why the initial heat transfer rate for tfilm=5 nm (case 6) > tfilm=4 nm (case 5) > tfilm=3 nm (case 2). For tfilm=2 nm (case 4), the evaporation rate (i.e. heat used up in phase change) is dominant making it have the highest initial heat transfer rate. The evaporation rate and heat flux rate exponentially drop with time, being more prominent for thinner films as they have lesser number of atoms. 108 PAGE 109 Table 33 shows final equilibrium properties, net heat flux and evaporation flux for cases 2, 46. The thickness of obtained no nevaporating film increases with increasing initial film thickness, and the interface becomes more spre ad out. The final vapor pressures increase accordingly for cases 2, 4 and 5, which justifies th e calculated Hamaker constants to increase as well. It is interesting to note here that for a constant nanochannel height, increasing the thickness of liquid film beyond a certain point does not signif icantly increase the final vapor pressure as can be seen for cases 5 and 6. Net heat flux incr eases with increasing film thickness as there are more atoms in the nanochannel to evaporate and/or store energy. Net evaporation flux first increases with tfilm and then decreases, due to the eff ect of two competing parameters: total number of liquid atoms and vapor volume in the nanochannel. Thus, to obtain net maximum evaporation flux there is an optimum initial liq uid film thickness for any specific height of nanochannel. Table 33. Case 2, 46 results for film evaporation in nanochannel Nonevaporating film thickness (nm) Equimolar liquidvapor position (nm) Final vapor pressure (MPa) Total heat flux (J/m 2 ) Total evaporation flux (kg/m 2 ) x 10 6 Hamaker constant A (J) x10 20 Case 4 0.99 2.125 1.6122 0.163 0.966156 5.229 Case 2 1.61 3.258 1.6798 0.215 1.00314 6.613 Case 5 2.97 4.392 1.7490 0.266 0.956312 7.023 Case 6 4.67 5.525 1.7525 0.321 0.846324 7.026 109 PAGE 110 Figure 318. Heat flux rate with time for varying film thickness (cases 2, 46) Figure 319. Evaporation rate with time for varying film thickness (cases 2, 46) 110 PAGE 111 Vapor PressureDensity Plots The average initial and final vapor pressures and vapor densities we re evaluated for the simulation runs for all six cases (16). Vapor pres sure vs. vapor density for initial and final state is plotted for varying nanochannel height ( cases 13) in figure 320, and for varying film thickness tflim (cases 2, 46) in figure 321. The pl ots show that the states lie on the saturation curve, which is obtained from thermodynamics table of argon, implying that the vapor is saturated. This further validates the existence of the nonevaporating thin f ilm in all six cases (16) as the final vapor condition is saturated and thus no additional evaporation from the film can occur. The vapor cannot increase beyond saturation values in this work since the only mode of evaporation is from the thin film which beco mes nonevaporating after a certain time period. Since there is no other physical process occurring in the system to add more vapor atoms, the vapor pressure remains saturated. This define s the limiting condition for the formation of the nonevaporating thin film at the th reephase contact line for a pa rticular combination of vapor pressure, substrate temperature and solidliqui d interaction. It should be noted that the evaporating thin film region of the threephase contact line has a curved interface, unlike the nonevaporating film which has a flat in terface. The YoungLaplace equation states 2vlPPr ( Pvvapor pressure, Plliquid pressure, surface tension, r radius of curvature). For the evaporating thin film region which has a curvature Pv > Pl. If the pressure and temperature in the liquid and vapor are (Pl, Tl) and (Pv, Tl) respectively (liquid and vapor are in thermal equilibrium) and the vapor phase is in a saturated state, the vapo r pressure on a curved interface at equilibrium is given by5 ()exp2vsatlllPPTrRT, where Psat is the saturation pressure corresponding to Tl on a flat surface. Thus, Pv < Psat, which implies that curvature of the interface lowers the vapor pressure compared with that above a planar interface for the same 111 PAGE 112 liquid temperature25. It follows from above that Pl PAGE 113 Figure 321. Pressuredensity plot of initial and final states fo r varying film thickness (cases 2, 46) Case 7 Simultaneous Evaporation and Condensation It has been widely speculated in literature that wallfluid molecular interaction is the dominant factor for the formation of a nonevapor ating film, the other parameters being pressure and temperature. In each of the cases (16) studi ed so far, a nonevaporating film was obtained. Wallfluid interaction and wall temperatures were the same, with evapor ation occurring from both walls. The dependence of Hamaker constant of the nonevaporating f ilm on vapor pressure was shown. It is easy to visualize the effect of nominal increase in temperature of the wall which will lead to more energy imparted to atoms, more evaporation and hence thinner nonevaporating films, if at all they have to form. Thus, case 7 studies as to how impor tant wallfluid molecular interaction is to the formation of nonevaporating film. Similar to case 2, two thin films each of thickness 3 nm is placed on upper and lower wa lls in a nanochannel of dimensions 6.117 nm6.117 nm25.5 nm. In this case, only the upper wa ll temperature is increased to 130 K 113 PAGE 114 while the lower wall temperature is kept the sa me as the initial domain temperature of 90 K. Figure 322 shows the variation of vapor pressure with time. The vapor pressure increases when evaporation is dominant and then decreases as more condensation occurs. Figure 323 shows the snapshots of the domain during the evaporationc ondensation process. Nonevaporating film is not obtained on the upper wall even though the wall is at the same temperature of cases 16. Wallliquid molecular interaction is important to th e extent that it determines the strength with which the solid molecules attract liquid molecule s i.e. governs wettability. But this case shows that pressure is the dominant factor for the form ation of nonevaporating f ilms. In other words, for a certain wallliquid pair, nonevaporating film will form or not form depending on the vapor pressure of the domain and the temperature of the wall. Figure 322. Vapor pressure variation with time for case 7 114 PAGE 115 Figure 323. Snapshots of domain at different time intervals for case 7 115 PAGE 116 CHAPTER 4 NANODROPLET IMPACT ON A HOMOGENEOUS SURFACE Impact of argon nanodroplet on homogene ous platinum wall is simulated using molecular dynamics. Wall temperature (90 K, 300 K) and surface wettability (high, low based on ArPt parameter) are varied to study four cases. An external forc e is applied on the drop for a short time period to achieve velocities in the range of 13 m/s 15 m/s. The simulation domain is a cuboid of dimensions 15.721 nm 15.721 nm13.936 nm. A liquid droplet, initially in the form of a cube of side 6.2 nm, is placed in the domain with the drop center at z=8.455 nm. The remaining volume is o ccupied by argon vapor. Bo th liquid and vapor are initially at their respective saturation st ates at 90 K. The simulation domain initially comprises of 4913 liquid droplet argon atoms, 315 argon vapor atoms and 11229 Pt atoms. The time step is 5 fs. The argon atoms are subjected to six boundaries, two in each of the x, y and z directions. The boundaries in the x and y directions are pe riodic throughout the simulation. The simulation starts off with a velocityscaling period (0500 ps) followed by the equilibration period (5001000 ps). The forceapplication period follows the equilibration period where a constant external force of magnitude N is applied to each liquid drop atom in the negative z direction during the time period of 10001301.275 ps in or der to attain a final drop velocity in the range of 13 m/s to 15 m/s. The boundaries in the zdi rection were also periodi c in nature so far. 153.508310At the end of the forceapplication peri od, the periodic boundary conditions in the zdirection are removed and the Pt wall is inserted into the simulation at the desired temperature. Now, the mirror boundary condition is applied at the upper boundary, wh ile the lower boundary in the zdirection becomes the Pt wall. The simulation is continued to run till 3000 ps. 116 PAGE 117 Maruyama and Kimura (2000) concluded from their work that the wettability of the surface can be change by varying the ArPt parameter in the LJ potenatial. Making use of this fact, droplet impact is studied for high and low wettabl e surfaces with variation in wall temperature. Four cases which are looke d upon in this work are: Case 1: Wall temperature at 90 K and 21 0.89410ArPtJ Case 2: Wall temperature at 90 K and 210.357610ArPt J Case 3: Wall temperature at 300 K and 21 0.89410ArPtJ Case 4: Wall temperature at 300 K and 210.357610ArPt J Evaluation of Drop Temperature At each time step, the total number of liquid atoms, N, constituting the drop can be found out. The drop temperature is determined using th e following equations, also used by Hanasaki et al. (2004): 1 N 2 x xi iTv (4.1) 1 N y iTv2 yi (4.2) 2 *2 ,, 1 N zz iz iTvv i (4.3) ***1 3xyzTTT N *T (4.4) Evaluation of Interface Markers, Interfacial Fit and Contact Angle The contact angle which the drop makes with th e Pt wall is determined in 2D in the xz plane. The drop is divided into squares of dimension 1ArAr 1ArAr and the number of atoms in each square are determined. Once all the at oms are allotted to their squares, the average 117 PAGE 118 number density is found out. The square with th e maximum number of atoms is determined and a check is performed from this square in the ou tward direction such that if the number of atoms in a square falls below 0.5 times the average num ber density, an interface ma rker is placed at the center of that square. Thus, in terface markers are placed around the drop using this procedure. After all interface markers are determined, an elliptical fit of nature 22()()ccxxzz1 is used to fit all the markers such that where xc and zc are the x,z coord inates of center of mass of the drop. The contact a ngle is the angle made by the tange nt to the elliptical fit at the point where it intersects with the wall i.e. z=0. Th e contact angle is evaluated in time intervals of 5 ps. 20.999 R Droplet Impingement Case 1 Figure 41 shows the variation of average dr op temperature and ener gy with time. The drop temperature fluctuates about a mean valu e of 89.13 K and the drop mean energy is 0.03035 eV/atom, during the equilibration period. The average drop temperature and energy evaluated in the final 1000 ps is 89.95 K and 0.02498 eV/ato m respectively. Figure 42 shows the drop zvelocity curve and the center of drop position with time. The drop attains a velocity of 13.69 m/s at the end of the forceappl ication period and travels at that velocity for a short time period. When the drop reaches within a distance rc from the wall, the long rang e attractive force from the wall atoms act of the drop atoms and thus increas e the velocity of the drop. The average velocity of the drop on impact is 26.63 m/s. The increase in the velocity depends on the size of the drop since in largersized drops; only the atoms at th e lower end of the drop are influenced by the wall atoms and the majority of the atoms forming the drop do not interact with the wall during impact. Thus the wettability of the wall is found to be significant in determining the actual impact velocity of nanosized drops. 118 PAGE 119 Figure 41. Variation of drop te mperature and energy with time for drop impact (case 1) Figure 42. Variation of dropcente r height and drop velocity with time for drop impact (case 1) Figure 43 shows the snapshots of the xz plane of the simu lation at different timesteps. The interface markers, interfacial fit and the contact angle values are shown in the snapshots. It can be seen that due to the high wettability of the wall, the drop spreads completely on impact. 119 PAGE 120 Please note that the wall is not actually present in Figure 43a (wall has not been inserted yet) and is just shown for convenience. Figure 43. Snapshots of xz pl ane for drop impact (case 1) Figure 44 shows the variation of contact a ngle with time during impact and spreading. The drop impacts on the wall at about 1445 ps and the contact angle decreases as the drop spreads. As the spreading edge of the drop reaches x and y bounda ries, the edges interact with each other due to the periodic nature of th e boundaries thus causing rapid spreading and formation of a thin film. This physically signif ies the drop merging with four other identical drops. The contact angle becomes zero once the fi lm forms. A larger computational domain 120 PAGE 121 needs to be used to notice comple te individual droplet sp reading (it was not used in this current work due to computational restrictions). Figure 44. Variation of contact angle with time for drop impact (case 1) Droplet Impingement Case 2 Figure 45 shows the snapshots of xz plane at different tim esteps. Figure 46 shows the variation of drop velocity and drop center height with time for the less wettable surface. The drop achieves a zvelocity of .64 m/s at the end of forceapplication period. Similar to case 1, the average drop velocity increases lower end of the drop comes in the long attractive force of wall atoms and the actual velocity during impact is 19.49 m/s, which is less than that obtained in case 1 as the wettability is lower in this case. The drop impacts on the wall at about 1465 ps. It can be seen that the drop does not spread on this lesswettable su rface. The drop center height in figure 46 is seen to oscillate thus showing that the drop os cillates on a nonspreading surface after impact, and the contact angle varies. The plot of contact angle with time (Figure 47) confirms the fact. The mean value of the contact angle during the final 500 ps is 126.7o. The curve fit for drop oscillation is shown from 25003000 ps. 121 PAGE 122 Figure 45. Snapshots of xz pl ane for drop impact (case 2) Figure 46. Variation of drop cente r height and drop velocity with time for drop impact (case 2) 122 PAGE 123 Figure 47. Variation of contact angle with time for drop impact (case 2) Droplet Impingement Case 3 and Case 4 Figures 48 and 49 show the snapshots of the yz plane for case 3 and case 4 respectively. For the highwettability surface (Fig ure 48), leidenfrost effect is observed and a layer of vapor is formed as soon as the drop imp acts on the wall which is at 300 K. It can be seen that the drop does not actually touch the wall (unlike case 1 and case 2) but rests on the vapor layer throughout the eva poration process. Complete drop eva poration is obtained. For the lowwettability surface (Figure 49), the leidenfrost effect is again obser ved but to a greater extent as the vapor layer which forms during impact pushes the drop away from the wall. Due to the low attraction from the wall atoms, the drop does not hol d to the surface as was seen in case 3. The temperature of the drop increases only for the time period it is in contact w ith the wall. Thus, the evaporation rate significantly reduces for low wettable surfaces. The interface markers are not evaluated for these cases since evaluation of cont act angle becomes irrelevant. The light colored circles represent vapor atoms and the dark colored circles repr esent liquid atoms. 123 PAGE 124 Figure 48. Snapshots of yz pl ane for drop impact (case 3) Figure 49. Snapshots of yz pl ane for drop impact (case 4) 124 PAGE 125 CHAPTER 5 MENISCUS EVAPORATION AT NANOSCALE The moving contact line is a threephase phenomenon where vaporliquid interface meets the solid surface. It has been widely accepted in lit erature that the contact line at the base of a bubble in nucleate boiling can be divided into three regions: n onevaporating thinfilm region, evaporating thinfilm region and intrinsic menisc us region. Thus nucleate boiling is a multiscale problem. The nonevaporating region is of the or der of nanometers, and no mass/heat transfer from it. Maximum evaporation and heat transfer occurs from the evap orating thinfilm region and the liquid is fed from the bulk liquid through the intrinsic meniscus region. There is a region, called the interline regio n, between the nonevaporating thin f ilm region and the evaporation thin film region over which the film varies in thickne ss and curvature to accommodate the transition between the two regions. This is the thinnest portion of the meniscus over which vaporization can occur, and since it is the thinnest, it is also the location where the evaporation rate is the highest (Carey 1992). The thinfilm regions are also characteri zed by capillary and disjoining pressures. Researchers have studied the heat tr ansfer characteristics in the thinfilm region analytically and experimentally. However, measur ing properties at nanoscale is still a challenge and experiments conducted were on microto millim eter scale. Since the thinfilm region starts off at nanoscale, it would be more befitting to use molecular dynamics to study the nonevaporating and interline region in a nanoscale meniscus. Simulation Model A concave meniscus was formed by placi ng liquid argon between a lower wall and an upper wall, with an opening in th e upper wall as shown in Figs. 51 and 52. The walls are made of three layers of platinum (Pt) atoms arranged in fcc (111) structure. The space above the meniscus is occupied by argon vapor. The domai n consists of a total of 14172 argon atoms and 125 PAGE 126 7776 platinum atoms. The total width, hei ght and depth are 34 nm, 59.5 nm and 4.3 nm respectively. The length of th e opening in the upper wall is 28.22 nm. The separation between the lower wall and the upper wall is 6.2 nm. The initial equilibrium temp erature is 90 K. The time step is 5 fs. The atomic interaction is gove rned by the modified LennardJones potential. The cutoff radius is set as 4cutArArr All the boundaries in x and y directions are periodic. The length of the boundary above the upper wall in the xdirection is restricted to the width of the opening. Any argon atom which goes above the upper wall does not interact with the wall atoms anymore. The top boundary in the zdirection is the mirror boundary condition where the argon atom is reflected back in the domain without any loss of ener gy or momentum i.e. the boundary is adiabatic and elastic in nature. The fluidwall thermal equilibrium model is used to numerically simulate heat transfer between wall and fluid atoms. Equations and parameters implemented are nondimensional (Rapaport 1995). The algorithm used to calculate the atomic force in teractions is the linkedcell algorithm which is a cellbased method and invol ves data organization (Allen and Tildesley 1987). The integrator method used here is the VelocityVerlet method (Sadus 1999). Two cases are simulated in this work: Case I: After the equilibrium period, the temperature of platinum wall underneath the opening (shown in red color in Fig. 2) is changed to 130 K while the rest of the wall (shown in blue color in Fig. 2) is kept atsame initial temperature of 90 K. Case II: After the equilibrium period, the temper atures of all walls are changed to 130 K. 126 PAGE 127 6.20 nm 53.30 nm Periodic Boundary Condition (PBC) Periodic Boundary Condition PBC PBC Figure 51. Schematic showing configur ation of the computational domain 127 PAGE 128 Figure 52. 3D schematic of the computational domain Liquid atoms are distinguished from vapor atoms based on the minimum number of neighboring atoms within a certain radius, as di scussed ealier. For argon at 90 K, the threshold number density was obtained as 7 and threshold radius as 0.53 nm i.e. if an argon atom has 7 or more neighboring atoms within a radius of 0.53 nm it is said to be a liquid atom or else a vapor atom. As discussed earlier, the Hamaker constant for nonevaporating films is obtained from molecular dynamics simulations. Starting from the LennardJones potential, which is the model of interaction between ArPt, th e following equation is derived: 128 PAGE 129 66 228811 () 12()3030()ArPtArPt LJA Uz dzdz (5.1) where, A is a Hamakerconstant parameter, d is the gap between Ar and Pt slabs, z is the total thickness of Ar slab (incl uding the gap thickness) and ULJ is the total interaction energy between ArPt slabs from molecular dynamics using LJ poten tial. This equation was used to calculate the Hamaker constant for nonevapora ting film with varying pressure and temperature earlier, and an average value of J is used in this work. The disjoining pressure is calculated as: 206.1310 A 6 3() 8 2 12()30()LJ ArPt ddUz A P dzzz 9 (5.2) The capillary pressure is the product of interfacial curvature K and surface tension coefficient as follows: 1.5 2,1cPKK (5.3) where and are, respectively, the first and second derivatives of thickness with respect to length x The variation of meniscus thickness is determ ined in 2D in the xz plane. The meniscus, formed from liquid argon atoms, is divi ded into square bins of dimension 1ArAr 1ArAr and the number of atoms in each square are determ ined. Once all the atoms are allotted to their squares, the average number dens ity is found out. A check is performed from the Pt wall in the positive zdirection such that if the number of atoms in a bin falls below 0.5 times the average number density, an interface marker is placed at the center of that bin. Thus, interface markers are placed at the meniscus interface using this pr ocedure. A fourthorder polynomial fit of the interface markers is used to obtain the function (x) The evaporating meniscus, formed of liquid at oms, is divided into 10 equal regions as shown in Fig. 53 along the xdi rection. Regions 0 and 11 are denoted for nanochannel region. 129 PAGE 130 The length of each region is 2.822 nm. The number of liquid atoms, along with the average temperature and energy of each region is tracked over time. Figure 53. Schematic depicting division of nanochannels and liqui d meniscus into 11 regions Consider a liquid nanomeniscus on a surface surrounded by vapor, similar to this work, as shown in Fig. 54. The meniscus is divided into two regions, A and B, with equal surface area. At time t1, as shown in Fig. 54a, region A of the liquid meniscus has a total of N1a atoms with an average energy of E1a. Similarly, region B of the liquid meniscus has a total of N1b atoms with an average energy of E1b. The vapor region above the meniscus has N1v atoms with an average energy of E1v. Evaporation takes place from the me niscus and at some later time t2, the number of atoms and average energy in region A, region B and vapor changes to (N2a, E2a), (N2b, E2b) and (N2v, E2v) respectively, as depicted in Fig. 54b. The goal here is to find the average heat and evaporation flux rates from regions A and B. Please note that the number of atoms in the domain is constant. Also, as there is no flow in the dire ctions parallel to the su rface, the net amount of atoms transferred between regions A and B is assumed to be zero. 130 PAGE 131 Figure 54. Schematic showing two regions of an evaporating meniscus at times t1 and t2 Let the total number of atoms eva porated in the entire domain be Nev, and the atoms evaporated from regions A and B are Neva and Nevb respectively. Thus, 1 evaaaNNN2 2 (5.4) 1 evbbbNNN (5.5) 21 1122 evvvevaevbababNNNNNNNNN (5.6) Total energy of domain at time t1 111111 totalaabbvvENENEN1E Total energy of domain at time t2 222222 totalaabbvvENENENE2 Thus, net energy transferred to th e fluid from the heated surface: 21 nettotaltotalEEE (5.7) 221122112211 aaaabbbbvvvvNENENENENENE 222 1222 11 211aaaevaabbbevbbvevaevbvvvNENNENENNENNNENE Rearranging: 131 PAGE 132 221 21 2121 221 1neteva aaa va vvvv bbbevbb E NEENEENEENEENEE (5.8) The last term accounts for the part of energy tran sfer used to heat up the initial vapor atoms from E1v to E2v. This is done indirectly by th e higher energy atoms which are evaporated from the meniscus interface and heat up the vapor atoms already present in the vapor space. In order to divide this term for region s A and B, we must know how many initial vapor atoms have been heated up by evaporating atom s from region A and from region B, which is very complex to evaluate accurately. Thus, as an approximation, the last term is divided for regions A and B by introducing a weighting factor based on the fraction of atoms evaporated from each region as follows: 121 vvvNEE 121121121eva evb ev evvvvvvvvvvNN NNNEENEENEEE (5.9) Thus, the net heat transf er term can now be written as the su mmation of heat transfer from region A and heat transfer from region B: netABEE (5.10) where, 221 21121eva A eveva aaa vavvvN NENEENEENEE 21 2 221 1evb B evvv bbbevbbN NENEENEENEE 1 v v This can be generalized to obtai n the heat transfer from a meni scus divided into N regions. The heat and evaporation flux rates, between time t1 and t2, from a region R can be evaluated using the following equations: 132 PAGE 133 Heat flux rate = 21221 21 121evR evtotalevR RRR vR vvvN N AreattNEENEENEE (5.11) Evaporation flux rate = 12 21aaatomNNm A reatt (5.12) Results and Discussion Case I Figure 55 shows the snapshots of the computational domain at different time intervals. In Fig. 55a, the meniscus has a smooth curvatur e as evaporation has not yet started and surface tension assists in the formation of the interface. Vigorous evaporation is seen in Fig. 4b which results in an uneven meniscus interface. As evaporation rate slows down comparatively, a uniform but increased curvature is obtained in Fig. 55c. The decrease in evaporation rate is due to the increase in vapor pressure. With continuous evaporation taking place, the thinnest part of the meniscus at the center continues to decreas e in thickness (Figs. 55d, 55e) until a uniform nonevaporating film forms. The cu rvature of remaining part of the meniscus is sharper than before evaporation began. The formation of the nonevaporating film is again confirmed from Fig. 56. This figure shows the variation of liquid atoms with time fo r regions 5 to10. Regions 5 and 6, which are at the center of the meniscus, show no further ev aporation after 2000 ps t hus resulting in a nonevaporating film. Liquid atoms in regions 7 and 8 continue to evaporate, similar to the results obtained from analytical solutions of a macrosc opic evaporating meniscus. However, this work shows a clearer picture of the physics occurring at and in the vicinity of the nanoscale nonevaporating film which cannot be accurately si mulated using macroscopic theories. Regions 7 and 8 can be labeled as the interline region. Re gion 10, which is adjacent to the nanochannel, shows an increase in the number of liquid atoms. The reason is th at the average temperature of 133 PAGE 134 this region is lower than the vapor and thus c ondensation in this region causes the increase in liquid atoms (as also seen in Fig. 55). Figure 55. Snapshots of xz plane at different time steps for Case I 134 PAGE 135 Figure 56. Variation of li quid atoms with time for regions 510 for Case I The pressure in the vapor space is calculated at every time step, and is averaged over a time period of 50 ps (10,000 time steps). The pressu re variation with time is shown in Fig. 57. As the computational domain is closed, the pressu re increases from an initial value, averaged from 8001000 ps, of 0.154 MPa to a value of 0.630 MPa in the final 200 ps. The saturation temperature corresponding to this final pressure is Tsatpres = 109.17 K. Figure 58 shows the average temperature of each region from 8001000 ps and 28003000 ps. For regions 5 and 6, although the temperatures are hi gh, there is no evaporation afte r 2000 ps. A nonevaporating thin film can exist in equilibrium on a solid surface even when the vapor pressure in the surrounding gas is below the normal saturation pressure for the film temperature (Carey 1992). This reiterates that the film in regions 5 and 6 is nonevaporating. The net heat flux and evaporation flux rates are shown in Fig. 59 for regions 1 through 10. The heat flux values are averaged over t1 = 8001000 ps, t2 = 15001700 ps and t3 = 27002900 ps. The overall average heat flux rates are maximum for regions 4 and 7 since evaporation 135 PAGE 136 continues to take place. The init ial heat flux rates from regions 5 and 6 are higher until the nonevaporating film is formed, after which no furt her evaporation occurs. Regions 1 and 10 show negative heat and mass flux rates as the average temperature of these regions are lower than vapor and condensation occurs in th ese regions. As it can be seen from the values, high heat flux rate of the order of 100 MW/m2 and high evaporation rates of the order of 1000 kg/m2s are achievable from nanoscale evaporating meniscus Based on these values, very high heat/mass transfer rates exist during the early stages of bubble growth during the formation of the nonevaporating film at the bubble ba se; however as the time required for its formation is only on the order of nanoseconds, the high flux rates are onl y sustained for a very short period of time. Engineered manipulation of extending the time period of formation of nonevaporating film, or its breakage after formation can result in signifi cantly higher heat fl ux rates and delaying the critical heat flux limit at macroscale. Chen et al. (2004) performed experiment s, on a heating surface consisting of microheaters of individu al heater area of 270 m x 270 m to measure the power consumed throughout the cycle of bubble growth, detachment and departure. Our aim here is to substantiate the heat flux values obtained from molecular dyn amics simulation, and gain better understanding of the heat transfer characteristics. Although FC72 was chosen as the boiling fluid by Chen et al. (2004) and argon is chosen in this work, the Jakob number is one of the key parameters that can affect the contribution of heat transfer mechanisms in nucleate boiling. Jakob number, defined as lplwsatvfgJacTTh, evaluated for both the cases based on their respective fluid conditions and wall superheat of 40 K result in the following values: and For a wall superheat of 44 K, Chen et al. (2004) obtained a heat flux of ~ 0.5 MW/m2 during the bubble nucleation and growth pro cess. The average heat flux in regions 4 7256.25FCJa 52.29ArJa 136 PAGE 137 Figure 57. Variation of vapor pr essure with time for Case I Figure 58. Average temperature of regions at initial and fi nal time periods for Case I through 7 in this work is 258 MW/m2.. Thus, this provides a missing link to the heat flux values obtained during the formation of the nonevaporating film and inte rline region in bubble growth 137 PAGE 138 process, which cannot be measured via expe riments yet. Figure 510 shows a schematic suggesting the heat flux distribut ion along the nanoand microlaye rs at the base of the bubble during the nucleation and early gr owth process. The formation of the nonevaporating film and the interline region result in very high heat flux rates of 1001000 MW/m2. The heat flux rate drops down significantly as we move away from the interline region into the evaporating film region, although heat transfer occurs for a much longer time in this region. The heat flux rate drops down further into the intrinsic meniscus re gion. The overall average h eat transfer rate is the range of 0.11 MW/m2 for a wall superheat of 40 K. Moore and Mesler (1961) studied during nucleate boiling by measuring the surface temperature of a microscale area beneath the nucleating bubble. The thermocouple measured the average temperature around a circle of diameter 0.005 inch (area = m2). The fluid boiled in the experiments was water at atmos pheric pressure, and it was maintained at its saturation temperature by auxili ary heaters. The authors found that the surface temperature occasionally dropped 20oF to 30oF. They calculated the heat flowing out of the heater (from the thermocouple region) during these temperature dr ops, and obtained heat flux rates as high as 200,000 Btu/hr.ft2 (~ 6.31 MW/m2). The authors conclude that th e high heat flux rates are due to the vaporization of the microlayer at the base of the bubble when it is growing on the surface. For comparison purposes, the critical heat flux of water at atmospheric pressure is ~ 1 MW/m2. Thus, high heat fluxes for a shor t interval of time are possible from the base of the bubble, similar to what we have obtained from MD simulations. 81.26710 138 PAGE 139 Figure 59. Average heat and ev aporation flux rates for region s 110 during the total heating period for Case I on the order of 100 MW/m2 average value of the order of 0.1 MW/m2 nonevaporating film and interline region evaporating film region intrinsic meniscus region Figure 510. Heat flux distribu tion along nanoand microlayers at the base of the bubble 139 PAGE 140 Figure 511 shows the disjoining pressure vari ation along the length of the meniscus at three different time steps. During the initial equi librium period, the disjoining pressure is small owing to the thicker meniscus height. As eva poration occurs, the thickness at center of the meniscus begins to decrease resulting in increase of disjoining pressure. It increases significantly upon the formation of nonevaporating film, and at the end of the simulation period (t=2500 ps) its value in the nonevaporating film region is 4.34 MPa. The di sjoining pressure quickly goes down to nearzero values as the meniscus thic kness increases away fr om the nonevaporating film region. The capillary pressure shows an interesting trend for nanoscale meniscus in Fig. 512. At t=1500 ps, the capillary pressure is higher than the disjoining pressure at the center of the meniscus since the meniscus has increased curv ature and nonevaporating film has not formed. At the end of the simulation, the capillary pressure is zero in the nonevaporating region as expected. However, negative capilla ry pressure values are seen at the ends of the meniscus, and the absolute values are significantly higher toward s the end of the simulation than compared to during initial equilibrium. The reas on for the negative pressure is th at the liquid film at the ends of the meniscus is being pulled by the liquid in the na nochannel, and also towards the center of the meniscus due to high disjoi ning pressure. Since there is no additional feedin of atoms to negate the pulling effects in this simulation setup, pressure at the ends of the meniscus reaches negative values. Usually a negative capillary pres sure would result in cav itation in the liquid. A sample calculation in region 10, where the negative pressure is greatest, shows that for LV = 0.009 N/m and P = 1.87 MPa (averaged over th e region), the value of Rc ~ 9.63 nm which is greater than the possible capillary height of 7.96 nm (average height of interface in the region). Thus, the meniscus can exist in me tastable state in this case. 140 PAGE 141 Figure 511. Variation of disjoining pressure along me niscus length for Case I Figure 512. Variation of capillary pre ssure along meniscus length for Case I 141 PAGE 142 Results and Discussion Case II Figure 513 depicts the comput ational domain for Case II at different time intervals. Unlike Case I, all walls are at a higher temper ature and liquid argon in th e nanochannels is also heated up. The liquid in the nanochannel expands decreasing its density due to which some atoms (~ 300 atoms from each nanochannel) are pushed towards the meniscus and in turn get evaporated. Evaporation flux is high er in Case II as the meniscus is at a greater temperature than Case I. As can be seen from Fi gs. 513df, meniscus curvature is less steep and more concave in shape compared to Case I. Figure 514 shows the variation of liquid atoms for regions 5 to 10. Similar to Case I, nonevaporating film forms at the center of the meniscus in regions 5 and 6. However, the nonevaporating film is formed at a much later st age, at around 2500 ps, and is also greater in thickness compared to Case I. Evaporation also occurs from regions 9 and 10 since those regions are at a higher temperature than in Case I. The pr essure in the vapor space is calculated at every time step, and is averaged over a time period of 50 ps (10,000 time steps). The pressure increases from an initial value, averaged from 8001000 ps, of 0.153 MPa to a value of 1.168 MPa in the final 200 ps. The saturation temperature corres ponding to this final vapor pressure is Tsatpres = 119.32 K. Figure 515 shows the average temperature of each region from 8001000 ps and 28003000 ps. Temperatures of all regions of the meniscus are greater than Tsatpres and thus evaporation occurs from all meniscus regions, except the nonevaporati ng regions. As mentioned earlier, a nonevaporating thin film can exist in equilibrium on a solid surface even when the vapor pressure in the surroundi ng gas is below the normal satu ration pressure for the film temperature (Carey 1992). 142 PAGE 143 Figure 513. Snapshots of xz plane at different time steps for Case II 143 PAGE 144 Figure 514. Variation of liquid atoms with time for regions 510 for Case II The net heat flux and evaporation flux rates are shown in Fig. 516 for regions 1 through 10. The heat flux values are averaged over t1 = 8001000 ps, t2 = 15001700 ps and t3 = 27002900 ps. The overall average heat flux rates are minimum for regions 5 and 6, unlike Case I, since nonevaporating film forms in these regions wh ereas evaporation conti nues to take place in other meniscus regions due to higher temperature. High heat flux values ar e seen in the first 600 ps of evaporation period from regions near th e nanochannels due to higher temperature, and addition of atoms from the nanochannel to these regions. During the later part of the evaporation period (16002800 ps), evaporation from those re gions significantly re duce as the remaining atoms have greater resistance to evaporation due to the attractive forces by the atoms in the nanochannels. In this period, maximum heat flux occu rs from regions 3, 4, 7 and 8; formation of nonevaporating film in region 5 and 6 reduce their heat flux ra tes. Overall, as expected, maximum evaporation occurs from region 4 and 7 i.e. regions surrounding the nonevaporating film (interline regions). 144 PAGE 145 Figure 517 shows the disjoining pressure vari ation along the width of the meniscus at three different time steps for Case II. The di sjoining pressure is greater for Case I (Pd = 4.34 MPa) than Case II (Pd = 1.31 MPa) as expected. For Case I, the meniscus is pulled towards the nanochannel walls as the temperatur e is lower there compared to at the center of the meniscus. This causes the nonevaporating fi lm to be thinner at the center However, in Case II the nanochannel walls are also heated causing more atoms to flow to the meniscus, and since the meniscus is not pulled with the same intensity to wards the nanochannels like in Case I, and also due to the higher temperature, the atoms in th e meniscus have more freedom to rearrange resulting in a more uniform curvature. This incr eases the thickness of the nonevaporating film at the center of the meniscus for Case II, and hence the disjoining pressure is lower. The capillary pressure variation along the leng th of the meniscus is shown in Fig. 518 for Case II. At t=1500 ps, the capil lary pressure is higher than the disjoining pressure at the center of the meniscus since nonevaporating film has not formed. At t=2500 ps, the capillary pressure is zero in the nonevaporating region as the nonevaporating film has a flat curvature. However, negative capillary pressure values are seen, similar to Case I, at the ends of the meniscus near the nanochannels. Since the menisc us is at a higher temperature in Case II along with lower disjoining pressure at center, pressu re at the ends of the meniscus reach higher negative values for Case I (Pc ~ 2.5 MPa) compared to Case II (Pc ~ 1 MPa). To verify that cavitation cannot occur, the critical cavitation radius is calculated in region 10. LV = 0.0044685 N/m and P = 0.79 MPa (averaged over the region), the value of Rc ~ 11.31 nm which is greater than the possible capillary height of 6.61 nm (a verage height of interf ace in the region). 145 PAGE 146 Figure 515. Average temperature of regions at initial and final time periods for Case II Figure 516. Average heat and evaporation flux rates for regions 110 during the total heating period for Case II 146 PAGE 147 Figure 517. Variation of disjoining pressure along me niscus length for Case II Figure 518. Variation of capillary pre ssure along meniscus length for Case II 147 PAGE 148 Meniscus Boundary Condition As mentioned earlier, the curved meniscus was formed in this simulation by using nanochannels on either side of the meniscus. Th e liquid atoms inside th e nanochannels provide the necessary attractive force for meniscus form ation. If the nanochannels were to be removed, the meniscus would change shape to minimize its free energy and result in a flat thin film. However, using nanochannels to form the meni scus is computationally demanding and also provides physical interference in puremeniscus analysis. Hence, a meniscus boundary condition is required. The boundary condition should provide external force, which should be a function of distance, to the meniscus atoms so as to form a curvature. Two methods are proposed in this section, the first one being an analytical solution and the second derived using curvefitting from this present work. Analytical Meniscus Boundary Condition Figure 519a shows the schematic of a meni scus. The meniscus is made up of liquid atoms surrounded by vapor. Since molecular dynamic s computations are very exhaustive, taking into account the complete meniscus would be impractical. Hence, the domain is restricted between two boundaries (named 1 and 2) as shown. In order to maintain the meniscus curvature, the liquid atoms within the boundaries should still experience force from the extended meniscus regions as shown. An analytical function is de rived here to serve th is purpose. The extended meniscus region is approximated by rectangular regions as shown. The force exerted by a rectangular region on an atom, say atom P, is evaluated as shown in Fig. 519b. Atom P has coordinates x0, y0, z0. The rectangular re gion is of height t thickness x infinite in ydirection, and its center is at coordinates x1, y1 and z1. Consider a infinitesimal element in the recta ngular region as shown of constant thickness x width dw and height dh. Let x = x1x0, y = y1y0 and z = z1z0. Using LennardJones potential for 148 PAGE 149 argonargon interaction, the energy of interaction, dU between P and the infinitesimal element is: 1264avN dU xdwdh rrM (5.13) where avNM2 is the number of atoms per unit volume in the rectangular region, and r is the defined as: 22 2 101010 22 2rxxywyzhz xywzh (5.14) Substituting in equation 5.13 and integrating as follows: 12 6 22 22 0 224t avN U x dwdh M xywzhxywzh 11 5 12 6 22 0 2263 131 4 256 8t avN Ux M xzhxzh d h 12 6 1263 3 4( ) 256 8avN Uxf x zf M ( )x z (5.15) where, 82644628 1 9/2 22100.00317463158401008576128 (,)hzt hzhxhxhxhxh fxz xhx (5.16) 22 2 3/2 22432 (,) 3hzt hzhxh fxz xhx (5.17) In order to obtain the net force acting on atom P from the rectangular region in x, y and z directions, the energy of interaction U has to be differentiated as follows: 149 PAGE 150 ;;xyzdUdUdU FFF dxdydz (5.18) As it can be seen from the above derivation, the final analytical solution turns out to be complicated, cumbersome and computationall y demanding to implement in a molecular dynamics simulation, and all combinations of atoms and rectangular regions will have to be considered at each time step. A simpler approx imation is presented in the next section. Figure 519. (a) Schematic of a meniscus to compute forces acting by the extended meniscus region on liquid atoms in the bounded meniscus region, and (b) configuration for force interaction between an atom and a rectangular region CurveFitted Meniscus Boundary Condition In this present work, the meniscus was formed using the presence of nanochannels as shown in schematic of Fig. 520. During initial equilibrium, the meniscus was divided into equal 150 PAGE 151 length slices of 0.5 (shown in Fig. 12 as 1, 2, 3, 4, 5 N) and the liquid atoms in the meniscus were assigned to their respect ive slices. The force exerted by the nanochannel, which includes atoms of argon as well as platinum wall, on al l the atoms present in each slice was computed using molecular dynamics simulation in x, y and zdirections. Figure 521 shows the variation of forces acting on the slices as a functi on of distance from the nanochannel. Figure 520. Schematic showing meniscus and nanochannels used in the present work to determine the force distribution along the meniscus length by the nanochannels Force in the y direction is nearly zero, as e xpected, since the boundary conditions are periodic. Forces in x and z directions follow a trend si milar to LennardJones potential. However, the force curves here are significantly steeper. A function similar to the LennardJones potential,B x D F AxCx2x2, is used for curve fitting. This external force can be applied as a meniscus boundary condition for a LennardJones fl uid. The following curve fits, as shown in Figs. 521a and 521b, ma tch the data points: 32 30.10421.179;0;0.023650.2144xy zFxxFFx (5.19) where x is in nm and F is in 1014 N. A general force function can be derived based on the above trends. Thus, a meniscus can be formed by applying a force, 3n F AnCn (n = x, y or z; n is in nm; F is in 1014 N), at the boundaries of the fluid. Th e coefficients A and C can be chosen 151 PAGE 152 based on the strength and depth of the force requir ed, which will in turn affect the curvature of the meniscus. Figure 521. Force distribution on the meniscus by the nanochannel as a function of the distance from the nanochannel (a) Fx and Fy, (b) Fz 152 PAGE 153 CHAPTER 6 CONCLUSION AND FUTURE WORK Conclusion Proposed Thermal Wall Model A novel fluidheated wall thermal equilibri um model for numerically simulating heating/cooling of fluid atoms by wall atoms is proposed. In order to validate the proposed model, the thermal conductivity of argon and the change in internal energy of the system are obtaine d from simulation and compare well to those from the thermodynamic properties of argon. As the two properties evaluate d are statedependent only, ad ditional validation cases are run to test the timedependent nature of the model by pl acing liquid argon between two platinum walls with simultaneous heating/cooling of the walls. The temperature gradient, along the height of the channel, is evaluated and compared to the analytical solution obtained by solving 1D heat equation and excellent agreement was found between the predicte d and theoretical results. Additional simulations are c onducted with heating at both walls by increasing the wall temperatures to 120 K for two different nanocha nnel heights. It is concluded that heat transfer occurs at a faster rate than pr edicted by the heat equation at nanoscale. The error between the molecular dynamics simu lation and the heat equations analytical result depends on the length scale of the problem and reduces with the increase in channel height. The results compared well in the previous case (with heating and cooling at the walls) since higher heating rate at th e upper wall is matched with higher cooling rate at the lower wall. Th ese comparisons show the physical soundness of the proposed model. The Maxwell speed distribution (MSD) was used to provide another verification of the proposed fluidwall thermal equilibrium model. At high temperatures argon vapor would behave like an ideal gas. The proposed heat transfer model is responsible for imparting kinetic energy to the argon atoms and evaporation. The simulation results of an evaporating argon thin film on platinum su rface matched well with the MSD function; this verifies that the veloci ties imparted by the proposed fluidwall heat transfer model follow expected physics and at high temperatures argon vapor behaves like an ideal gas. 153 PAGE 154 Transient Film Evaporation The fluidwall thermal equilibrium model is used to capture the nanoscale physics of transient phase transition of a thin liquid ar gon film on a heated platinum surface and the eventual colloidal adsorption phenomenon as the evaporation is dimini shing using molecular dynamics. The objective was to provide microsc opic characterizations of the dynamic thermal energy transport mechanisms during the liquid film evaporation and also the resulting nonevaporable colloidal adsorbed liquid layer at the end of the evaporation process. The simulation domain is constructed of platinum wall atoms with argon as the working fluid. Phase change process is studied by simulating evapor ation of a thin liquid argon film on the platinum wall. High heat transfer and evaporation rates are obtained. Gradual vaporization of the film takes pl ace exponentially decreasing with time and a nonevaporating film is obtained. A Hamakertype analysis is done for the nonevaporating film and the Hamaker constant is obtained using molecular dynamics whic h has not been attempted before using molecular dynamics simulation. The value of the Hamakertype constant falls in the typical range. This analysis is quantifies the nonevapor ating film and makes an attempt to use molecular dynamics simulation results in continuum mechanics. Simulation of Seven Cases in Nanochannel Film Evaporation The transient film evaporation simulation is extended further to gain better understanding on the occurrence of the nonevapor ating film, and the parameters which govern its formation. Seven cases were studied to determine the effect of variation of nanochanne l height and variation of liquid film thickness on heat tran sfer rate and evaporation rate. Nonevaporating film was obtained in cases 16. The starting heat transfer rate depends on the evaporation rate (due to phase change) and the thickness of liquid film (more at oms implies more energy stored). 154 PAGE 155 Net heat flux increases with increase in bot h nanochannel height and liquid film thickess, and et evaporation flux increases with an increase in nanochannel height. However, an optimum film thickness exists for a fixed nanochannel height as total number of liquid atoms & vapor vo lume compete with each other. Nonevaporating film thickness, liquidvapor equimolar position, vapor pressure and Hamaker constant of the nonevaporating films were evaluated. The value of Hamaker constant increases with an increase in vapor pressure. Vapor pressure was plotted against vapor dens ity for initial and final states, and all the states were found to lie on the vapor saturation curve confirming the occurrence of nonevaporating film. Case 7 was designed to simulate evaporati on from the upper wall of the nanochannel and condensation on the lower wall. All liquid f ilm from the upper wall was evaporated and no nonevaporating film was obtained. Study of all cases show that pressure and temperature are the dominant f actors governing the formation of nonevaporating film, wher eas wallliquid molecula r interaction is a lesser factor. Wallliquid molecular interacti on only governs the wettability of surface, and a nonevaporating film may or may not form on a specifi c wallliquid pair depending on the vapor pressure and wall temperature. NanoDroplet Impact on a Homogeneous Surface The dynamics of nanodroplet impact on a homogeneous surface was also simulated using molecular dynamics. The wettability and te mperature of the surfac e was varied to study four cases. An elliptical interfacial fit was obt ained from the computed interface markers with the value of Contact angle made by the drop w ith the wall was obtained using the interfacial fit. 20.999 R The nanodrops accelerate significantly on reaching close to the wall and the actual impact velocity in higher than that expected to reach by the application of the external force. Complete spreading was obtained for the high wettable surface at 90 K, whereas drop after impact on the low wettable su rface at 90 K oscillated about a mean contact angle of 126.7o. 155 PAGE 156 Liedenfrost phenomenon was observed at wall temperature of 300 K for both cases of wettability, although it was seen that the evaporation rate reduced significantly for the low wettable surface. Meniscus Evaporation at Nanoscale The nonevaporating region and the in terline region of a threephase contact line are on the order of nanometers, and hence the heat and mass transf er characteristics are studied using molecular dynamics. A concave meniscus was formed by placing liquid argon between a lower and an upper platinum wall, with an opening in the upper wall. After reaching initial equilibrium, the liquid meniscus was heated by the platinum wa ll and evaporation was observed. Two cases are studied where in one only a part of the platinum wall is increased in te mperature while in the other all the walls temperatures are elevated to obtain evapor ation. The liquid meniscus was divided into 10 equal regions, and variation of temperature, energy, total number of liquid atoms in each region was studied over the period of simu lation. A nonevaporating film was formed at the center of the meniscus. Heat and mass transfer continued to occur from the interline regions (regions surrounding the nonevapora ting film). Very high heat flux rate of the order of 100 MW/m2 and high evaporation rates of the order of 1000 kg/m2s are achievable from nanoscale evaporating meniscus. Capillary and disjoining pressures of the liquid meniscus were also computed. The disjoining pressure increased significantly after the formation of the nonevaporating film. High negative capillary pressures were obtained towards the end of the simulation. The reason for the negative pressure is th at the liquid film at the ends of the meniscus is being pulled by the liquid in the nanochannel, and also towards the center of the meniscus due to high disjoining pressure. It is shown that cavitation cannot occur as the capillary size is smaller than the critical cavitation radius, and th us the meniscus can exist in metastable state. Negative capillary pressure is not uncommon, and has been experimentally measured by researchers. A meniscus boundary condition was developed using analytical and curvefitting 156 PAGE 157 methods. The analytical solution turns out to be complicated, cu mbersome and computationally demanding to implement in a molecular dynamics simulation. A simpler approximation, using curvefitting over the statistical data computed from this work, leads to a force function of the form which can be applied at the boundari es of a liquid film to create a meniscus. 3 nFAnCn2 Future Work For the future work, we will follow the mu ltiscale approach with an innovative threeregion model (molecular dynamics, coarsegrained continuum) to ensure a smooth matching between the nanoscale and the continuum. As sh own in Figure 61, a multiscale threeregion model is proposed where an additional micros cale region will be sandw iched between the nano region and the continuum macro region to ensure a smooth transition for accurate results. The threeregion model will be applied in the liqui d region only, whereas direct nano to continuum matching will be done in the vapor region. The re ason is that the density of vapor region is significantly less compared to liquid, and creatin g coarsegrained model would not be practical. The numerical scheme for each of the three regions is explained next. Molecular Dynamics Simulation (Nano Region) As explained in Chapter 2, all MD simulations follow Newtons second la w for simple atomic systems or in a more generalized form (Lagrangi an/Eulerian formulation) for complex geometry systems. The force on each molecule by all other molecules in the system can be determined by using the molecular interaction potential function. From this force, the equation of motion can be integrated over time to obtain the new positi ons and velocities of each molecule. This information can be statistically analyzed by space/time averaging to provide macroscopic physical properties, and t hus the evolution of the system can be studied. 157 PAGE 158 Figure 61. Schematic showing multiscal e method of solution for proposed work 158 PAGE 159 CoarseGrained Method (Micro Region) The coarsegrained (CG) approach is basically st ill a discrete particle method where the particles are not the actual atoms or molecu les, instead, each CG particle (a quasimolecule) represents a group of real particles. The number of real par ticles represented by a CG particle depends on the coarseness level specified. After the details in an actual atomistic description are reduced, the accessible timeand lengthscales of molecula r simulations can be increased. CG models involve two steps: the firs t step is to group the degrees of free dom (molecules) in the system into fewer structural units of CG par ticles. The second step is the c onstruction of an effective force field to describe the interac tion between the CG particles. The forcematching (FM) method (Chu et al. 2006) will be used to derive equivale nt CG force fields directly from the underlying atomistic forces. The FM procedure applied to the CG trajectories and forces will yield the effective interaction between CG particles as is present in the underlying atomistic simulation. The coarsegrain concept has been extensive re searched and developed by the group led by Prof. S. Izvekov. A full review of this met hod is given in Izvekov and Voth (2005). Finite Volume Method (Macro Continuum Region) The numerical algorithm of the finite volume method (FVM), a wi dely used numerical solution technique in computational fl uid dynamics (CFD) for a continuum domain, consists of the following steps: Formal integration of the governing differential equations of fluid flow and heat transfer over all the finite control vol umes of the solution domain Discretization involves the substitution of a variety of finitedifferencetype approximations for the terms in the integrated equation representing flow processes such as convection, diffusion and sources. This convert s the integral equations into a system of algebraic equations. Solution of the algebraic equations (matrix inversion) by an iterative method 159 PAGE 160 The first step, the control volume integration, distinguishes the finite volume method from all other CFD techniques. The resulting statements express the exact cons ervation of relevant properties for each finite size cell. This clear relationship between the numerical algorithm and the underlying physical conservation principle forms makes the con cepts of FVM more effective and accurate than finitedifference, finite element and spectral methods. MultiScale Model and Matching Procedure Figure 61 above shows a schematic of the multis cale model based on the above methods applied to the bubble growth and departur e dynamics. Figure 61(a) repres ents a typical bubble on a flat heater surface. Figure 61(b) depicts the zoomedin region at the threephase contact line and showcases the multiscale nature of the problem. It has been widely accepted that the threephase contact line region is constituted of: nonevaporating adsorbed thin film region, evaporating thinfilm region (micro layer) and bulk liquid regi on (macro layer). It is important to model the regions at their respective lengthscales to obtain accurate solutions and understand the underlying physics. Figure 61(c) shows the MD, AACG and FVM met hods implementation in the proposed work. MD simulation and MDCG matching (nanomicro matching) The matching between the nano and micro regions will be accomplished by a very attractive approach that is based on creating a mixed allatom and coarseg rained particle (AACG) model (Shi et al. 2006), in which the most interesting pa rt of the system is represented in full atomistic detail, and the remaining parts are modeled at the CG level. The mixe d AACG mode l will be used and implemented in the proposed work for ma tching the two parts. The interactions in the mixed AACG model can be cla ssified into: atomatom, atomCG and CGCG interactions. Forces on all atoms and CG sites are given as: 160 PAGE 161 ,atom atomatom atomCG ii j jatomsji jCGFFF i j; ,CG atomCG CGCG ii ji j jatoms jCGjiFFF As explained in the previous section, CG particles will be constructed and their force fields will be derived using the FM method. Figure 61(c) shows the CG particles as large dots. The AACG method accounts for all atoma tom, atomCG and CGCG inter actions, and thus MD and CG lengthand timescales will be matched intuitively. CG and FVM matching (micromacro matching) Domain decomposition technique will be used to match the CGFVM for liquid domain, and MDFVM for vapor domain. The schematic of the technique is shown in Figure 62. Initially a steady state solution is assumed in eith er of the domain, say in CG domain, CG. The solution in CG becomes the overlapped boundary c ondition for the continuum domain, contIBC. Now, we obtain a solution in cont using contIBC and nonoverlapped, contEBC. The resulting solution in cont gives a new boundary condition CGIBC for the CG domain. A new solution in CG is obtained using the new boundary conditions. These st eps are repeated to march forward in time and study the system evolution. The key element of this method is the scheme for coupling the CG and continuum solutions to ensure continuity of fluxes across the overlap region. Liu et al. (2007) used the same scheme to couple MDconti nuum simulations to stu dy heat transfer in a fluid layer confined between flat and rough su rfaces. Their results matched well with pure MD results. Our proposed work differs in the way that we have added an additional CG model between the MD and continuum levels to ensu re smoother transitioning between the nano and micro length scales in liquid domain. This w ill lead to more accurate results, as well as application to longer length scales in liquids as density is higher. Constrained dynamics will be used to couple the CG and c ontinuum velocities in momentum equations, and temperatures and 161 PAGE 162 energy will be coupled by applying velocity rescaling. Mass fluxes will be coupled by adding/removing CG particles. As can be seen in Figure 61, blue boxe s depict the staggered grids for continuum simulation. The propertie s will be averaged from the CG sites: ux defined at crosses, uy defined at triangles, and p and T defined at circles. These matching techniques will ensure coupling at different time and lengthscales, and passing of information between them while ensuring conservation of mass, momentum and energy. Overlap Region CG/MD Domain Continuum DomaincontCGcontEBCCG/MDEBCCG/MDIBCcontIBCCG/MD = CG/MDIBC CG/MDEBCcont = contIBC contEBCEBC External boundary condition IBC Iterative boundary condition Figure 62. Domain Decomposition method 162 PAGE 163 LIST OF REFERENCES Abraham FF. 1978. The interfacial density profile of a LennardJones fluid in contact with a (100) lennardjones wall and its relationship to idealized flui d/wall systems: A monte carlo simulation. J. Chem. Phys. 68 : 37136 Ado MH, de Ruijter M, Vou M, De Coninc k J. 1999. Droplet spreading on heterogeneous substrates using molecular dynamics. Phys Rev E. 59 : 746 Allen MP, Tildesley DJ. 1987. Computer simulation of liquids, Oxford England; New York: Clarendon Press; Oxford Un iversity Press. 385pp Angell CA. 1988. Approaching the limits. Nature. 331 : 2067 Blake TD. 2006. The physics of moving wetting lines. J. Colloid Interface Sci. 299 : 113 Briggs LJ. 1955. Maximum superheating of water as a measure of negative pressure. J. Appl. Phys. 26 : 10013 Bussmann M, Mostaghimi J, Chandra S. 1999. On a threedimensional volume tracking model of droplet impact. Phys. Fluids. 11 : 140617 Butt H, Cappella B, Kappl M. 2005. Force measurements with the atomic force microscope: Technique, interpretation and applications. Surface Science Reports. 59 : 1152 Carey VP. 1992. Liquidvapor phasechange phenomena, Taylor and Francis. 672pp Chen T, Chung JN. 2003a. Heattransfer effects of coalescence of bubbles from various site distributions. Proceedings: Mathematical, Phys ical and Engineering Sciences. 459 : 2497527 Chen T, Chung JN. 2002. Coalescence of bubbl es in nucleate boiling on microheaters. Int. J. Heat Mass Transfer. 45 : 232941 Chen T, Chung JN. 2003b. An experimental study of miniaturescale pool boiling. J. Heat Transfer. 125 : 107486 Chen T, Klausner JF, Chung JN. 2004. Subcooled bo iling heat transfer and dryout on a constant temperature microheater. Int J Heat Fluid Flow. 25 : 27487 Chen W, Fish J. 2006. A generalized spacetime mathematical homogenization theory for bridging atomistic and continuum scales. Int J Numer Methods Eng. 67 : 25371 Chen X, Wu NJ, Smith L, Ignatiev A. 2004. Th infilm heterostructure solid oxide fuel cells. Appl. Phys. Lett. 84 : 27002 Chu J, Izveko S, Voth GA. 2006. The multiscale challenge for biomolecular systems: Coarsegrained modeling. Molecular Simulation. 32 : 211 163 PAGE 164 Daiguji H, Hihara E. 1999. Molecular dynamics study of the liquidva por interface of lithium bromide aqueous solutions. Heat and Mass Transfer. 35 : 2139 Davidson MR. 2000. Boundary integral prediction of the spreading of an inviscid drop impacting on a solid surface. Chemical Engineering Science. 55 : 115970 de Andrade J, Stassen H. 2004. Molecular dynamics studies of th ermal conductivity time correlation functions. Journal of Molecular Liquids. 110 : 16976 DelgadoBuscalioni R, Cove ney PV. 2003. Continuumparticle hybrid coupling for mass, momentum, and energy transfer s in unsteady fluid flow. Phys Rev E. 67 : 046704 Dhir VK, Liaw SP. 1989. Framework for a unifi ed model for nucleate and transition pool boiling. ASME J. Heat Transfer. 111 : 73946 Dhir VK. 2006. Mechanistic prediction of nucle ate boiling heat transferachievable or a hopeless task? J. Heat Transfer. 128 : 112 Dhir VK. 2001. Numerical simulations of poolboiling heat transfer. AICHE J. 47 : 81334 Ding H, Spelt PDM. 2007. Inertial effects in dr oplet spreading: A comp arison between diffuseinterface and levelset simulations. J. Fluid Mech. 576 : 28796 Ding H, Spelt PDM. 2007. Wetting condition in di ffuse interface simulations of contact line motion. Phys Rev E. 75 : 046708 Drazer G, Khusid B, Koplik J, Acrivos A. 2005. Wetting and particle adsorption in nanoflows. Phys. Fluids. 17 : 017102 Dussan V. EB. 1976. The moving contact line: The slip boundary condition. Journal of Fluid Mechanics Digital Archive. 77 : 66584 Dussan V. EB, Davis SH. 1974. On the motion of a fluidfluid interface along a solid surface. Journal of Fluid Mechanics Digital Archive. 65 : 7195 Dussan EB. 1979. On the spreading of liquids on solid surfaces: Static and dynamic contact lines. Annu. Rev. Fluid Mech. 11 : 371400 Fedorchenko AI, Wang A, Wang Y. 2005. Effect of capillary and viscous forces on spreading of a liquid drop impinging on a solid surface. Phys. Fluids. 17 : 093104 Fisher JC. 1948. The fracture of liquids. J. Appl. Phys. 19 : 10627 Flekkoy EG, Wagner G, Feder J. 2000. Hybrid model for combined pa rticle and continuum dynamics. EPL (Europhysics Letters). 52 : 2716 Francois M, Shyy W. 2003. Co mputations of drop dynamics with the immersed boundary method, part 2. Drop imp act and heat transfer. Numerical Heat Transfer, Part B: Fundamentals: An Interna tional Journal of Com putation and Methodology. 44 : 119 164 PAGE 165 Fuentes J, Cerro RL. 2007. Surface forces and inertial effects on moving contact lines. Chemical Engineering Science. 62 : 323141 Fukai J, Zhao Z, Poulikakos D, Megaridis CM Miyatake O. 1993. Modeling of the deformation of a liquid droplet impinging upon a flat surface. Phys. Fluids A. 5 : 258899 Fukai J, Tanaka M, Miyatake O. 1998. Maximum spreading of liquid droplet s colliding with flat surfaces. J. Chem. Eng. Japan. 31 : 45661 Genske P, Stephan K. 2006. Numeri cal simulation of heat transfer during growth of single vapor bubbles in nucleate boiling. International Journal of Thermal Sciences. 45 : 299309 Gentner F, Rioboo R, Baland JP, De Coninck J. 2004. Low inertia impact dynamics for nanodrops. Langmuir. 20 : 474855 Gunjal PR, Ranade VV, Cha udhari RV. 2005. Dynamics of drop impact on solid surface: Experiments and VOF simulations. AICHE J. 51 : 5978 Hadjiconstantinou NG. 1999. Hybrid Atomistic Continuum formulations and the moving contactline problem. Journal of Comput ational Physics. 154 : 24565 Hamaker HC. 1937. The Londonvan der waals at traction between spherical particles. Physica. 4 : 105872 Han DG, Choi GM. 1998. Computer simulation of the electrical c onductivity of composites: The effect of geometrical arrangement. Solid State Ionics. 106 : 7187 Hanasaki I, Nakatani A, Kita gawa H. 2004. Molecular dynamics study of ar flow and he flow inside carbon nanotube junction as a molecular nozzle and diffuser. Science and Technology of Advanced Materials. 5 : 10713 Haramura Y, Katto Y. 1983. A new hydrodynamic model of critical heat flux, applicable widely to both pool and forced convection boiling on submerged bodies in saturated liquids. Int. J. Heat Mass Transfer. 26 : 38999 Hirschfelder JO. 1954. Molecular theory of gases and liquids, New York: Wiley. 1249pp Huh C, Scriven LE. 1971. Hydrodynamic model of steady movement of a solid/liquid/fluid contact line. J. Colloid Interface Sci. 35 : 85101 Ingebrigtsen T, Toxvaerd S. 2007. Contact angles of lennardjone s liquids and droplets on planar surfaces. The Journal of Physical Chemistry C. 111 : 851823 Israelachvili J, Israelachvili JN. 1992. Intermolecular and surface forces, Academic Press London. 450pp. 2nd ed. Izvekov S, Voth GA. 2005. A multiscale coarsegraining method for biomolecular systems. The Journal of Physical Chemistry B. 109 : 246973 165 PAGE 166 Johnston HL, Grilly ER. 1942. Viscosities of carbon monoxide, helium, neon, and argon between 80 and 300K. coefficients of viscosity. J. Phys. Chem. 46 : 94863 Kandlikar SG, Kuan WK, Mukherjee A. 2005. E xperimental study of heat transfer in an evaporating meniscus on a moving heated surface. J.Heat Transfer. 127 : 24452 Kimura T, Maruyama S. 2003. A molecular dynamics simulation of water droplet in contact with a platinum surface. Presented at 6th ASME/JSME Thermal Engg. Conf. Hawaii, 2003, A1183. Klein SA. 2007. EES manual for windows, Fchart software, Koplik J, Banavar JR, Willemsen JF. 1989. Molecu lar dynamics of fluid flow at solid surfaces. Physics of Fluids A. 1 : 781 Koplik J, Banavar JR, Willemsen JF. 1988. Mole cular dynamics of poiseuille flow and moving contact lines. Phys. Rev. Lett. 60 : 1282 Liao J, Mei R, Klausner JF. 2004. The influe nce of the bulk liquid thermal boundary layer on saturated nucleate boiling. Int J Heat Fluid Flow. 25 : 196208 Lienhard JH, Dhir VK. 1973. Extended hydrodynamic theory of the peak and minimum pool boiling heat fluxes. Rep. CR2270, NASA, Liu J, Chen S, Nie X, Robbins MO. 2007. A conti nuumatomistic simulation of heat transfer in microand nanoflows. Journal of Comput ational Physics. 227 : 27991 Lunkad SF, Buwa VV, Ni gam KDP. 2007. Numerical simulations of drop impact and spreading on horizontal and inclined surfaces. Chemical Engineering Science. 62 : 721424 Mareschal M, Baus M, Lovett R. 1997. The local pressure in a cylindrical liquidvapor interface: A simulation study. J. Chem. Phys. 106 : 64554 Markvoort AJ, Hilbers PAJ, Ned ea SV. 2005. Molecular dynamics study of the influence of wallgas interactions on h eat flow in nanochannels. Phys. Rev. E. 71 : 066702 Martn ME, Aguilar MA, Chalmet S, RuizLpez MF. 2002. An iterative pr ocedure to determine lennardjones parameters for th eir use in quantum mechanics/molecular mechanics liquid state simulations. Chem. Phys. 284 : 60714 Maruyama S, Kimura T. 2000. A molecular dynam ics simulation of a bubble nucleation on solid surface. Int. J. Heat Tech. 18 : 6974 Maruyama S, Kimura T. 1999. A study on therma l resistance over a solidliquid interface by the molecular dynamics method. Therm. Sci. Eng. 7 : 63 166 PAGE 167 Maruyama S, Kurashige T, Matsumoto S, Yamaguchi Y, Kimura T. 1 February 1998. Liquid droplet in contact with a solid surface. Microscale Thermophysical Engineering. 2 : 49,62(14) Michels A, Wijker H, Wijker H. 1949. Isothe rms of argon between 0c and 150c and pressures up to 2900 atmospheres. Physica. 15 : 62733 Moffatt HK. 1964. Viscous and resistive eddies near a sharp corner. Journal of Fluid Mechanics Digital Archive. 18 : 118 Moore FD, Mesler RB. 1961. The measurement of rapid surface temperature fluctuations during nucleate boiling of water. AICHE J. 7 : 6204 Mukherjee A, Dhir VK. 2004. Study of lateral merger of vapor bubbles during nucleate pool boiling. J. Heat Transfer. 126 : 102339 Mukherjee A, Kandlikar SG. 2007. Numerical study of single bubbles with dynamic contact angle during nucleate pool boiling. Int. J. Heat Mass Transfer. 50 : 12738 Mukherjee A, Kandlikar SG. 2006. Numerical study of an eva porating meniscus on a moving heated surface. J. Heat Transfer. 128 : 128592 Nagayama G, Cheng P. 2004. Effects of interf ace wettability on microscale flow by molecular dynamics simulation. International Journal of Heat and Mass Transfer,. 47 : 50113 Nagayama G, Tsuruta T, Che ng P. 2006. Molecular dynamics simulation on bubble formation in a nanochannel. Int. J. Heat Mass Transfer. 49 : 443743 Nie X, Robbins MO, Chen S. 2006. Resolving singular forces in cav ity flow: Multiscale modeling from atomic to millimeter scales. Phys. Rev. Lett. 96 : 134501 Nijmeijer MJP, Bakker AF, Bruin C, Sikkenk JH. 1988. A molecular dynamics simulation of the lennardjones liquidvapor interface. J. Chem. Phys. 89 : 378992 NIST chemistry webbook. http://webbook.nist.gov/chemistry/fluid/ Nosonovsky M, Bhushan B. 2008. Phase behavior of capillary bridges: Towards nanoscale water phase diagram. Physical chemistry chemical physics : PCCP. 10 : Nukiyama S. 1934. The maximum and minimum values of heat transmittal from metal to boiling water under atmospheric pressure. J. Japan Soc. Mech. Eng. 37 : 367 O'Connell ST, Thompson PA. 1995. Molecular dynamicsconti nuum hybrid computations: A tool for studying complex fluid flows. Phys Rev E. 52 : R5792 Ohara T, Suzuki D. 2000. Intermolecular en ergy transfer at a solidliquid interface. Nanoscale and Microscale Thermophysical Engineering. 4 : 189 167 PAGE 168 Panchamgam SS, Chatterjee A, Plawsky JL, Wa yner Jr. PC. 2008. Comprehensive experimental and theoretical study of fluid flow and heat tr ansfer in a microscopic evaporating meniscus in a miniature heat exchanger. Int. J. Heat Mass Transfer. 51 : 536879 Park H, Carr WW, Zhu J, Morris JF. 2003. Single drop impaction on a solid surface. AICHE J. 49 : 246171 Park SH, Weng JG, Tien CL. 2001. A mol ecular dynamics study on surface tension of microbubbles. Int. J. Heat Mass Transfer. 44 : 184956 PasandidehFard M, Qiao YM, Chandra S, Mosta ghimi J. 1996. Capillary effects during droplet impact on a solid surface. Phys. Fluids. 8 : 6509 Phan CM, Nguyen AV, Evans GM. 2006. Combini ng hydrodynamics and molecular kinetics to predict dewetting between a sm all bubble and a solid surface. J. Colloid Interface Sci. 296 : 66976 Priezjev NV. 2007. Ratedependent slip boundary conditions for simple fluids. Phys. Rev. E. 75 : 051605 Qian T, Wang X, Sheng P. 2006. A variational approach to moving contact line hydrodynamics. J. Fluid Mech. 564 : 33360 Rapaport DC. 1995. The art of molecular dynamics simulation, Cambridge ; New York: Cambridge University Press. 400pp Reznik SN, Yarin AL. 2002. Spreading of an ax isymmetric viscous drop due to gravity and capillarity on a dry horizontal wall. Int. J. Multiphase Flow. 28 : 143757 Sadasivan P, Unal C, Nelson R. 1995. Perspectiv e: Issues in CHF modelingthe need for new experiments. J. Heat Transfer. 117 : 55867 Sadus RJ. 1999. Molecular simulation of fluids, Amsterdam: Elsevier Science. 552pp Sakashita H, Kumada T. 19930815. A new model for CHF in pool boiling at higher pressure. JSME international journal.Ser.B, Fluids and thermal engineering. 36 : 4228 Schoen PAE, Poulikakos D, Arcidiacono S. 2005. Ph ase change of a confined subcooled simple liquid in a nanoscale cavity. Phys Rev E. 71 : 041602 Shi Q, Izvekov S, Voth G. 2006. Mixed atom istic and coarsegrained molecular dynamics: Simulation of a membranebound ion channel. The Journal of Physical Chemistry B. 110 : 150458 Shikhmurzaev YD. 2006. Singularities at the m oving contact line. mathematical, physical and computational aspects. Physica D. 217 : 12133 168 PAGE 169 Sikalo S, Wilhelm H, Roisman IV, Jakirlic S, Tropea C. 2005. Dynamic contact angle of spreading droplets: Expe riments and simulations. Phys. Fluids. 17 : 062103 Son G. 2001. A numerical method fo r bubble motion with phase change. Numerical Heat Transfer, Part B: Fundamentals: An International Journal of Computation and Methodology. 39 : 509 Spijker P, ten Eikelder, H. M. M., Markvoort AJ, Nedea SV, Hilb ers PAJ. 2008. Implicit particle wall boundary condition in molecular dynamics. Proceedings of the Institution of Mechanical Engineers Part C Journal of Mechanical Engineering Science. 222 : 85564 Stoddard SD, Ford J. 1973. Numerical experiments on the stochastic behavi or of a lennardjones gas system. Phys. Rev. A. 8 : 1504 Stokes RJ, Evans DF. 1996. Fundamentals of interfacial engineering, New York: VCH Publishers. 701pp Sutmann G. 2002. Classical molecular dynamics, quantum simulations of complex manybody sustems: From theory to algorithms lecture notes. NIC Series. Rep. 10, Tas NR, Mela P, Kramer T, Berenschot JW, van dB. 2003. Capillarity induced negative pressure of water plugs in nanochannels. Nano Letters. 3 : 153740 Thompson PA, Brinckerhoff WB, Robbins MO. 1993. Microscopic studies of static and dynamic contact angles. J. Adhes. Sci. Technol. 7 : 535,554(20) Thompson PA, Robbins MO. 1989. Simulations of contactline motion: Slip and the dynamic contact angle. Phys. Rev. Lett. 63 : 766 Thompson PA, Troian SM. 1997. A general boundary condition for liquid flow at solid surfaces. Nature. 389 : 3602 Tomar G, Biswas G, Sharma A, Agrawal A. 2005. Numerical simulati on of bubble growth in film boiling using a coupled levelset and volumeoffluid method. Phys. Fluids. 17 : 112103 Tyree MT. 2003. Plant hydrau lics: The ascent of water. Nature. 423 : 923Voronov RS, Papavassiliou DV, Lee LL. 2006. Boundary slip and wetting properties of interfaces: Correlation of the contact angle with the slip length. J. Chem. Phys. 124 : 204701 Wagner G, Flekkoy E, Feder J, Jossang T. 1 August 2002. Coupling molecular dynamics and continuum dynamics. Comput. Phys. Commun. 147 : 670,673(4) Wayner Jr. PC. 1999. Intermolecular forces in phasechange heat tr ansfer: 1998 kern award review. AICHE J. 45 : 205568 169 PAGE 170 Wemhoff AP, Carey VP. 2005. Mole cular dynamics exploration of thin liquid films on solid surfaces. 1. monatomic fluid films. Nanoscale and Microscale Thermophysical Engineering. 9 : 331 Weng J, Park S, Lukes JR, Tien C. 2000. Molecu lar dynamics investigation of thickness effect on liquid films. J. Chem. Phys. 113 : 591723 Wu CYW. 2003. A molecular dynamics simulatio n of bubble nucleation in homogeneous liquid under heating with constant mean negative pressure. Nanoscale and Microscale Thermophysical Engineering. 7 : 137 Wu J, Dhir VK, Qian J. 2007. Numerical simu lation of subcooled nuc leate boiling by coupling levelset method with movingmesh method. Numerical Heat Transfer, Part B: Fundamentals: An Interna tional Journal of Com putation and Methodology. 51 : 535 Wu YW, Pan C. 2006. Molecular dynamics simulati on of thin film evaporation of LennardJones liquid. Nanoscale and Microscale Th ermophysical Engineering. 10 : 157 Xu JL, Zhou ZQ. 2004. Molecula r dynamics simulation of li quid argon flow at platinum surfaces. Heat and Mass Transfer. 40 : 85969 Xu J, Li Y. 2007. Boundary conditions at th e solidliquid surface over the multiscale channel size from nanometer to micron. International Journal of Heat and Mass Transfer,. 50 : 257181 Yamamoto M, Kinoshita M, Kaki uchi T. 2000. Structure of the pt(111)/liquid interface: A firstprinciples/RHNC calculation. Electrochim. Acta. 46 : 16574 Yang S. 2006. Effects of surface roughness and interface wettability on nanoscale flow in a nanochannel. Microfluidics and Nanofluidics. 2 : 50111 Yang SH, Nosonovsky M, Zhang H, Chung K. 2008. Nanoscale water capillary bridges under deeply negative pressure. Chemical Physics Letters. 451 : 8892 Yi P, Poulikakos D, Walther J, Yadigarogl u G. 2002. Molecular dynamics simulation of vaporization of an ultrathin liquid argon layer on a surface. International Journal of Heat and Mass Transfer. 45 : 2087100 Zhang ZM. 2007. Nano/microscale heat transfer, New York, NY: McGrawHill. 479pp Ziarani AS, Mohamad AA. 2008. Effect of wall roughne ss on the slip of fluid in a microchannel. Nanoscale and Microscale Thermophysical Engineering. 12 : 154 Zuber N. 1959. Hydrodynamic aspects of boiling heat transfer. Rep. AEC Report No. AECU4439 170 PAGE 171 BIOGRAPHICAL SKETCH Shalabh Chandra Maroo was born in 1982, in Bh ilai, India. He attended high school at Deepika E. M. School in Rourkela, India and graduated in 1999. He competed in Indian Institute of TechnologyJoint Entrance Examinations and was selected for b achelors program in mechanical engineering at Indian Institute of Technology (IIT), Bombay, India. He received his Bachelor of Technology degree from IIT Bomb ay in 2003. Following graduation, he was admitted to the graduate program in mechanical engineering at the University of Florida in August, 2003. In December 2005, he completed his research work fo r a Master of Science degree with Dr. D. Yogi Goswami on Theoretical Analys is of Solar Driven Flash Desalination System based on Passive Vacuum Genera tion. Shalabh joined Dr. Jac ob N. Chung's research group in mechanical engineering at the University of Fl orida in January 2006 for this doctoral research work. He completed his research work for th e Doctor of Philosophy degree in Fall 2009. 171 