<%BANNER%>

Effect of Vapor Dynamics on Interfacial Instabilities


PAGE 1

EFFECT OF VAPOR DYNAMICS ON INTERFACIAL INSTABILITIES By OZGUR OZEN 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 2004

PAGE 2

Copyright 2004 by Ozgur Ozen

PAGE 3

This doctoral thesis is dedi cated to the teachers all around the world who dedicated their lives to future generations.

PAGE 4

ACKNOWLEDGMENTS First of all, I would like to thank Prof essor Ranga Narayanan for his support and advice. He has been both a mentor and a friend. His enthusiasm for his work and dedication to teaching have inspired me to accept every challenge. Many thanks and love go to my wife, Je n, for her constant support and generous love throughout my graduate education. I also have the highest appreciation for my parents and my sister for their love and nurturing throughout my educational career. Despite the large distance, they have alwa ys been there for me and have always encouraged me to overcome every obstacle. The members of my PhD committee, Dr. Ja son Butler, Prof. Oscar Crisalle and Prof. Loc Vu-Quoc, also deserve my gratitude. I have really enjoyed taking classes from Prof. Vu-Quoc and being a teaching assist ant for Prof. Crisalle. Their teaching philosophies have deeply influenced me. Many thanks go to my friends Mehmet A li Elicin, Eric Theisen, Kerem Uguz and Berk Usta for always being there wh en I needed help and friendship. Special thanks go to Dr. Charles Martucci whose advice change d the direction of my life in order to pursue a career in academics. I cannot repay him for his mentoring and support. My final acknowledgments are to my pare nts-in-law, Betty and Dave Hackworth, who have looked out for Jen and me in every pos sible way. Without their help this thesis would have taken longer to write.

PAGE 5

v TABLE OF CONTENTS page ACKNOWLEDGMENTS...................................................................................................4 LIST OF TABLES...........................................................................................................viii LIST OF FIGURES...........................................................................................................ix ABSTRACT....................................................................................................................... xi CHAPTER 1 INTRODUCTION........................................................................................................1 Physics of Various Instabilities Conn ected to Vapor-Liquid Phase Change................1 Evaporative Instability...........................................................................................1 Rayleigh Convection or Buoya ncy-Driven Convection........................................6 Pure Marangoni Convection or Inte rfacial Tension Gradient-Driven Convection.........................................................................................................8 Rayleigh-Taylor Instability.................................................................................12 Plan of the Thesis........................................................................................................13 2 THE MATHEMATICAL MODEL............................................................................15 Nonlinear Equations...................................................................................................16 Dimensionless Nonlinear Equations...........................................................................19 Some Observations.....................................................................................................21 The Linear Model.......................................................................................................22 Perturbation Equations........................................................................................23 The Base State Solution and the Perturbed Equations........................................24 Endnote: Challenges in phase-change problems: Interface conditions......................28 3 INSTABILITY DUE TO EVAPORATION – THE LINEAR MODEL....................30 The Pure Phase-Change Problem...............................................................................30 Stabilizing Effect of the Vapor Flow..........................................................................38 Endnotes.....................................................................................................................42 The Exchange of Instability when Heated from the Vapor Side.........................42 The Assumption of Flat Surface in the Stability Analysis of the Phase-Change Problem....................................................................................42

PAGE 6

vi 4 WEAKLY NONLINEAR ANALYSIS......................................................................44 The Pure Phase-Change Problem...............................................................................44 Perturbation Expansions.............................................................................................47 Overview.............................................................................................................48 Perturbed Equations.............................................................................................50 Results of Calculations...............................................................................................57 Endnote: The Determination of .............................................................................62 A different approach to the weak nonlinear analysis..................................................64 5 THE RAYLEIGH-TAYLOR PROB LEM WITH PHASE-CHANGE.......................73 Introduction.................................................................................................................73 Physical and Mathematical Model..............................................................................74 Analytical Results for Infinitely Deep Fluids.............................................................78 Numerical Results for Finite Depth Fluids.................................................................84 6 EFFECT OF VAPOR DYNAMICS ON THE RAYLEIGH-MARANGONI PROBLEM.................................................................................................................89 Introduction.................................................................................................................89 Experimental Apparatus and Procedure.....................................................................91 Experimental Results..................................................................................................95 7 CONCLUDING REMARKS AND FUTURE SCOPE............................................102 Active Vapor Layers.................................................................................................102 Conditions at the Interface of an Already Evaporating Liquid.................................104 Coupling of Evaporative In stabilities with Other Phase-Change Phenomenon.......104 APPENDIX A THERMODYNAMIC EQUILIBRIUM RELATION AT THE INTERFACE........105 B A SIMPLE DERIVATION TO OBTAIN THE CRITICAL CONDITION FOR ONSET OF CONVECTION IN A FLUID LAYER................................................108 C A SIMPLE DERIVATION FOR THE RA YLEIGH TAYLOR INSTABILITY....111 D THE BOUSSINESQ APPROXIMATION...............................................................113 E THE UNIT NORMAL OF A SURFACE.................................................................116 F DERIVATION OF THE SURFACE SPEED...........................................................117 G THE INTERFACIAL ENERGY BALANCE..........................................................119

PAGE 7

vii H THE EXPANSION OF A DOMAIN VARIABLE AND ITS DERIVATIVES ALONG THE MAPPING.........................................................................................121 I NUMERICAL METHOD........................................................................................126 LIST OF REFERENCES.................................................................................................129 BIOGRAPHICAL SKETCH...........................................................................................131

PAGE 8

viii LIST OF TABLES Table page 3-1 Physical properties of water and water vapor at 100 C...........................................31 3-2 Definitions of R1, R2 and R3.....................................................................................31 3-3 Comparison of critical temperature differences.......................................................32 3-4 Definitions of R4, R5 and R6.....................................................................................36 6-1 Physical properties of silicone oil and air at 20 C...................................................96

PAGE 9

ix LIST OF FIGURES Figure page 1-1 Physical model of th e phase-change problem............................................................2 1-2a Heat input from the liquid side...................................................................................3 1-2b Heat input from the vapor side...................................................................................4 1-3 Physics with fluid dynamics.......................................................................................5 1-4 Rayleigh convection...................................................................................................6 1-5 Marangoni convection................................................................................................8 3-1 1 R vs. wavenumber..................................................................................................33 3-2 2R and 3R vs. wavenumber....................................................................................36 3-3 4R through 6R vs. wavenumber.............................................................................38 3-4 The velocity vectors in both phases at the onset of instability. z0 represents the interface position in the base state............................................................................40 3-5 The critical temperature difference ve rsus the dimensionless wavenumber for water-water vapor system where d3 mm and *d1.5 mm ...................................41 3-6 Evaporation number vs. dimensionle ss wavenumber for two water-water vapor systems.....................................................................................................................43 4-1 Heat transfer across the bottom plate.......................................................................59 4-2 The scaled pressure perturbations in the vapor domain for the case where d1 mm 0.1 and 1 z0 represents the interface position in the base state..........60 4-3 The scaled pressure perturbations in the vapor domain for the case where d1 mm 0.3 and 1 z0 represents the interface position in the base state..........61 4-4 The scaled pressure perturbations in the more viscous vapor for the case where d1 mm 0.3 and 1 z0 represents the interf ace position in the base state.......................................................................................................................... .62

PAGE 10

x 5-1 The physical model..................................................................................................75 5-2 The dependence of the phenomenologi cal coefficient on the wavenumber............83 5-3 Comparison of models ............................................................................................85 5-4 Effect of shallow fluid layers...................................................................................87 5-5 Critical temperature difference vs. dimensionless wavenumber..............................88 6-1 Different types of convection mechanisms..............................................................89 6-2 Schematic of the experimental setup........................................................................91 6-3 Experimental setup...................................................................................................93 6-4 z-component of velocity in a silicone oil-air system at the onset of instability where the silicone oil and air depths are 5mm and 3mm, respectively....................97 6-5 z-component of velocity in a silicone oil-air system at the onset of instability where the silicone oil and air depths are 5mm and 5mm, respectively....................97 6-6 Flow profiles in a silicone oil-air sy stem at the onset of convection where the silicone oil and air depths are 5mm and 3mm, respectively....................................98 6-7 z-component of velocity in a silicone oil-air system at the onset of instability where the silicone oil and air depths are 5mm and 9mm, respectively....................99 6-8 The critical temperature difference vs. air height for a bilayer system where the silicone oil depth is 8mm and the aspect ratio is 2.5..............................................100 6-9 Pictures of experiments taken by IR camera at the onset of instability.................100 B–1 Critical Rayleigh number vs. dimensionless wavenumber....................................110

PAGE 11

xi 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 EFFECT OF VAPOR DYNAMICS ON INTERFACIAL INSTABILITIES By Ozgur Ozen December 2004 Chair: Ranga Narayanan Major Department: Chemical Engineering Instabilities at fluid-fluid interfaces are of interest to researchers on account of their role in many industrial and na tural processes. These instabilities are often accompanied by phase-change between liquids an d their vapors, and in some instances are caused by it. Examples of processes involving such instabil ities are film boiling and drying processes. An everyday occurrence is seen in the pa ttern formation caused by the drying of coffee drops on a kitchen counter. It turns out that the inst abilities are accompanied by fl uid flow in both phases but many past studies of the phenomena have been done by considering only the liquid fluid dynamics while the role of the vapor phase has been either left out of the study or relegated to heat transfer only while its fluid dynamics has been ignored. In this study, evaporative instabilities are analyzed theoretically and experimentally at the interface between a liquid and its vapor where bot h the liquid and the vapor are fluid dynamically active in order to inve stigate and understand the role of vapor dynamics on evaporative and interf acial instabilities. In the theoretical model, the Navier

PAGE 12

xii Stokes and energy equations in the liquid and vapor domains were using linear and nonlinear perturbation methods w ith the assumption that both fluids are incompressible using the Boussinesq approximation. The linear theory tells us the conditions at the onset of instability and the critical point at wh ich the interface becomes unstable. The weak nonlinear method tells us the behavior of th e system slightly beyond the onset of the instability. The linear and weak nonlinear theo ries are of importance for the experimental reasons as they provide us with the info rmation on how and when the instability will occur. In a separate study of non-evaporative syst ems, experiments were performed in a liquid-gas bilayer system in a circular contai ner with the intention of studying instabilities in the absence of evaporation. The depth of the gas layer was changed to observe the effect of the fluid dynamics in the gas laye r on the onset flow pa ttern and conditions. These experiments were then compared with a corresponding theore tical model and it was concluded that the gas dynamics played an ever-increasing role as the gas layer height was increased. Past workers had a ssumed the vapor to be a conductive, fluid dynamically passive medium. However, th e experiments show us that deep enough vapor layers can cause the bilayer system to convect earlier than expected by thermally coupling with the liquid at the interface. In summary, from theoretical studies it wa s found out that active vapor layers play an important role on the evaporative instab ilities even in the absence of natural convection. Moreover, it is discovered th at the fluid dynamically active vapor plays a very strong stabilizing role on the evaporation problem. In addition to that, from the experimental studies in the absence of phase-c hange it was revealed that the vapor layers

PAGE 13

xiii can and do affect the onset flow pattern and conditions. Both theoretical and experimental studies also indicate that th e commonly used infinite ly deep vapor layer assumption is not always true and can lead to errors.

PAGE 14

1 CHAPTER 1 INTRODUCTION Evaporation and condensation involve the transformation of a substance from the liquid to the vapor state or th e other way around by the addition or depletion of energy. It turns out that these processes can lead to an interfacial instability. This instability is typically manifested by the sudden undulati on of an erstwhile flat surface between a liquid and a vapor as a control variable, such as the temperature drop across the bilayer, exceeds a critical value. Ther e is a need to understand the phenomenon of the instability in the light of the application of evaporation to many technologi es and this is the focus of our work. Evaporation arises in a variety of industrial and natural processes, such as in heat pipe technology, the drying of oxide film s, the drying of lakebeds and in dry-eye syndrome. However, evaporation comes accompanied by fluid flow and natural convection due to buoyancy and interfacial tension gradients. In order to have an understanding of the various phenomena associated with evap oration and the first order effects that promote instability during eva poration we advance some “pic ture arguments,” first for the pure phase-change problem and later for connected instabilities such as convection and gravitational instability. Physics of Various Instabilities Conn ected to Vapor-Liquid Phase Change Evaporative Instability Evaporative instability is best understood by considering an evaporating liquid in contact with its own vapor. To get an idea of the physical model the reader’s attention is

PAGE 15

2 drawn to Figure 1-1, which depicts a liquid-vapor bilayer where a temperature difference is imposed across it. The input variables ar e the fluid depths, the temperature difference and the pressure of the vapor at the upper plat e. As a result of fixing these variables the evaporation rate, the fluid be havior and the interface morphology can be determined. It may be noted that evaporation can be obtained in the absence of a temperature difference by manipulating the pressure at the upper plate. However, in this work we shall imagine that a temperature differen ce is always imposed. The interface is originally thought to be flat a nd it is the stability of this state that is in question. The liquid and vapor depths may be fixed by suitably adjusting the liquid feed and vapor removal rates. In fact, gi ven a temperature drop the upper plate pressure may be adjusted in such a way that the evapor ation or condensation rate is zero. We will consider much of the analysis in this thes is for a case where initially there is no mass transfer across the interface as it simplifies the calculations an d is sufficient to reveal the physics behind this problem. We now pro ceed in a stepwise manner to consider evaporative instability and return to Figure 1-1. Figure 1-1. Physical model of the phase-change problem First, imagine that the base flow rate is zero and the fluid dynamic s is left out. By imagining that the fluid mechanics is omitted, we are also pretending that pressure

PAGE 16

3 perturbations do not arise. Th is will obviously change when the perturbed fluid motion is taken into account. Interesting conclusi ons can be drawn by simply looking at the temperature profiles in the base state. The schematics in Figures 12a and 1-2b depict an evaporation problem at a zero evaporati on rate. Figure 1-2a shows the heating arrangement where the heat needed for evapor ation is supplied from the liquid side while Figure 1-2b shows heat input fr om the vapor side. The arrows in the figures show the direction of the front motion if the base stat e were to have some evaporation though this is not central to the physical arguments presented here. Figure 1-2a. Heat input from the liquid side In both figures the dotted wave represents a perturbed interface while the dotted lines represent the temperature pr ofiles in the perturbed state at the interface. Note that the temperature at the perturbe d interface is the same as the unperturbed state because the temperature is assumed to be at equilibriu m with the pressure at the interface, an assumption that is justified at low evaporati on rates. When the fluid dynamics is ignored the pressure is not perturbed and consequently also the temperature is unperturbed at the interface (see Appendix A for a description of the thermodynamic equilibrium relation at the interface and its development upon the impos ition of an infinitesimal perturbation). When heated from the liquid side, a crest is further away from the heat source and also closer to the heat sink. Therefore the temper ature gradient on the vapor side gets sharper

PAGE 17

4 whereas the temperature gradient on the liquid side becomes less sharp (i.e., less heat is transferred to that point and more heat is ta ken away). The rate of evaporation from the crest is thus decreased which makes the in terface unstable. On the other hand, when heated from the vapor side, the trough gets cl oser to the heat sink and is further away from the heat source. Consequently, less heat is transferred to the fr ont and more heat is taken away from it. The speed of the trough is thus decreased and as a result the front is stabilized. Note that the picture argument presented here suggests that when the interface is unstable, it is unstable for disturbances of all periodici ties or wavenumbers and when it is stable, it is stable for all wavenumbers. Figure 1-2b. Heat input from the vapor side Now include the fluid dynamics in both phases in the thought experiment, but do not consider any other convectiv e effect. Figure 1-3 is similar to Figure 1-2.a except that fluid dynamics is now taken into account in both phases. Observe that a deflected surface does not have the same temperature as an erst while flat interface for a trough is expected to be at a higher temperature than a crest as it is closer to th e heat source. Remember also that we started out with a base state wher e there was thermodynamic equilibrium at the interface and no phase-change was taking place. Now that the trough is at a temperature higher than the saturation temperature ther e will be evaporation taking place there and consequently, when perturbed, the vapor will condense into its own liquid at the crest

PAGE 18

5 where the temperature is colder than the saturation temperat ure. There will be upward flow at the troughs and downward flow at the cr ests as indicated in Figure 1-3. Again, instability occurs at every wavenumber. Figure 1-3. Physics with fluid dynamics Now surface tension acts on pr essure differences in fl uid dynamics problems via the surface curvature and in many phase-cha nge instability problems it is the surface tension that comes in to stabi lize the situation. Its effect is felt most strongly at small wavelengths and weakly on lazy perturbati ons. But, surface tension comes attended by the complication that it depends upon temperatur e and we will take this into account as we continue to study the problem in stages. There is also another complication and this time it is due to gravity If the vapor were at the botto m as in pool boiling then gravity would play the destabilizing role of a h eavy fluid on top of a light fluid. The destabilization would take place in that case when gravitational potential energy overwhelms surface potential energy. This is the familiar Rayleigh-Taylor instability problem and we give a brief discussion of this instability later in this chapter. However, even if the vapor were on the top and the li quid were heated, gravity still plays a role via buoyancy and to understand the nature of natural convect ion due to buoyancy and interfacial tension gradient we tu rn to more picture arguments.

PAGE 19

6 Before we move on to a discussion of thes e other effects, it might be noted that evaporative instability can take place even if the interface we re assumed to be flat. In this case, the instability would be manifested by fluid motion. This case will not be of interest to us. Rayleigh Convection or Buoyancy-Driven Convection Buoyancy-driven convection, often term ed natural convection or Rayleigh convection, occurs when a fluid is subject to a temperature gradient in a gravitational field and when there is a variation of density with respect to temperat ure. Imagine a layer of liquid bounded vertically by two horizon tal surfaces, with the lower surface at a temperature greater than the uppe r one as in Figure 1-4. Figure 1-4. Rayleigh convection As density typically decreases with an in crease in temperature, the fluid near the top surface is heavier than the fluid at the bot tom plate, creating a t op-heavy arrangement, one that is gravitationally uns table. For sufficiently small temperature differences, the fluid simply conducts heat from the lower surface to the upper one, creating a linear temperature field across the fluid. Now imagine a mechanical perturbation being imparted to an element of fluid so that it is projected downward. As the density of the “fluid packet” is greater than its environm ent it will proceed downward. Due to mass conservation, fluid from below will be disp laced and will move upward. This motion

PAGE 20

7 will continue unabated unless the fluid’s kinematic viscosity and thermal diffusivity are large enough that the mechanic al perturbations and thermal perturbations die out quickly and the original mechanically quiescent st ate is re-established. From this thought experiment it is immediately clear that unlike many classical problems of fluid mechanics, the problem of the onset of convec tion can not be considered for fluids in the limit of vanishing viscosity or thermal diffusivi ty. In other words, disturbances to the fluid are dampened by the kinematic viscosity and thermal diffusivity of the fluid and enhanced by the buoyancy forces. Under stea dy conditions there is a balance that is established between these effects, represen ted by the dimensionless Rayleigh number, a group that arises from the scaling of th e modeling equations. It is given by 3T Rad g 1.1 Here, is the negative thermal expansion coe fficient and is usually positive, g is the magnitude of gravity, T is the vertical temperature difference across the fluid layer and is taken positive when the fluid is heated from below, d is the depth of the fluid, is its kinematic viscosity, and is its thermal diffusivity. If the temperature difference is increased beyond what will be called the critical temperature difference, then the gravitational instability overcomes the viscous and thermal damping effects and the fluid is set into motion, causing buoyancy-driven convection. (See Appendix B for an elementary derivation of the critical condition for the onset of buoyancy driven convection.) If the upper surface adjoins a va por phase then convection can initiate in either phase. If it begins in the lower phase it will propagate to the upper phase by way of viscous dragging. If convection is si multaneous in both phases they can become

PAGE 21

8 thermally coupled with the in terface becoming an isotherm. If it initiates in the upper phase, it cannot propagate to the lower phase except by way of generating interfacial tension gradients, which in turn will generate weak flow in the lower liquid. For that reason, we turn to interfacial tension gr adient-driven or Mara ngoni convection in the following section. Pure Marangoni Convection or Interfacia l Tension Gradient-Driven Convection Interfacial tension gradient-driven convection or Marangoni convection is unlike buoyancy-driven convection and can occur in a fluid even in the absence of a gravitational field. Imagine a layer of fl uid which is bounded below by a rigid plate and whose upper surface is in contact with a passi ve gas as depicted in Figure 1-5. Figure 1-5. Marangoni convection A passive gas is not one without density or viscosity. It is merely a gas wherein one tacitly assumes that therma l or mechanical perturbations cannot occur. Pretend that above the passive gas, there is another rigid plate. For the sake of consistency with the earlier argument, let the lower plate be at a temperature hi gher than the upper plate's temperature. Now, imagine that the inte rface between the lower liquid and the passive gas is momentarily disturbed. The regions of the interface that are pushed up experience a cooler temperature. Likewise, the re gions of the interface that are pushed down x y

PAGE 22

9 increase in temperature. Typically, interf acial tension decreases with an increase in temperature. Therefore, the regions of th e interface, which are pushed up increase in interfacial tension, which pulls on the interface, while the re gions of the interface, which are pushed down, decrease in surface tens ion. When the fluid is pulled along the interface, warmer fluid from the bulk replaces the fluid at the interface near `x’ enhancing the interfacial tension gradie nt-induced flow. If the temp erature difference across the liquid is sufficiently small, then the thermal diffusivity of the fluid will conduct away the heat faster than the disturbance can amplify. The dynamic viscosity will also resist the flow causing the interface to become flat ag ain and the interfacial tension to become constant. As in the case of buoyancy-driv en convection, there exists a critical temperature difference when the surface tensi on gradient-driven flow is not dampened by the thermal diffusivity or viscosity, and the fluid is set into motion. Scaling of the equations that model the interfacial moment um balance leads to a dimensionless group that exhibits the balance between dissipat ive effects of the viscosity and thermal diffusivity and the promoting effects of interf acial tension gradient. This dimensionless group is called the Marangoni number and is given by TT Mad 1.2 Here T is the change in the interfacial tensi on with respect to the temperature, and is the dynamic viscosity of the fluid. It is noteworthy that while the above argument assumes the interface to be deflectable it is not necessary for the interface to undulate in order for Marangoni convection to take place. To see this a little more clearly, look at Figure 1-5, but pretend that the interface is flat. Points ‘x’ and ‘y’ will therefore be on

PAGE 23

10 the same horizontal level. Imagine an el ement of fluid projected downward from the interface at point ‘y’. Fluid n earby must replace the lost flui d at ‘y’ and this fluid comes from the region near ‘x’. Consequently, hot fluid from below moves to ‘x’. If the interfacial tension decr eases with temperature, as is th e case in most substances, then fluid flows from warm surface regions to cold surface regions. Therefore, fluid moves from ‘x’ to ‘y’ and the convec tion continues unless large thermal diffusivity and viscosity dampen out the disturbances. This mean s, theoretically speaking, that Marangoni convection can be analyzed assuming a flat fr ee surface, i.e., a surface with an infinite surface tension. This fact was observed by Pearson (1958) who gave the first theoretical analysis of Marangoni convection using a me thodology similar to that used by Rayleigh for the case of buoyancy driven convection. The physical argument presented above can easily be extended to a bilayer system where the gas layer is fluid dynamically active. In this case, the gas layer will also convect due to surface tension gr adients and will bring in cold er fluid from the cold plate to ‘x’ opposing the destabilizing nature of the convection in the liquid; however, the effect of the convection in the gas is le ss because of its high kinematic viscosity. The extent of either Rayleigh convecti on or Marangoni convec tion is primarily a function of the fluid depths. By examini ng equations 1.1 and 1.2, we notice that Rayleigh convection is proporti onal to the cube of the fl uid depth and that Marangoni convection is directly proportional to the flui d depth. From these scaling arguments, we can conclude that for deep fluids, buoyancy-dr iven convection is more prevalent, whereas for shallow depths, interfacial tension gradie nt-driven convection is more prevalent, and at intermediate fluid depths, both Rayleigh and Marangoni convection can occur. If the

PAGE 24

11 interface is allowed to deflect one might an ticipate that in the case of pure Marangoni convection, warm fluid flows towards a trough while in the case of Rayleigh convection cold fluid flows away from a trough leadi ng us to believe that Marangoni motion can counter Rayleigh driven motion. This counter effect was demonstrated computationally by Sarma (1987) and could not be seen in th e study by Nield (1964) wherein flat surfaces were assumed. When a layer of fluid is heated from above it creates stable density stratification. Therefore the buoyancy force does not cause co nvection. Marangoni convection, though, may still occur in fluids being heated from a bove. For example, consider a single layer of fluid superimposed by a passive gas and heat ed from above. One might imagine that purely steady Marangoni convection would not occur. For instance, suppose a random perturbation causes some part of the surface to become warmer than the rest of the surface. The interfacial tens ion will decrease in this re gion and the higher interfacial tension elsewhere will pull fluid away from th is hot spot. By continuity, fluid lying below the hot spot will move up to replace the displaced fluid. As the lower fluid is cooler than the surface, this region now cool s off and the interfacial tension increases, thereby either stabiliz ing the region or causing an overshoot in temperature of that spot and generating oscillatory motion. Only cal culations can tell precisely which case will occur. In real systems the upper fluid is ne ver truly passive. Given the same scenario, fluid movement along the interface will also dr ag warmer fluid from above and cooler fluid from below. Depending upon the fluid pr operties and depth ratios the instability will or will not be reinforced.

PAGE 25

12 Oscillatory convection can occur if the fluid is heated from above and buoyancy effects are included. This wa s shown by Rednikov et al. (1998) for the case when the upper fluid is passive. The explanation is as follows. If a small volume of fluid is displaced within the bulk of the fluid, the dens ity stratification acts as a restoring force. The Marangoni force acts similarly and the two effects reinforce each other causing an overshoot and sustained oscillati ons. Indeed, as was demonstr ated in their paper, this only occurs in certain fluids, at certain dept hs, where the buoyancy a nd interfacial tension forces are comparable. There is another phenomenon associated with surface tension gradient-driven convection, often called the long wavelength Marangoni instability. This instability typically occurs when either the surface tension or the depth of the fluid layer is very small. The initiation of this instability is similar to the description given above. However, in the long wavelength scenario, the convection cells are much larger than the regular, or short wavelength, Marangoni conve ction. As the convection propagates, it causes large-scale deformation in the interface, which can actually cause the interface to rupture; that is, the interface deforms to such an extent that it comes in contact with the lower plate. This phenomenon occurs in th e drying of films and in coating processes (Schatz). Rayleigh-Taylor Instability So far we have talked about the liqui d and gas systems where the liquid in consideration underlies the gas phase. There ar e however situations, such as film boiling, where the opposite arrangement is encountered When the vapor pha se underlies a liquid phase, this is a top-heavy arrangement cal led Rayleigh-Taylor instability and it is potentially unstable. The instability occu rs when the gravitational potential energy

PAGE 26

13 exceeds the stabilizing e ffect of surface potential energy. In the classical treatment of the problem, it is assumed that there is no heat and mass transfer, and it is well known that the system is unstable to disturbances w ith a wavenumber smaller than the critical wavenumber given by 1 2 cg where is the surface tension, and are the liquid and vapor de nsities, respectively. A derivation for the classical Rayleigh-T aylor critical wavenumber from the book by Chandrasekhar (1960) for infinitely deep flui d layers is given in Appendix D. This instability is also discussed in connection to the phase-change problem in Chapter 5. What makes the Rayleigh-Taylor problem of sp ecial interest to us is the unavoidable existence of both the heat and mass transfer in film boiling. Plan of the Thesis The rest of the thesis cont ains the modeling equations, their scaling and the method of solution to determine the onset conditions of the instabilities. This is followed by a discussion of the results of calculations for the onset of instability in the evaporation problem, the Rayleigh-Taylor problem and the pure convection problem. In each of these discussions, the role of vapor dynamics is st ressed: a role that has not been emphasized up to this point in this chapter. The vapor is a source of heat tran sfer in all of these problems. To assume that the vapor plays a passive role is then ordinarily incorrect. When the vapor moves it does so on account of several mechanisms. It may move because pressure perturbations are imparted to it or because the thermal profiles generate phase-change at least upon perturbation, or it can move on account of viscous dragging.

PAGE 27

14 Whatever the mechanism for its motion, it must change the mode of heat transport from conduction to convection. This change in the transport must feed b ack to the interface and either assist the undulations or hinder them. Our goal is to understand the role of the feedback, the modeling and the physics of instability mechanisms. This thesis is different from the traditi onal type of dissertations because of the absence of a chapter on the litera ture search. I have decided to refer to the literature in the context of the work presented.

PAGE 28

15 CHAPTER 2 THE MATHEMATICAL MODEL This chapter includes all of the equations used to analyze the problems described in Chapter 1. The equations are divided into th e field equations that are meant to describe the behavior of the state doma in variables and the interfacial equations that describe the relation of the domain variables with surface variables such as surface speed and surface curvature which, in turn, depend upon surf ace position. The domain equations and interfacial conditions are written in a gene ral manner without rega rd for distinction between the various problems. We shall draw the attention of the reader to the differences as we go along and thereby indi cate how each problem belongs to a special case of the general equations. What may not have been apparent from th e discussion in the earlier chapter is that the instabilities are related to a nonlinear property of the modeling equations. The nonlinearities arise on acc ount of more than one reason. One of the reasons to expect nonlinear characteristics in the modeling e quations is that the interface position is coupled to the transport behavior and the tw o depend on each other. Another reason is that the domain variables are coupled. For example, the interface temperature distribution affects the pressure distribution, which in turn a ffects the flow profiles. The flow profiles then affect the temperature prof iles, either reinforcing them or causing the disturbances to subside. We now turn to th e nonlinear equations a nd shall point out their characteristics as we go along.

PAGE 29

16 Nonlinear Equations The various problems we will deal with ar e all fluid flow problems; hence, the equations that model the fluid physics in uns caled form are given by the Navier Stokes equations in each phase, i.e., by 2v vvPT g v t 2.1 and *******2*v vvPT g v t 2.2 Each of the problems is connected with temperature gradients. These affect the energy transport. Therefore the energy equa tions in each phase are needed. They are given by T T v t T2 2.3 and 2 * *T T v t T 2.4 Mass conservation in each phase requires the continuity equations. They are v0 2.5 and *v0 2.6 In the equations above, v P and T are velocity, pressu re and temperature fields, respectively, and the aste risk denotes the vapor ph ase. In the equations, , and are the density, viscosity and ther mal diffusivity, respectively. A major assumption that has tacitly been made is that the densities of the fluids are treated to be constant in all of the term s save those associated with gravitational

PAGE 30

17 acceleration. This assumption is fair provided the temperature drop is not very large or at least as long as the total density change divi ded by its mean value in each phase is a small fraction. This is the celebrated Boussinesq approximation. It can be justified by a perturbation series in a scaled thermal expansion coefficient. This justification, given by Mihaljan (1962), has proven to be invaluable in studies in ther mal convection greatly simplifying the calculations. More on this approximation is given in Appendix E. We continue with the modeling equations and assume that the bottom and top plate temperatures are kept constant and therefore b ottomTdT and ** topTdT hold. The top and bottom walls are impermeable, therefore, z-dvn0 and ** zdvn0 hold. The no-slip condition applie s along the top and bottom plates, and gives rise to z-dvt0 and ** zdvt0 At the interface, the mass balance equation is given by **vunvun 2.7 where 2 1 2x Z 1 x Z n k i is the unit outward normal (Appendix F) and 2 1 2x Z 1 t Z n u is the interface speed (Appendix G). When there is no phasechange at the interface, we se t the interfacial mass balance equal to zero, and we obtain two interface conditions from the in terfacial mass balance equation: vun0

PAGE 31

18 and **vun0 At the interface, the temperatures and the tangential components of velocities of both fluids are equal to each other, so *T T and t v t v* hold where 2 1 2x Z 1 x Z t k i is the unit tangent vector. The interfacial tension and its gradient with respect to the temperature at the interface come into the picture through the fo rce balance. By ta king the dot product of the stress balance with the unit normal and unit tangent vectors separately, we get the normal and tangential stress balance equa tions. The force balance is given by **** s v vuTn2H n t v vuTn 2.8 where S I -P T is the total stress, is the interfacial tension and 2 3 2 2 2x Z 1 x Z H 2 is the surface mean curvature. In addition we need the inte rfacial energy balance, viz. 2 2 *******1 ˆ vuvuqSvun 2 1 ˆ vuvuqSvun 2 H H 2.9 Here H ˆ is the enthalpy per unit mass, S is the extra stress, and q is the heat flux. Now, if there is phase-change across the interf ace, we are short of one equation as the

PAGE 32

19 interfacial mass balance (Equation 2.7) gives us only one relation. To get this additional equation we observe that when the evaporation rate is very small we may assume that there is thermodynamic equilibrium at the phase-change boundary, that is *** base b aseT PPln T 2.10 where b aseP and b aseT are the base state pressure a nd temperature of the fluid, and is the latent heat of evaporation per unit mass. Dimensionless Nonlinear Equations The nonlinear equations in dimensionless fo rm, if necessary, are given below. The following scales are used: length, d; velocity, d v ; time, v d; pressure, v P d and temperature, coldTT T T where cold hotT T T The instability in these problems is caused by the temperature gradients, so it is reas onable to scale the ve locities with respect to the thermal diffusivity since it is the dissipa tion of the temperature perturbations that is important. The temperature is scaled with respect to the total temperature difference across the bilayer system and the reference temp erature is the cold temperature, which in many systems we cover will be the top plate. This will prove to be of convenience when we move to higher order expansions of the problem. Throughout the course of our study we will keep the top and bottom plate temperatures constant and by scaling the temperature with respect to our choice of re ference temperature and the total temperature difference, we will significantly simplif y the top and bottom temperature boundary conditions. Hereafter, we will drop the t ilde on the scaled temperature. Using the expansion for the density given in Appendix E, the momentum equations in each phase become

PAGE 33

20 21v +vvP Ra T v PrtzzGii 2.11 and *****2* ***1v +vvP Ra T v PrtzzGii 2.12 where Pr and 3gdG are the Prandtl and Gravity numbers, respectively. The energy equations in each domain are given by 2T T+vT t 2.13 and 2*** *T T+vT t 2.14 The mass balance at the interface, (,) zZxt now becomes *vunvun 2.15 The normal stress balance in dimensionless form is given by ** ***11 Cav vuT:nn2HCav vuT:nn PrPr 2.16 where d Ca is the Capillary number. The tangential stress balance becomes sS:ntS:ntMaT0 2.17 In addition, the interfaci al energy balance becomes

PAGE 34

21 2 2 ** PC * PC11k 1KvuvuvunETnTn 22k VSvuSvun0 2.18 Here, 2 PC 2K d, pC T E and PC 2V d E is the Evaporation number. At the interface, the temperatures of both fluids are equal to each other, and the thermodynamic equilibrium at the phase-change boundary holds, that is * KEBASE b aseT PPln T 2.19 where KE *2 d Some Observations We can make several observations on th e above equations. First, as noted, the interfacial mass balance is simplified when there is no phase-change and the thermodynamic equilibrium relationship is then ignored. It is also possible to consider phase-change with a zero temperature drop acro ss the bilayer. However, in this case a non-zero evaporation base state is required. This case is not considered in this study. Second, when phase-change takes place and the thermodynamic equilibrium is taken into account, the effect of curvatur e on latent heat also needs to be taken into account. However, the effect is weak and this fact is proven in Appendix A. Third, by imposing a restriction of a flat undulated interface and infinite interf acial tension, the equations simplify considerably. The normal component of the stress balance is then used only to determine the value of the curvature times th e interfacial tension but not to close the problem. Fourth, if the Ma rangoni effect were excluded, then the ta ngential stress

PAGE 35

22 balance at the interface simplifies by rendering the interfacial tension gradient to zero but no singularity arises on account of this. Fina lly, the interfacial ba lances are all frame invariant. This is seen in the force balan ce if the mass balance is taken into account and it is seen in the energy balan ce in the form that it is pres ented. The energy balance is derived for the reader’s benefit in Appendix H. The Linear Model To see the equations that determine the onset of instability requi res us to linearize the equations about a base state. First, a wo rd about the base state: The base state is always a solution to the nonlinea r equations and often it might seem defeating to look for a base state if it means solvi ng a full set of intricate coupled nonlinear equations. But in practice for a large class of problems the base state is seen almost by inspection or by guessing it. For example, in the problems consid ered in this thesis an obvious base state is the quiescent state with a flat interface. It is not in global thermodynamic equilibrium as a temperature gradient is imposed but it is a state whose stability is called into question. To study the stability of the interface, ar bitrary disturbances of infinitesimal amplitude are applied and for a given set of input variables, the growth and decay time constants of these disturbances or equivalen tly the critical temper ature difference across both layers that result in ma rginal stability are determined. The reader might wonder why the disturbances are taken to be small and what this restriction implies. If a state is unstable to infinitesimal disturbances it must be unstable to all disturbances. The same is untrue of finite amplitude disturbances. T hus, stability to small disturbances usually means very little but instability implies a lot. Using small disturbances also allows the dropping of terms with quadratic interactions and higher degree interactions. Therefore,

PAGE 36

23 the plan is to linearize the above equations about a known base state and to find the onset of the interfacial instability from the perturbed model. Perturbation Equations Hereon, the variables of the base state are denoted by the subscript 0, and the variables of the perturbed st ate by the subscript 1. Thus, the temperature, when perturbed, is described by ) ( z dz dT T T T2 1 0 1 0 O where is a small perturbation parameter repres enting the deviation from the base state and 1z is the mapping of the pertur bed configuration onto the refe rence configuration. Its meaning is explained in the Appendix I and, at the interface, it is simply the perturbation of the surface deflection to first order, Z1, a variable to be determined during the course of the calculation. We can further expand 1T and other subscript ‘one’ variables using a normal mode expansion. For example, x i te e z 1 1 ) ( T ˆ T Here, is the inverse time constant while is a wavenumber associated with the given perturbation. A wavenumber arises because th e system is infinite in lateral extent. However, the wavelength, which is proportional to reciprocal of th e wavenumber, may be related to the spacing between the sidewalls in cont ainers with periodic boundary conditions. The same expansion as above is used for both components of velocity and pressure as well as temperature in both phases. It al so holds for the interface deflection variable

PAGE 37

24 Z1. The first priority is to get the base state, as this is the state whose stability is in question and about which the expansions must be done. The Base State Solution and the Perturbed Equations In the base state there is no flow in either the liquid or the gas. For the phasechange problem, assume that there is no phase -change. In other words, assume that the upper plate pressure is fixed so that the phase-c hange rate in the base state is zero. This of course does not mean that the phase-c hange rate will be zero upon perturbation. Consequently, both fluids are in a conductive state, thus 0 v v* 0 0 Hereon, we will drop the overbar of the scaled temperature. The temperature profiles in the liquid and its vapor in dimensionless form become ** T kk z kkkk and ** T kk z kkkk where d d*. After we perturb the domain and boundary equations in the manner given in the previous section, we arrive at the following equations in the domain in dimensionless form: 2 2 x11x1 2d1 v P v dzPr wiw 2.20 2 2 1 z11z1 2P d1 vRa T v dzzPrw 2.21

PAGE 38

25 1 z1 0 1 2 2 2T v dz dT T dz d w 2.22 and 0 dz dv vz1 x1 iw 2.23 for the lower phase, and 2 2*** x11x1 2**d1 v P v dzPrwiw 2.24 2* 2*** 1 z11z1 2***P d1 vRa T v dzzPrw 2.25 1 * z1 0 * 1 2 2 2T v dz dT T dz d w 2.26 and 0 dz dv v* z1 x1 iw 2.27 for the upper phase. At the lower plate, 1 z, the liquid is subject to 0 vz1, 0 vx1 and 0 T1 while the conditions at the top plate become 0 v* z1, 0 v* x1 and 0 T* 1. At the interface the mass balance turns into 1 * z1 z1 Z 1 v v 2.28 which, in the absence of phase-change, becomes z11v Z and z11v Z. The continuity of temperature a nd the no-slip conditions become

PAGE 39

26 1 0 1 1 0 1Z dz dT T Z dz dT T 2.29 and x1 x1v v. 2.30 The normal and tangential stress balances assume the following dimensionless form: * *2 z1z1 111dvdv Ca PP2 Ca Bo Z0 dzdzw 2.31 and 0 Z dz dT T Ma v dz dv v dz dv1 0 1 z1 x1 z1 x1 iw iw iw 2.32 where 2 *d g Bo and again d Ca and d T MaT are the Bond, Capillary and Marangoni numbers, respectivel y. The interfacial energy balance is 1 1 1 z1 Z dz dT dz dT E v k k 2.33 where E stands for the Evapor ation number and is given by pCT E pC is the liquid heat capacity. Again, one has to be careful in order to avoid a mistake here. The interfacial energy balance in the abse nce of phase-change simply becomes * 11dTdT dzdzk k The equilibrium condition at the interface, when there is phase-change across the interface, is * 0 11 KE1PE1 basedT TZ dz P Z T 2.34

PAGE 40

27 where PEg d is a dimensionless parameter from the linearized Clapeyron equation, and b aseT is the scaled base state temperature of the vapor at the flat interface and it is equal to 0T. There are several observations to be made in the perturbed equations. First, the growth constant appears in the domain as well as interfacial equations. Dropping them in the domain equations lead to an inte rfacial mode, however dropping them may not always be justifiable but when the growth constant is dropped from the domain its calculation becomes straightforward. It must th en also be real as it is given by the ratio of two determinants. Second, the growth c onstant is ordinarily an eigenvalue for a sensible control parameter such as E. Third, it is conceivable that solutions exist such that Z1 is zero but they may not be the first to become unstable as a control parameter, such as E, changes. Fourth, when we pay cl oser attention to some of the dimensionless numbers we recognize that their values are very small. A typical value for Capillary number is of the order of 810 and for KE is of the order of 1010 This indicates that for systems where there is a high surface tens ion between the two fl uids it may not be a bad idea to set Ca equal to zero, which conse quently will imply that the surface deflection 1Z is equal to zero. That is to say that we can assume the interface to be flat at all times and the instability will still set in because the temperature and velocity will still affect each other through the domain equations. The same goes also for the KE term. In the next chapter the pert urbation problem is solved for the special case of pure evaporation in the absence of natural convection, and from the calculations a physical interpretation of the problem is sought.

PAGE 41

28 Endnote: Challenges in phase-change problems: Interface conditions There is a crucial difference between the classical Rayleigh-Marangoni-Bnard problem and the evaporation problem. The di fference lies in the inte rfacial conditions. In the Rayleigh-Marangoni-Bnard, the inte rfacial mass balance gives us 2 interface conditions as a result of the no mass transfer condition at the interface. However, in the evaporation problem, we only get one conditi on from the interfacial mass balance, and consequently, we are 1 interface condition shor t. The interface conditions are completed with our choice of equilibrium conditions for the temperature and pressure at the interface. The choice of conditions at th e interface between the liquid and the vapor poses a challenge in phase-change problems. Or dinarily, the continuity of temperature is used along with some other thermodynamic equilibrium condition (Huang and Joseph, 1992) or a relation derived from the kinetic th eory (Oron et al., 1997) is used instead. Ward and Stanga (2001) have investigated the assumption of continuity of temperature and they have verified the existence of a temperature discontinuity at the interface. Margerit et al. (2003) studied the role of such interfaci al nonequilibrium effects on Bnard-Marangoni instability. However, Ward and Stanga report that the magnitude of the temperature discontinuity increases with the evaporation flux and note that when equilibrium exists between the liquid and vapor phases, the temperature is the same in each phase. Consequently, when the evaporati on rate is very small we may assume that there is thermodynamic equilibrium at the ph ase-change boundary and the temperature of both fluids are equal to each ot her at the interface, that is *T Pf and *T T.

PAGE 42

29 It is no problem to consider non-zero evaporation in the base state. The equilibrium relation is replaced by one that involves a phenomenol ogical coefficient, which may, in fact, depend upon the wavelength of the disturbance.

PAGE 43

30 CHAPTER 3 INSTABILITY DUE TO EVAPORATION – THE LINEAR MODEL In this chapter we shall concentrate on the solution of the linearized equations given in the last chapter to find the effect of active va por layers on the evaporative instability without and with convection. The calculations are presented first for the phase-change problem without the comp lications of Rayleigh and Marangoni instabilities. The effects of Rayleigh and Marangoni convect ion are added later and each step towards the complete phase-change problem is discussed. The results for different cases are presented in the form of plot s and physical explanations are given. The Pure Phase-Change Problem To understand the physics of the pure phase-change problem, we first display the results of calculations using the full set of phase-change problem equations as mentioned in Chapter 2. The only terms left out of the equations are the terms that have gravity and surface tension-gradient in them. The first goal of the calculations is to determine the conditions at the onset of instabi lity. In order to do this, we set the dimensionless inverse time constant, to zero. This implie s that we seek only the conditions for steady bifurcation from the base state. Another set of calculations, where we have calculated the values of was also performed for the case of fl uids heated from the liquid side where it was discovered the leading to be always real. Upon setting to zero, the unscaled eigenvalue in the linearized problem is the temperature drop. The Chebyshev spectral tau method (Johnson, 1996) is used to solve the re sulting eigenvalue prob lem that determines the critical temperature difference. The inte rested readers are referred to Appendix J for

PAGE 44

31 more details on the method. The calculati ons are done assuming the physical properties of water and water vapor at the satura tion temperature under 1 atm, i.e., at 100 C as given by Rubinstein et. al. ( 2002). The properties used in the calculations are given in Table 3.1. Table 3-1. Physical properties of water a nd water vapor at 100C 3m kg 960 3 *m kg 6 0 sec m kg 10 9 24 sec m kg 10 3 15 C sec m J 10 8 61 k C sec m J 10 5 22 k sec m 10 7 12 7 sec m 10 0 22 5 C 1 10 0 64 C 1 10 0 63 kg J 10 3 26 m N 10 8 52 C m N 10 0 24 T To make a proper comparison we consider different physical problems and present graphs of the ratio of the critical temperat ure differences for the different cases against the wavenumber. The definitions of these ratios are shown in Table 3.2 where the subscripts pc, Ma, g, and Ra stand for pha se-change, Marangoni convection, gravity and buoyancy-driven convection, respectively. Table 3-2. Definitions of R1, R2 and R3 1R Ma Ma pcT T 2 R Ma pc g Ma pcT T 3 R Ma pc Ra Ma pcT T

PAGE 45

32 For example, p c-MaT denotes the critical temperat ure difference for the case where there is phase-change and Marangoni convection, and 1 R therefore is the ratio of p c-MaT to the critical temperature difference for the case of pure Marangoni convection, MaT First, we present the results for the case where the critical temperature difference for a pure phase-change problem is compared to the critical temperature difference for a case where the Marangoni effect is added to the phase-change problem. This comparison is given in Table 3-3 and the readers can s ee for themselves that the addition of the Marangoni effect does little to change the threshold of evaporative instability. Table 3-3. Comparison of cr itical temperature differences Wavenumber, C Tpc C TMa pc 0.1 0.03 0.03 0.2 0.47 0.47 0.3 1.78 1.78 1.1 6.70 6.69 1.2 6.96 6.95 1.3 7.25 7.24 2.1 10.64 10.63 2.2 11.24 11.23 2.3 11.89 11.88 3.1 19.21 19.19 3.2 20.44 20.42 3.3 21.75 21.73 Marangoni convection is driven by the temp erature variations along the interface. In order to understand the eff ect of phase-change on temperat ure gradients at a deflected interface, let us first consider such an interface in the absen ce of gravity and observe that the deflection must make the temperature vary along the interface. As discussed earlier, the temperature at a trough is expected to be high er than that at a crest as it is closer to the heat source. When perturbed, the liquid starts to evaporate into its own vapor at a trough since the temperature at the trough is now hi gher than the base state temperature.

PAGE 46

33 Likewise, condensation takes place at the cr ests. When the liquid evaporates at the trough, the otherwise hot trough temperature goe s down since part of the energy input is now used for phase-change whereas at the cres t the heat is released, and the otherwise low crest temperature goes up. In other words this phase-change action reduces the temperature perturbations along the interface in phase-change problems, which, in turn, reduces the surface tension gradient-driven convection. Marangoni convection therefore ultimately never plays a significant role on th e stability of the phase-change problem. We can also predict the reduction of the te mperature perturbations at the interface by considering the linearized thermodynamic equilib rium relation at the interface. They are seen from the terms in the equation 2.34 in Chapter 2, which are very small for most fluids, thereby making the temperat ure perturbations very small. 0.3 0.35 0.4 0.45 0.5 0 1 2 3 4 5 6 7 Dimensionless wavenumber, R R1 Figure 3-1. 1 R vs. wavenumber. In Figure 3-1, we observe that at larg e wavenumbers the pure Marangoni problem is considerably more unstable than the phase-change problem with Marangoni

PAGE 47

34 convection. At low wave numbers the reverse is true. Before we go on and explain these observations let us make clea r the restrictions behind the calculations that generated Figure 3-1. In producing Figure 3-1 we have perform ed calculations by taking two similar systems where a liquid underlies its own vapor at zero gravity, and in each system heat is supplied from the liquid side. In one case, we pretend that there is no phase-change across the interface between the liquid and its vapor. In this case, the flow can only be caused by surface-tension gradients along the in terface. In the other case, we allow phase-change across the interface. To make a fair comparison, as there is initially no flow in the pure Marangoni problem, we adju st the pressure on th e vapor side in the second case such that the phase-change rate is zero. As a result, phase-change is only possible through perturbations at the interface. Accordingly, we have two systems with no-flow in their base states. They are two similar systems with only one difference; one has the possibility of phase-change upon pe rturbation while the other has no such possibility at all. We then plot the ratio, 1 R versus wavenumber, in Figure 3.1. At small wavenumbers, not at or very ne ar the zero wavenumber, the pure Marangoni problem is more stable than the phase-ch ange problem with Ma rangoni convection for a water and water vapor system. Small wavenumbers imply that there are reduced transverse thermal gradients thereby stabi lizing the Marangoni eff ect, yet allowing the phase change instability to continue unabate d. At large wavenumbers, we observe the strong action of thermal and momentum diffusion. In many physical situations and laboratory experiments, gravity plays an important role. Gravity comes equipped with the burden of buoyancy-driven convection.

PAGE 48

35 However, the effects of Rayleigh convection can be avoided in appli cations that involve very thin liquid and vapor layers even though gravity continues to play a role. Gravity simply pulls the perturbed interface back to its original position. In other words, gravity plays a stabilizing effect at all wavenumbers. Even so, its effect on the stability of the interface is more obvious at small wavenumbe rs, as small wavenumbers imply weaker surface curvature, and that explains why the critical temperature difference is higher at small wavenumbers when gravity is present. This is depicted by the curve labeled 2R in Figure 3-2. Working with very thin fluid layers is di fficult in a laboratory environment. Thus, eventually we have to acc ount for the effect of buoyancy on the problem. For that reason, Rayleigh convection is now introduced into the problem. It is well known from the Rayleigh problem that, in the absence of phase-change, buoyanc y-driven convection can easily be stabilized at very small and ve ry large wavenumbers. Consequently, in Figure 3-2 we observe that when 2R and 3R are plotted against the wavenumber, both plots understandably merge at small wavenumb ers. At large enough wavenumbers, in the absence of buoyancy-driven convection, the st ability offered by gravity is swept away by the stabilizing effect of the phase-change mechanism and the curve, depicted by 2R, merges with the unity line. This occurs b ecause the surface curvature and, in turn, the phase-change becomes more important compared to gravity at large wavenumbers. In contrast, 3R diverges to become less stable because of the destabilizing character of buoyancy-driven convection. Yet, at even larger wavenumbers 3R ultimately merges with the unity line, no t shown in Figure 3-2.

PAGE 49

36 1 2 3 4 5 6 7 8 9 0 0.5 1 1.5 2 2.5 Dimensionless wavenumber, R R2R3 Figure 3-2. 2R and 3R vs. wavenumber. In order to demonstrate the cr ucial role of the ac tive vapor layer in the stability of a bilayer phase-change problem, we vary the dept h of the vapor layer for a constant liquid depth and calculate the critical temperature difference for a phase-change problem with Marangoni and Rayleigh convecti on. The results are plotte d in Figure 3-3, where the ratio of the critical temperature differen ce for systems with different values of is graphed against wavenumber. The ratios are defined in Table 3-4. Table 3-4. Definitions of R4, R5 and R6 4 R 0.1 with system a of T 0.2 with system a of Tcritical critical 5 R 0.1 with system a of T 0.3 with system a of Tcritical critical 6 R 0.1 with system a of T 0.4 with system a of Tcritical critical Figure 3-3 shows us that increasing the vapor depth stabilizes the interface. The depth of the vapor layer plays a dual role wh en we account for the Rayleigh convection in our problem. The deeper the vapor layer, th e more unstable it becomes due to Rayleigh convection which is scaled as the cube of the depth; however, the important stabilizing

PAGE 50

37 role of the increased vapor de pth on the interfacial stability due to phase-change is still present. The stabilization comes about becau se of the increase in the vapor flow; meanwhile the destabilizing effect of heat tr ansfer is also weakened. To understand the stability mechanism in Figure 3-3, let us firs t look at the problem from a heat transfer point of view. The heat that is transported to the inte rface is used for phase-change and kinetic energy with the balance being conducted away from the interface through the vapor. In our model, the interface temperat ure in the base state is al ways equal to the saturation temperature. Its value depends on the pressure at the top plate and the top plate pressure is set to make the evaporation rate in the ba se state equal to zero. Hence, a change in fluid depths must influence th e pressure required at the top plate to keep the phasechange rate equal to zero a nd thence the saturation temp erature at the interface. Therefore, the proximity of the top plate to the interface affects the temperature gradients at the interface in both phases. The further away the top plat e, or in other words, the greater the ratio of the vapor depth to th e liquid depth, the smaller the temperature gradients in the base state. We learned from our preliminary analysis of the phasechange problem in the absence of fluid dyna mics that the heat transfer mechanism destabilizes the interface. Thus from a heat transfer point of view, the deeper the vapor depth the less the destabilizing effect of the heat transfer. From a fluid dynamic point of view, the deeper the vapor laye r the less the top plate impedes the stabilizing flow in the vapor. Hence, the interface becomes more st able for these two reasons. It is this stabilizing role of the vapor flow that ma kes the passive vapor assumption questionable.

PAGE 51

38 1 2 3 4 5 6 7 8 9 0 5 10 15 20 25 Dimensionless wavenumber, R R4R5R6 Figure 3-3. 4R through 6R vs. wavenumber. Stabilizing Effect of the Vapor Flow In order to understand how the vapor flow st abilizes the system, it is best to look at the velocity profiles at the onset of instab ility of a pure phase-change problem in the absence of gravity a nd Marangoni effects. In Figure 3-4 the velocity profiles are show n for a case where the vapor depth to the liquid depth ratio is taken to be equal to 0.4. The position, 0z0 represents the liquidvapor interface in the base state. To determin e the direction of the flow at a trough (or at a crest) we calculate the value of the vertical component of th e velocity in either phase at the interface relative to the computed value of 1Z. This is done using the eigenvector that corresponds to th e critical temperature differen ce. Then, we find the sign of 1 z1Z v (or 1 z1Z v) where z1v (or z1v) is the liquid (or vapor) ve locity at the interface and 1Z is the surface deflection. This ratio turns out to be negative, implying that when the deflection

PAGE 52

39 is negative, the flow is in the positive directi on. In other words, th ere is up flow at the trough and down-flow at the crest. The flui d flows upset the temperature gradients in both phases at the interface and mo re so in the vapor because of the stronger flow at the interface in the vapor. Once the phase-change takes place upon perturbation, the flow in the liquid will bring the warmer liquid from below attempting to make the trough hotter whereas the stronger flow in th e vapor tries to conv ect this heat away from the trough. Observe the difference in the magnitudes of th e velocities and the natu re of the profiles between the two phases. The hor izontal component of the velo city in the vapor is very strong compared to the liquid velocity. Th e vapor brings the wa rm fluid from the trough to the colder crest, making its temperature go back up again. In ot her words, the liquid flow convects heat to the interface and de stabilizes the interf ace while the heat conduction from the bottom plate to the interface also encourages instability. The vapor flow offers stability though the heat conduction in the vapor encourages instability. The velocity profiles of the perturbed flow tell us that the liquid flow plays a destabilizing role whereas the vapor flow plays a stabilizing role. To give further credibility to the statem ent that the vapor flow is stabilizing whereas the liquid flow is destabilizing, we can perform a little test. By deliberately changing the viscosity of the vapor we can stre ngthen or weaken the vapor flow. Thus, if we increase the vapor viscosity by a small am ount, the vapor flow should be impeded and the system will become unstable. Indeed, an increase in the vapor viscosity destabilizes the interface at all wavenumbers. Accordingly, an increase in the liquid viscosity stabilizes the interface because the destabilizing liquid flow is impeded. It is useful to remember that in both phases the flow is caused by the temperature perturbations at the

PAGE 53

40 interface. This is obvious from the energy e quation at the interface wh erein one sees that perturbed temperature gradients must result in perturbed flows at the interface. This flow at the interface generates pressure pertur bations, which in turn feedback to the temperature through the thermodynamic equilibri um at the interface. To summarize, the stability of a phase-change problem is determined by the competitive nature of heat conduction and heat convection, where the relativ e values of physical properties of both the liquid and the vapor affect the magnitude of the heat transfer mechanisms. The phase-change problem in the absence of gravity becomes unstable at all wavenumbers given that a certain critical temperature di fference is applied across the fluid layers. 0 1 2 3 4 5 6 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 x z Figure 3-4. The velocity vectors in bo th phases at the onset of instability. z0 represents the interface position in the base state. Thus, if we were to conduct experiments in space, we would have to cut-off the small wavenumbers by limiting the size of the co ntainer so that the sy stem is not unstable from the beginning of the experiment. However, when gravity is present, there exists a critical wavenumber at which the interface become s unstable, as depicted in Figure 3-5.

PAGE 54

41 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 0 5 10 15 Dimensionless wavenumber, T (C) Figure 3-5. The critical temp erature difference versus the dimensionless wavenumber for water-water vapor system where d3mm and *d1.5mm. This figure shows the critical temperature drop in C vs. the dimensionless wavenumber, We immediately observe that the critical T is within an experimentally realizable range. In this chapter, we have attempted to explain the physics of the phase-change problem at the onset of instability in a comprehensive manner using a straightforward linear stability calculation. We conclude that the consid eration of the vapor fluid dynamics is vital for the study of stability of the phase-cha nge problem and leaving it out is an invalid assumption. The linear calculati ons tell us what the flow and temperature profiles are at the onset of in stability. This information is enough for us to understand the physics of the instability; however, one needs to suppor t the theory by conducting experiments. In order to know what to expe ct during the course of an experiment we need to turn to the nonlinear analysis. This analysis will give us the information on the type of the instability; in other words, it will tell us whether the instability is subcritical or

PAGE 55

42 supercritical. In Chapter 2, we defined th e input variables to be the fluid depths, the pressure at the top plate and the overall temperature drop between the hot and cold plates. The critical conditions for th e onset of the instability ar e determined along with the critical wavenumber. Advancing the temperat ure drop beyond the crit ical value will not only reveal the magnitudes of the flow and te mperature profiles, but will also give the liquid feed rate and the vapor removal rate, as well as the enhancement in heat transport. In order to find the answer to the behavior of the two-phase problem beyond the onset of instability, we turn to the nonlinea r problem in the next chapter. Endnotes The Exchange of Instability when Heated from the Vapor Side Another set of calcu lations was performed for th e evaporation problem when heated from the vapor side. It is not necessary for the onset mode to be steady. Sometimes the onset modes can also be os cillatory and one should always check the values of the inverse time constants with the leading real part. When the bilayer system is heated from the vapor side and only gravity effect, but not buoyancy-driven and surface tension gradient-driven convection, is ta ken into account it turns out that onset of instability is oscillatory. The Assumption of Flat Surface in the St ability Analysis of the Phase-Change Problem The phase-change problem can also be studied under the assumption that the interface does not deflect with the onset of in stability, i.e., the onset of instability is represented by the transition from the conductive state to the convect ive state. If the interface were artificially constr ained to be flat then surface curvature would be zero and interfacial tension would be infinity making the product a determinate value only upon

PAGE 56

43 introduction of the solved vari ables into the normal component of the force balance at the interface. Under the flat in terface restriction the pure pha se-change problem exhibits a minimum as seen in Figure 3-6. Weak transverse gradients stabilize on the low end, as one needs these gradients in order to have local evaporative regions and adjoining condensation regions. In the absence of tran sverse gradients we would not have local evaporation and condensation. Bo th are needed to have an un stable base state. At the high end of wave number stabilization is on ce again due to viscous and thermal diffusive effects. 0 0.5 1 1.5 0 0.005 0.01 0.015 Dimensionless wavenumber, Ea b Figure 3-6. Evaporation nu mber vs. dimensionless wavenumber for two water-water vapor systems. Curve (a) represents the results for calculations where the interface was assumed to be flat wher eas curve (b) represents calculations with a deflected interface.

PAGE 57

44 CHAPTER 4 WEAKLY NONLINEAR ANALYSIS The Pure Phase-Change Problem In Chapter 3, linear stability methods were used to analyze the phase-change between a liquid and its own vapor where th e vapor, in addition to the liquid, was assumed to be fluid dynamically active. We concluded there that the active vapor layer plays a significant role in determining the stability of the phase-change problem, and offered explanations for the physics behind th is stability. In this chapter, we will examine what happens once we pass the onset of instability. There ar e several reasons to study the post onset region. One reason is to construct the steady solution beyond the onset of the instability and to find the magn itude of the state variables at the given conditions as this is of experi mental importance. In a labora tory experiment the onset of instability can be detected by different m eans, such as IR-imaging, flow tracers or measuring the change in heat transport with respect to a change in the control variable once the critical state is passed. The use of tracers to observe the onset of the instability may cause problems since it could affect the su rface tension at the interface between the liquid and the vapor. Moreover, tracer accumula tion at the bottom of the liquid or at the liquid-vapor interface may present another expe rimental difficulty. IR imaging, though proven to be efficient in flow visual ization of bilayer convection by Johnson & Narayanan (1996), may be difficult because of the likelihood of the vapor condensing on the cold top plate when the bilayer is heated from the liquid side; thereby obstructing the imaging.

PAGE 58

45 The change in the heat transf er at the bottom plate is one way to detect the onset of instability in many convectiv e instability problems (Chandr asekhar, 1961), provided that the change is measurable. Therefore, anot her reason to study the post-onset region via nonlinear calculations is to find the magnitude of the cha nge in the heat transfer. However, there are other adva ntages of doing the nonlinear calculations. Knowing the magnitude of the perturbed domain variables will also tell us how the flows in both domains change when we change any one of th e input variables. Th e linear calculations only reveal the velocities of both fluid phase s relative to each other for a specific set of input variables, whereas th e nonlinear calculations give us the chance to compare different bilayer systems with each other. This way, we can better understand the effect of each input variable or physical property on the instability. The physical model is the same as describe d in Figure 1.1. To get an idea on the nonlinear features, we will investigate the phase-change problem in the absence of gravity and Marangoni convec tion in this chapter. The problem is studied in tw o spatial dimensions, x and z. The nonlinear equations in dimensionless form for the pure phase-cha nge problem are repeated below for the sake of convenience. The momentum equations in each phase are 21v +vv-Pv Prt 4.1 and ***2* **1v +vv-Pv Prt 4.2 Recall that the energy equations in each domain are given by

PAGE 59

46 2T T+vT t 4.3 and 2*** *T T+vT t 4.4 Also remember that the continuity equati ons in each phase hold assuming that the fluids are incompressible in the sens e of the Boussinesq approximation. We will complete the mathematical model by repeating the interfacial conditions at the phase-change boundary. The mass balance at th e interface, (,)zZxt is given by *vunvun 4.5 At the interface, the tangential components of velocities of both fluids are equal to each other, so t v t v* holds. The interfacial tension come s into the picture through th e force balance. This balance equation has normal and ta ngential components and is given by ** ***11 Cav vuTn2H nCav vuTn PrPr 4.6 In addition, we need the inte rfacial energy balance, viz. 2 2 ** PC * PC11k 1KvuvuvunETnTn 22k VSvuSvun0 4.7

PAGE 60

47 At the interface, the temperatures of both fluids are equal to each other and the thermodynamic equilibrium at the phase-change boundary is * KEBASE BASET PPln T 4.8 In the nonlinear analysis one goes beyond th e critical state by advancing a control variable such as the temperature drop acro ss the plates. The scaled temperature drop, however, appears as E, the Evaporation number, therefore it is E that will be advanced from its critical value. Once E is increased all of the state variables will respond and it is their response that interests us. Whether the response is linearly proportional to the change or not is also of interest and w ill be determined. To this end we introduce perturbation expansions a bout the base state for each response state variable. Perturbation Expansions The domain variables and their derivativ es are expanded along the mapping as follows: 2 2 2 000 1 0C11C2112 2 0000 23 2 3 23 000 211 C31212113 223 000000uuu u 1 uuuzu2zzz z2zzz uuu uuu 1 u3z3z3zz3zzz... 6zzzzzz and 232 2 2 2 0000 121 C1C112 2232 0000000 342 223 3 23 3000 211 C1212113 223342 0000000uuuu uuu u1 z2zzz zzzz2zzzz uuuu uuu 1 3z3z3zz3zzz... 6zzzzzzz In the above expansions, 1z, 2z and 3z specify the mappings from the current state in the (x,z)-domain to the refe rence quiescent state in the (x0,z0)-domain at first, second

PAGE 61

48 and third orders, respectively. High er derivatives with respect to z0 and derivatives in other directions follow in a like manner. It has been shown by Johns and Narayanan (2002) that the derivatives of the spatial mappings 1z,2z,3z etc. disappear due to algebraic cancellations and never need be of any concern. In the above expansions, is the dimensionless control variable, and C is its value at which the system becomes unstable, namely the critical value. In our problem will be replaced by the parameter E. We go on and set the t terms equal to zero and invest igate the nature of the steadystate solutions to the problem once we cross the threshold value. To do this, we will go beyond the critical point by an amount i.e., C The control variable is now defined as CEE We will substitute the expansions of the domain variables into the nonlinear equations and find the values of along with the amplitude of the perturbations that excited the instability, us ing the method of dominant balance at every order of the calculation as given by Grindrod (1996). The equations are very long and the process of balancing the orde rs is tedious; thus, in order to assist the reader, we give an overview of what we are going to do in the calculation Overview Suppose that the steady nonlinear problem is in the form of Nu=f, where N is a nonlinear operator containing the control variable u is the solution vector and b is an inhomogeneity vector that driv es the problem. The vector u contains the velocity components, the pressure and the temperature, in addition to the surface position, Z, which is the last element of the vector u. The base problem must satisfy the nonlinear equation and its solution is triv ially known. In our applicatio n, the base state is a pure

PAGE 62

49 conductive state, so in each flui d the velocity is zero and the temperature field is linear. The first order problem is required to be an eigenvalue problem and this means that cannot be determined at this order. Th e first order problem is written as: 1Lu=0. Here, the operator L derived from the O( ) equations by linearization of N and contains the control variable The eigenvector u1 is, of course, known only up to an arbitrary multiplicative constant, A, which is the amp litude of the disturbances that cause the instability. At first order we can determine u1 and the critical value of but not and A. It is our aim to determine what A, as well as is through the dominant balance method, and therefore we proceed to the next order. The second order problem is in the form 2111Lu=fu,u. The solvability condition is applied to this system of equations. But, in our problem the operator L is not self-adjoint. Hence, to make the algebra simpler, the elements of the vector u2 are solved in terms of its last element, i.e., Z2 and 2 1Z, where Z2 is the second order correction to the surface posi tion Z. Solvability is then applied to the differential equation for Z2, which is actually the normal stre ss balance at the liquid-vapor interface as the differential operator for Z at ev ery order is self-adjoint. The normal stress balances in the first and second orders look as follows; 2 1 1 2 0dZ Z0 dxa and

PAGE 63

50 2 2 2 221 2 0dZ ZZ dxab We find in this problem that, because Z1 is a multiple of cos(x) the solvability condition is automatically satisfied at the s econd order; thus we m ove on to the third order problem. In our problem, at third or der, we find that a balance can be struck between the various terms for a specific value of We then introduce the specified value of into the equations and apply the so lvability condition on the normal stress balance, which is now in the form of 2 3 3 33112 2 0dZ ZZ,ZZ dxab By doing so, we get an equation th at determines the value of A2. Hence, in the dominant balance method, we find the value of and A at the same order. The value of tells us what type of bifurcation we should expect to see once the instability sets in. For example, if is equal to one the bifurcat ion is trans-critical; and if is equal to it is a pitchfork. To determine whether the p itchfork is forward or backward we focus on the sign of A2. If it is positive the pitchfork is forward, else it is backward. In the latter case, E must be expanded as CEE The work showing how to determine is shown in detail as an endnote to this chapter. Perturbed Equations Using the method outlined above, we expa nd the nonlinear domain equations given earlier up to the third order. Thus far we have pretended that our physical system is infinite in lateral extent, but the pure phase -change problem is unlike the standard Bnard problem, which produces a nonzero critical wavenu mber. For that reason, in this chapter

PAGE 64

51 we will choose an arbitrary wavenumber, a nd will proceed with the calculations using that wavenumber. As a result, we have tacitly assumed that the fluids are confined within a laterally finite system with slippery and insulating walls. An analysis for long wavelength modes assuming a passive vapor has been considered by Oron et. al. (1997). Returning to the expansi ons of the governing equations, we first write the momentum equations in both phases as 2 22 C0101C0202101 3 2 C030310220112 EEPvEEPvvv 2Pr 13 EEPvvvvv0 6Pr 4.9 and 2 *2**2*** C0101C0202101 *** 3 *2***** C0303102201 **12 EEPvEEPvvv 2Pr 13 EEPvvvvv0 6Pr 4.10 The energy equations in each phase are 2 222 00C01100C02200101 3 2 C033002011021 TEETvTEETvT2vT 2 1 EETvT3vTvT0 6 4.11 and 2*2*** 00C01100 2 2***** C02200101 ** 3 2******* C03300201102 **TEETvT 1 EETvT2vT 2 1 EETvT3vTvT0 6 4.12

PAGE 65

52 The continuity equations become 23 C01C02C0311 EEvEEvEEv0 26 4.13 and 23 *** C01C02C0311 EEvEEvEEv0 26 4.14 The unit normal and the unit tangential vect ors, upon perturbation up to the third order, become 2 2 121 CC 000 3 3 3 112 C 0000ZZZ 1 nEEEE x2xx Z ZZZ 1 EE33 6xxxx kiik ik 4.15 and 2 2 112 CC 000 3 3 3 121 C 0000ZZZ 1 tEEEE x2xx Z ZZZ 1 EE33 6xxxx ikik ik 4.16 At the lower plate, 1 z, the boundary conditions vn0 vt0 and T1 are expanded as 23 Cz1Cz2Cz311 EEvEEvEEv0 26 4.17 23 Cx1Cx2Cx311 EEvEEvEEv0 26 4.18 and

PAGE 66

53 23 0C1C2C311 TEETEETEET1 26 4.19 while the conditions at the top plate *vn0 *vt0 and *T0 become 23 *** Cz1Cz2Cz311 EEvEEvEEv0 26 4.20 23 *** Cx1Cx2Cx311 EEvEEvEEv0 26 4.21 and 23 **** 0C1C2C311 TEETEETEET0 26 4.22 For reasons of brevity, we do not provide al l of the expansions of all interfacial conditions save three that are crucial to the an alysis of the problem. One of those three conditions is the interfacial en ergy balance and, after simplifi cations, it assumes the form

PAGE 67

54 * ** 00 11 Cz1 0000 z11 z21x1 00 2 2111 1 2 2 0000 C *2** 2111 1 2 0000TT TT kk EEEvE zkzzkz vZ v2Z2v zx TTTZ 2Z2 zzxx 1 EEE 2 TTTZ k 2Z2 kzzxx * z1z1 PC 00 2 2 2 1z2z1z1 z3z1121 2 0000 21x11 x1x21 0000 22 PCz1z1x1z 3 Cvv 2V zz Zvvv v3v3Z3Z3Z xzzz ZZvZ 3v3v6Z xxzx 3Kvvvv 1 EE 6 *2*2 1z1 223 2 3 211 121 223 0000 2 122111 1 0000000 2*2*3* 2 3 211 121 223 0000 **2* 12211 1 000000v T TTT 3Z3Z3Z zzzz TZTZTZ E336Z xxxxzxx T TTT 3Z3Z3Z zzzz k k TZTZTZ 336Z xxxxzx 1 0 2 z2z1x11z1 z11 2 00000 z11z1z1 PCx1z21 0000 *** z2z1x1 z1 00x vvvZv 6v22Z zxzxz vZvv V12v6v2Z zxzz vvv 6v2 zxz 2* 1z1 1 2 000 *** ** z11z1z1 x1z21 0000Zv 2Z xz vZvv 12v6v2Z zxzz 0 4.23

PAGE 68

55 Another one is the thermodynamic equilibrium relation between the vapor pressure and the vapor temperature; * 0 KE0sat sat ** 0 CKE111 sat0 * 1 KE21 0 2 C 22 ** ** 00 1 21211 sat00sat0 K 3 CT PPln T dT 1 EEPTZ Tdz P P2Z z 1 EE 2 dTdT T 11 T2ZZTZ TzdzTdz 1 EE 6 **2* *2 211 E3121 2 000 **2* *2 0 211 31213 2 sat0000 2 ** ** 00 1 21211 sat000 3 * 0 11 sat0PPP P3Z3Z3Z zzz dT TTT 1 T3Z3Z3ZZ Tzzzdz dTdT T 1 3T2ZZTZ Tzdzdz dT 1 2TZ Tdz 30 4.24 We also give the third order normal stress balance equation because we will use it to find the values of and A. The third order norm al stress balance is given by

PAGE 69

56 00 *2 * z1z11 C11 2 000 * 11 2121 00 2*2* 2 z2z1z2z1 C11 22 0000 2 z1z1z1CaPP vvZ EECaPP22 zzx PP P2ZP2Z zz Ca vvvv 1 EE22Z22Z 2zzzz Ca 2vvv Pr 2 2 0 2 2 211 3121 2 000 **2* *2 211 3121 2 000 2 223 2 z3 z2z1z1z11 121 223 000000 z 3 CZ x PPP P3Z3Z3Z zzz PPP P3Z3Z3Z zzz v vvvvZ Ca23Z3Z3Z12 zzzzzx v 2 1 EE 6 2 2*2*3** 2 3 z2z1z1z11 121 223 000000 z11 z21x1z1z1 00 z11z z21x1z2 00vvvvZ 3Z3Z3Z12 zzzzzx vZ Ca 3v2Z2vvv Przx vZv Ca 3v2Z2vv2 Przx * 21 1x1z1 00 2 2 2 3 11 22 0000 Z Z2vv zx Z ZZ 9 xxx 4.25 Remember that in all of the expansion formulas, E will be replaced by EC+ The first order problem is an eigenvalue problem a nd it was solved in an earlier chapter. The solution to the first order problem is used at every subsequent order along with the solvability condition. We learn from the first order equations that the eigenvector u1 in terms of Z1 is

PAGE 70

57 1110ˆ ˆ uA u(Z)cosx Recall that is an input and it is determined by the width of the compartment in the lateral direction. We cannot determine the value of either from the first or the second order equations. As noted in the ea rlier section, the solvability condition is automatically satisfied at the second orde r. The solution to the second order, u2, will go as 2222 22a202b12c10ˆˆˆ uB u(Z)cosx+Au(Z)Au(Z)cos2x It is at the third order that we learn that is equal to 2 1 At this order, there exists a value for A other than zero. A appears as A2, and the value of A2 is obtained using the solvability condition on the third order normal stress balance. The third order solution, u3, depends upon the form of u1, u2P as well as Z3, i.e., 33 33a303b1120 33 3c11203d103e120uC u(Z)cosxAu(Z,ZZ)cosx Au(Z,ZZ)cos3xA u(Z)cosxB Au(ZZ)cos2x1 Results of Calculations The Chebyshev spectral tau method is used to solve the eigenvalue problem that determines the critical temperature difference and the value of the unknown multiplicative constant A. Once we know what A is, we can calculate the actual change in the heat transfer rate at the bottom plate when we advance T a little beyond the critical point. Our system is finite in the lateral extent with periodic conditions. The heat transfer up to the second order is calculated as

PAGE 71

58 1 0 122 2 000 0000 00ˆ ˆˆˆ dT dTdTdT 1 qnkcosxcos2xdx dzdz2dzdz 4.26 Clearly, 0cosx and 0cos2x when integrated over the width, give zero. We are more interested in the percentage of change of the heat transfer, so we have plotted 0 00qqn qn against CE E in Figure 4-1 for a system where the liquid and vapor depths are 4mm and 0.4mm, and the wavenumber is 0.3 In nonlinear calculations the wavenumber of the disturbance is fixed. Accordingly, one can find the width of the box using the wavelength of the disturbance. The wavenumber, 0.3 corresponds to a cell width of about 84mm. We learn that once the interface is unstable, more heat is transferred through the bottom plate to the system. Thus, if we wanted to know when the system becomes unstable, we only have to measure the amount of heat that is needed to keep the bottom plate at a constant temperature. The heat ne cessary changes with the onset of instability. For a typical problem we can observe a 2% in crease for a supercriti cal change of about 10%.

PAGE 72

59 0.98 1 1.02 1.04 1.06 1.08 1.1 -0.5 0 0.5 1 1.5 2 E/EcPercentage increase in the heat transfer Figure 4-1. Heat transfer across the bottom plate. The nonlinear analysis results do more than giving the heat transfer enhancement. They can also be of use in helping us to better understand the physics of evaporative instability. To see how, let us first observe fr om the linear calcula tions that the deeper vapor layers make the interface more stable. Then apply the nonlin ear calculations to two cases with the same base state interface pr essure and temperature where in one of the systems 0.1 and in the other 0.3 and compare the magnitudes of the pressure perturbations in the vapor at the interface. Because 1 *** 2 01PPP at the interface and because we wish to compare the problems w ith each other for the same percentage of departure from the critical state, we will express *P as 1 2 1 *** C 2 0C1 CEE PPE P E

PAGE 73

60 In Figures 4-2 and 4-3, we have plotted th e magnitudes of pressure perturbations in the vapor, 1 2 C1E P, for the cases where 0.1 and 0.3 We learn from the calculations that even though the nature of the profiles of the vapor pressure perturbations are the same for both cases, the magnitudes are not. The vapor pressure perturbation at the interface for the case with the shallow vapor layer (0.1 ) is about three times more than the pressure perturbation at the interface for the case with 0.3 resulting in a more unstable configuration. In Chapter 3, during our discussion of the physics of the pure phase-change problem we had deliberately changed the value of the vapor viscosity to see its effect on the instability. The sy stem, where the vapor viscosity was increased, proved to be more unstable. -7.174 -7.172 -7.17 -7.168 -7.166 -7.164 -7.162 x 10-3 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 Ec 1/2 P1 *z Figure 4-2. The scaled pre ssure perturbations in the vapor domain for the case where d1 mm 0.1 and 1 z0 represents the interface position in the base state.

PAGE 74

61 -2.23 -2.22 -2.21 -2.2 x 10-3 0 0.05 0.1 0.15 0.2 0.25 0.3 Ec 1/2 P1 *z Figure 4-3. The scaled pre ssure perturbations in the vapor domain for the case where d1mm, 0.3 and 1 z0 represents the interface position in the base state. As shown in Figure 4-4, our nonlinear cal culations return a larger pressure perturbation value at the interface for the syst em with the more viscous vapor. Thus, the nonlinear calculations support the argument that the va por flow dynamics affect the stability of the phase-change problem.

PAGE 75

62 -3.28 -3.27 -3.26 -3.25 -3.24 -3.23 x 10-3 0 0.05 0.1 0.15 0.2 0.25 0.3 Ec 1/2 P1 *z Figure 4-4. The scaled pressu re perturbations in the more viscous vapor for the case where d1mm, 0.3 and 1 z0 represents the interface position in the base state. Endnote: The Determination of There is a lot of analytical work th at goes into determining the value of For that reason, we will demonstrate this cumb ersome procedure on a simpler problem. Assume the following domain equation 2 2du f(u)0 dx 4.27 where f(0)0, 'df f0(0)0 dx and 2 '' 2df f0(0)0 dx are given, and the domain equation is subjected to Dirichlet condition at x0 and x1 To determine the behavior of the system around the critical point, in this case c we will go a little beyond the critical point and replace by c When we expand the domain variable u around a base state 0u as a power series of the domain equation can be written as

PAGE 76

63 2 23 123 2 2 '23''23 c123123 '23'' 1231d11 uuu dx26 11111 f(0)uuuf(0)uuu... 26226 1111 f(0)uuuf(0)u 262 2 23 231 uu...0 26 The order problem is then 2 1 c1 2du f0u0 dx 4.28 subject to 1u00 and 1u00 When we solve for 1u, we find that 1uA sin(x) and 2 c 'f0 Now we will solve for 2u and 3u and so on. To do that we will move to the next order of equations, which is unknown since is yet unknown and needs to be determined. There exist 3 possibilities: 12 is one of the possibilities, and this implies that 1f0u0 Now, this can only be true if only if 1u0 so obviously this is not a plausible choice. Another possibility is that 12 When we set 1 we get the following equation as the next order equation: 2 '''2' 2 c2c11 2du 111 f0uf0uf0u0 2dx22 4.29 If we multiply equation 4.28 by 2u and equation 4.29 by 1u, then integrate both equations over the domain, and subtract them from each other we get

PAGE 77

64 11 ''3'2 c11 001 f0u dxf0u dx0 2 Using the relation above, we can determine a value for A which is non-zero, and it is 2 ''f0 3 A 4f0 So far, one of the three possibilities give s us a non-zero value for A. The next possibility is that 12 Now this gives us 1 ''3 c1 01 f0u dx0 2 This is not possible since 1 3 1 04 u dx 3 Thus, we conclude that 1 Note that the value of A is determined at the same order where we determined the value of A different approach to the weak nonlinear analysis There is more information one can actually get from the weak nonlinear analysis. To do that, we will have to modify our physical model. This time we will again set the evaporation rate in the base state equal to zero by adjusting the vapor pressure, however, now we will leave the feed and removal valves open, and as we increase the temperature difference applied across the bilayer system we will adjust the vapor pressure accordingly such that the evaporation rate is equal to zero in the base state even if we cross the critical point. Intuitively, one expects the net evap oration rate to change once we cross the critical point. Beyond the critical point the net evaporation rate, U, can be expanded as b aseextraUUU 4.30

PAGE 78

65 where extraU is extraflatUU 4.31 In equation 4.31, flatU is the change in the net evapora tion rate if the instability had not occurred, and is the correction to it because of the ons et of instability. In this analysis, we will adjust the vapor pressure as we in crease the temperature difference applied, accordingly, flatU is always zero even after we cross the critical point. Thus, the only new unknown in the problem is The sign of tells us what will happen to the net evaporation rate beyond the critic al point. If it turns out to be positive, then with the onset of instability the system will suck in liquid through the bottom valve and remove the excess vapor through the top valve. On e does need to solve the complete weak nonlinear problem to get this information. Moreover, one can even determine the sign by analytical calculations. From our earlier analysis we know that the solvability condition was automatically satisfied at the s econd order. That means that we can go ahead and solve the second order problem. By doing that we will find in terms of the amplitude of the deflection. The momentum equations in each phase at the second order are 2 02021012 Pvvv Pr 4.32 and *2*** 0202101 **2 Pvvv Pr 4.33 The energy equations in each phase are

PAGE 79

66 2 02200101TvT2vT 4.34 and 2***** 02200101 **TvT2vT 4.35 The continuity equations become 02v0 4.36 and 02v0 4.37 At the lower plate, 1 z, the boundary conditions vnU vt0 and T1 are expanded as z22vU x2v0 and 2T0 while the conditions at the top plate ** 2vnU *vt0 and *T0 become ** z22vU x2v0 and 2T0

PAGE 80

67 For reasons of brevity, we do not provide al l of the expansions of all interfacial conditions save three that are crucial to the an alysis of the problem. One of those three conditions is the interfacial en ergy balance and, after simplific ations, it assumes the form * 22z11 z21x1 0000 22*** ** 111111z1z1 11PC 22 00000000TTvZ k vE2Z2v zkzzx TTZTTZvv k E2Z22Z22V zxxkzxxzz 4.38 Another one is the thermodynamic equilibrium relation between the vapor pressure and the vapor temperature; * ** 0 1 KE2122 0base0 22 * 0 1 111 ** base0base0dT P 1 P2ZTZ zTdz dT T 11 2ZTZ TzTdz 4.39 The normal stress balance assumes the following form *2 * z2z22 22 2 000 *22* * 11z1z1 1111z1z1z1 22 0000vvZ CaPP22 zzx PPvvCa Ca2Z2Z4Z4Z2vvv zzzzPr 4.40 When one solves the first order prob lem, one determines that the domain variables are actually trigonometric functions in the 0x-direction. They are, in the liquid, x1x100ˆ vvzsinx z1z100ˆ vvzcosx 1100ˆ TTzcosx and

PAGE 81

68 1100ˆ PPzcosx The vapor domain variables assume similar fo rms. We plug them into the second order domain equation and find the form of the s econd order domain variables. What follows is the demonstration of how we find the form s of the second order domain variables. The second order energy equation is 2 02200101TvT2vT After substituting the first order domain variables we get 2 0 11 02z2x1z1 000dT TT Tv2vv dzxz which translates to 222 0 1 02z2x10100z1000 00ˆ dT dT ˆ ˆˆ Tv2vzTzsinxvzzcosx dzdz Using trigonometric manipulations we get 2 0 1z1 02z2z100100 000 1z1 z1001000 00ˆ ˆ dT dTdv ˆ ˆ TvvzzTzz dzdzdz ˆ ˆ dTdv ˆ ˆ vzzTzzcos2x dzdz From the manipulations above, one learns that 220200ˆ ˆˆ TTzTzcos2x Following the same method, one obtains x2x20x200ˆ ˆˆ vvzvzsin2x z2z20z200ˆ ˆˆ vvzvzcos2x and

PAGE 82

69 220200ˆ ˆˆ PPzPzcos2x We are trying to find the sign of At the second order, is actually equal to the x independent part of the correction to the net evaporation U, which is 2ˆ U. Thus, we only need to solve for the x-independent part of the particular solution for 2U. The domain equations in the liquid for x-independen t part of the second order problem are 2 x2 2 0ˆ dv 0 dz 2 2z2z1 z1 2 000ˆ ˆˆ dPdvdv 2 ˆ v dzdzPrdz 0 21z1 z2z11 0000ˆˆ ˆ dT dTdTdv ˆ ˆˆ vvT dzdzdzdz and z2 0ˆ dv 0 dz In the vapor domain, they are 2* x2 2 0ˆ dv 0 dz *2** 2z2z1 z1 *2* 000ˆ ˆˆ dPdvdv 2 ˆ v dzdzPrdz *** *** 0 21z1 z2z11 ** 0000ˆˆ ˆ dT dTdTdv ˆ ˆˆ vvT dzdzdzdz and z2 0ˆ dv 0 dz

PAGE 83

70 At the bottom plate, the boundary conditions are x2ˆ v0 z22ˆ ˆ vU 2ˆ T0 At the top plate, x2ˆ v0 ** z22ˆ ˆ vU 2ˆ T0 2ˆ P0 At the interface, we have * z2z2ˆˆ vv 4.41* x2x2ˆˆ vv 4.42 * 00 11 22221 1 0000ˆˆ dTdT dTdT ˆˆˆˆˆ TZTZZ dzdzdzdz 4.43 * x2x2 00ˆˆ dvdv dzdz 4.44 ** 0 22KE2 base0 22 ** 0 11 111KE1 ** base0base00dT 1 ˆˆˆ TZP Tdz ˆˆ dT dTdP 11 ˆˆˆˆ ZTZZ TdzTdzdz 4.45

PAGE 84

71 ** 0 22 z1 z2z1z11PCz1 0000ˆ ˆˆ ˆ dT dTdTdv k ˆ ˆˆˆˆ vEEvvZ2V1v dzkdzdzdz 4.46 2*2 *** * z2z22 22 z1 2 000ˆ ˆˆ dvdvZCa ˆˆ ˆ CaPP221v dzdzxPr 4.47 These equations allow an analytical soluti on for the problem. First, we learn that everywhere in each domain x20ˆ vz0 x20ˆ vz0 z202ˆ ˆ vzU z202 *ˆ ˆ vzU We also learn that 2 20z11 ˆ ˆ Pzv Pr 2* ** 20z11 ˆ ˆ Pzv Pr For reasons of brevity, I will omit giving the analytical solution for 20ˆ Tz and 20ˆ Tz. In conclusion, we find the corr ection to the ne t evaporation, 2ˆ U, and it turns out to be positive at each wavenumber. This information is useful if we were to do experiments since now we have an additional indicator, othe r than measuring the heat transfer at the bottom wall, to determine wh ether the instability has set in or not. Opening the valves and allowing the bilayer system such in liquid through the bottom

PAGE 85

72 wall with the onset of in stability is a preferable method si nce it is a visual method, unlike measuring the change in the heat transfer.

PAGE 86

73 CHAPTER 5 THE RAYLEIGH-TAYLOR PROB LEM WITH PHASE-CHANGE Introduction As we discussed earlier in Chapter 1, a vapor underlying a li quid results in a configuration that is potentia lly unstable. This problem, known as the Rayleigh-Taylor problem, is well known and an account of its underlying physics is given by Chandrasekhar (1961). It may be recalled from Chandrasekhar that the interface is unstable for wavenumbers smaller than a di mensionless critical wavenumber given by 1 2 lg2 RTdg w Here, l and g are the liquid and vapor densities, respectively, is the interfacial tension, d is the horizontal length scale and g is the coefficient of gravitational acceleration. The derivation of this result is also given in Appendix D. The classical Rayleigh-Taylor instability problem is ordina rily studied under the assumption that there is no mass transfer across the interface. The problem was reformulated by Hsieh (1972) who considered the Rayleigh-Ta ylor problem with mass and h eat transfer in inviscid fluids of finite vertical depths. He solved both the velocity and the temperature fields in both fluids and derived an analytical dispersion relati on i.e., a relation between and Later, Hsieh (1978) studied the same problem again for inviscid fl uids but introduced a phenomenological coefficient, which allowed him to decouple the velocity and temperature fields for both infini tely deep fluids as well as for fluid layers with finite

PAGE 87

74 depths. Using the same idea of a phenomenolog ical coefficient, Ho (1980) analytically obtained a result for viscous fluids with finite depths. In the same paper, Ho also derived the critical dispersion relation for viscous fl uids with equal kinematic viscosities for infinitely deep fluid layers. As an applica tion to pool boiling, Badr atinova et al. (1996) have found an analytical expression for the critical wavenumber for a system where a vapor with finite depth underlies its own in finitely deep liquid, solving for both the velocity and the temperature fields. Howeve r, we have seen in Chapter 3 that the assumption of large depth fluids may lead to conclusions that do not quite represent the physics. We might also like to know what the phenomenological coefficient used by Ho (1980) and other workers such as Adham-Khodapa rast et al. (1995) re ally means. These workers found a critical disper sion relation for infinitely deep fluid layers with equal kinematic viscosities to be given by l 22 RT40www Here, l is the kinematic viscosity of both the liquid and the gas phases. Note that upon setting the phenomenological coefficient equal to zero in the above relation, one recovers the classical Rayleigh-Taylor instability result. In this chapter we therefore investig ate the Rayleigh-Taylor problem with interfacial mass and heat transf er and present the results of our calculations for viscous fluids where both fluid layers ha ve finite depths. In additi on, we derive an analytical expression for the phenomenological coefficient, albeit for infinitely deep viscous fluid layers and show that this coefficient is actually a function of the wavenumber. Physical and Mathematical Model Unlike the earlier chapters this time the physical model is a system where a vapor of depth, dg, underlies its own liquid of depth, dl, and is heated from below (Figure 5-1).

PAGE 88

75 Figure 5-1. The physical model As before the base state is one where ther e is a flat interface between the liquid and its underlying vapor, and the stability of this base state is in question. Also as before the liquid and vapor depths are a ssumed constant by suitably adjusting the vapor removal and liquid feed rates. The phase-change rate at the interface can be controlled by adjusting the vapor pressure at the vapor-side wall. The phase-change rate is not constant as it can change upon perturbation. Consequently, in this problem as in the pure evaporation problem the input variables are the liquid and vapor depths, the temperature difference between the bottom and top plates and the pressure at the lower plate. Furthermore, the pressure at the lower plate can be adjusted su ch that the phase-change rate at the interface is equal to zero. In the analysis of the problem, we assume as usual that the liquid and its vapor are in thermodynamic equilibrium because we will only consider either zero or very small phase-change rates at the interface in the base state. The equations that model the physics are given by the Navier-Stokes, energy and continuity equations in each phase assuming that the fluids are incompressible. It is assumed that the bottom and top plate temperat ures are kept constant. When we assume that the liquid and its vapor are in equilibri um and that no phase change is taking place across the interface, we need not feed a nd remove fluid through the bottom and top plates. The no-slip cond ition applies along the top a nd bottom plates. Hereon, xv, zv,

PAGE 89

76 and T are taken to be the x and z-compone nts of velocity, and temperature fields, respectively, while the asteri sk denotes the vapor phase. The equations are brought into dimensi onless form using th e following scales: length, ll d ; velocity, lv d ; time, v d where l and l are the liquid viscosity and the liquid thermal diffusivity, respectively. The temperature is scaled as g 0 0T T dT d dz In the temperature scale, g 0 0dT dz is base state temperature gr adient in the vapor at the interface. Linearized stabili ty is performed in the standard way by expanding the state variables in the form of 00tx 0ˆ uuziwee leading to the following dimensionless equations of motion 22 22l z1 22dd v0 dzdzPrww 5.1 and 22 22g z1 22dd1 v0 dzdzPrww 5.2 where and Pr are the inverse time constant and liquid Prandtl number, respectively, is the ratio of the vapor kinematic viscos ity to the liquid kinematic viscosity and w is the dimensionless wavenumber of the perturbation. The perturbed energy equations in dimensionless form become 2 2ll z1z1 2d Tk v dzw 5.3

PAGE 90

77 and 2 2gg z1z1 2d1 Tv dzw 5.4 In the equations above, k is the ratio of the vapor thermal conductivity to the liquid thermal conductivity, and here is the ratio of the vapor thermal diffusivity to the liquid thermal diffusivity. At the lower plate, the no-flow, no-slip and constant temperature imply that the vapor is subject to g z1v0 g z1 0dv 0 dz and g 1T0 while similar conditions at the top plate become l z1v0 l z1 0dv 0 dz and l 1T0 At the interface the mass balance is given by lg z1z11v v 1 Z 5.5 where 1Z is the perturbed surface deflection, and is the ratio of vapor density to liquid density. The continuity of temperature and the noslip conditions at the interface become lg 1111Tk ZTZ and gl z1z1 00dvdv dzdz The normal and tangential stress balances assume the following dimensionless form: 3gg3ll 22222 z1z1z1z1 RT1 33 0000dvdvdvdv 1 33Z0 dzPrdzdzPrdzwwwww 5.6 and

PAGE 91

78 2g2l 2g2l z1z1 z1z1 22 00dvdv v v0 dzdzww 5.7 where is the ratio of the vapor viscosity to the liquid viscosity. The interfacial energy balance is lg l 11 z11 00dTdT v ZEk0 dzdz 5.8 where again E stands for the Ev aporation number and is given by l pCT E Here, l p C and are the liquid heat capacity and the la tent heat of evaporation per unit mass, respectively. The equilibrium condition at the interface becomes 3gg 222g z1z1 KEPE111 3 00dvdv 1 ZTZ dzPrdzwww 5.9 where again g KE 2 d and PEg d are two dimensionless parameters from the linearized Clapeyron equation, and is the inverse of the scaled saturation temperature at the flat interface. Analytical Results for In finitely Deep Fluids Ordinarily the fluid layers are assumed to be of finite depths but to provide the reader with an expression for the phenomenological coefficient, by way of an analytical procedure and avoid a numerical calculation we assu me that they are infinitely deep. The Rayleigh-Taylor problem is essentially an interfacial instability problem, so to find the leading inverse time constant, one can omit the ‘’s coming from the domain equations and domain vari ables, and keep the ‘’s that appear only along with the surface

PAGE 92

79 deflection. Doing this makes it easier to find an analytical di spersion relation for When we subject the vapo r and liquid domain equati ons to the conditions at 0z and 0z, respectively, we get the following expressions for the velocity and temperature fields. 00zz lll z10vA B z wwee 00zz ggg z10vA B z wwee 000l zzz llll2 100kB1 TC Az Bz 222wwweee ww 000g zzz gggg2 10011B1 TC Az Bz 222wwweee ww The steps taken to get the analytical resu lt are given below. Using the no slip condition and the ta ngential stre ss balance one gets llggABAB ww and llggABABww. We deduce from these two equations that llABw and ggAB w

PAGE 93

80 Now we can write two of the unknowns, Bl and Bg, in terms of Al and of Ag. We then plug in the expressions for l z1v and g z1v into the normal stress balance and we get lg22 RT12A2AZwwww. The mass balance gives us lg 1AA1Z We can use the relation obtained from the interfacial mass balance to eliminate Al from the normal stress balance and we get Ag in terms of the firs t order perturbation of the surface deflection, Z1, which is 22 RT g 121 AZ 2www w When we plug the expressions for l 1T and g 1T into the continuity of temperature and the thermodynamic equili brium relation we get gl 11CZCk Z and gg 1KEPE1CZ2A Zw. Now we are in a position to write lC and gC in terms of Z1. We still have not used the interfacial energy balance. When we pl ug everything into the energy balance at the interface we get 1....Z0 For this equality to hold either Z1 or the term in the parenthesis must be equal to zero. Z1 equal to zero is not an option since then there would be no solution to the

PAGE 94

81 problem solved. Setting the terms in the parenthesis equal to zero we get the following expression as the dispersion relation: FG H 5.10 where 22 2 RT KE1 3k 2 1k E 24 ww w F ww PE(1k) E2k Gw and 2 KE1 3k 12 1k 3k E1E 44 w H ww When we set 0 we get an expression for the cr itical wavenumber in terms of the vapor temperature gradient at the interface and the physical properties of both fluids. We mentioned earlier that Ho (1980) and Adham-Khodaparast et al. (1995) derived a dispersion relation for the phenomenological coefficient, by decoupling the temperature and velocity fields. Once is introduced only the ve locity fields in both fluids were obtained. To obt ain the result of Adham-Khoda parast et al., one only needs to drop all the energy equations from our model and replace the interfacial energy balance by the following balance equation in dimensional form: l z111vZZ

PAGE 95

82 Simply put, the phenomenological coefficient is introduced in the interfacial energy balance where it replaces the heat conducti on terms. The dispersion relation of AdhamKhodaparast et al. using dime nsionless equations becomes lg 22 RTd 20 www When we assume that the kinematic visc osities of both fluids are equal to each other, we get the familiar result of Adham-Khod aparast et al. in dimensionless form. It is l 22 RTd 40www Going back to analytical dispersion rela tion, when we rearrange the equation for the critical wavenumber, we get the expr ession for the phenomenological coefficient, at the onset of instability. It is gl 2 2 KEd 1 3k 2 1k E 4 G w w We see that at the onset of instability, itself is a function of the wavenumber as well as the physical prop erties of the fluids. In Figure 5-2, the dimensionles s phenomenological coefficient, is plotted against the dimensionless wavenumber to se e the shape of the curve and investigate whether there is any range of w in which remains roughly unchanged. The calculations are done assuming th e physical properties of wa ter and water vapor at its

PAGE 96

83 saturation temperature under 1 atm, i.e., at 100 C. The physical properties used in the calculations are given in Table 3-1. Figure 5-2 shows the strong dependence of the phenomenological coefficient on the wavenumber in the plotted wavenumber ra nge. The classical Rayleigh-Taylor dimensionless critical wavenumber for the same system is 73.4 e We observe that at exactly the classical Rayleigh-Taylor critical wavenumber the value of the phenomenological coefficient is equal to zero. When we look at the relation derived by Adham-Khodaparast et al., it is easier to see why 0 when RTww 0 0.5 1 1.5 2 2.5 3 3.5 x 10-7 0 1 2 3 4 5 6 7 8 9 10 x 106 Dimensionless wavenumber, Phenomenological coefficient, Figure 5-2. The dependence of the phenom enological coefficient on the wavenumber The graph therefore tells us that the use of a phenomenological coefficient to decouple the temperature and velocity fields is not a good assumption as the phenomenological coefficient its elf depends strongly upon the wavenumber. In fact, the coefficient even turns negative for large wave numbers.

PAGE 97

84 Numerical Results for Finite Depth Fluids To see whether the assumption that one or bot h fluids are infinite ly deep is a valid assumption, we have plotted the critical te mperature gradient in the vapor against the wavenumber in Figure 5-3 for different system s. In one of the cases, the fluids are assumed to be infinitely deep whereas in the others both fluid layers have finite depths. Positive temperature gradients mean that the bi layer system is heated from the liquid side, namely from the top. We know from the classical Rayleigh-Taylo r problem that the system is unstable for wavenumbers smaller than the critical wavenumber, however, when there is phasechange taking place at the liquid-vapor inte rface, the interface can be stabilized by heating from the liquid side ev en at small wavenumbers. In the plots, the region under the critical curve is the unstable region. As we increase the liqui d and vapor depths, the critical values of the vapor temperature grad ients approach the values for the infinitely deep layers. We observe that for large wavenumbers the curves overlap whereas for small wavenumbers the critical temperature grad ients in the vapor at the interface diverge from the values for the infinitely deep layers. From Figure 5-3 we can tell that the assumption of infinitely deep layers is al so not a good one for all wavenumbers. For a water-water vapor system, the finite layers, where lgdd0.1 m. deep, give the same results as the infinitely deep layers at al l wavenumbers, however, a 0.1-meter depth for a vapor layer is a relatively large flui d depth under many practical circumstances.

PAGE 98

85 0.5 1 1.5 2 2.5 3 3.5 x 10-7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless wavenumber, Dimensionless critical vapor temperature gradientc b a Figure 5-3. Comparison of models a) Infinitely deep fluid layers b) dl = 0.05 m and dg = 0.05 m c) dl = 0.01 m and dg = 0.01 m In many applications of the Rayleigh-T aylor instability, the vapor layer is shallower than the liquid layer. For th at reason, in Figure 5-4 we plotted the dimensionless critical temperature gradient against the wavenumber for two different cases to see the effect of th e vapor layer depth on the inst ability. Observe that when realistic fluid depth ratios are considered the effect of the fluid de pths is considerably more important. The larger the ratio of the liquid depth to the vapor depth is, the easier it is to stabilize the bilayer system by heating from the liquid side. That is, the minimum temperature gradient required to stabilize th e system is less. For shallow vapor layers, the flow in the vapor is impeded by the bottom wall and the flow in the liquid starts to play a role on the stability. The liqui d flow brings warmer liquid to the troughs encouraging evaporation, thus preventing them from growing deeper. More thorough explanations of the stabilization mechanism through the flow generated by phase-change are given in chapter 3. Thus far, in all the cases we considered, the instability was delayed by heating from the liquid side. Again, when we turn to industrial applications,

PAGE 99

86 we recognize that the system is mostly heat ed from the vapor side, namely heated from below. Figure 5-5 gives a plot of the critical temperature di fference across the bilayer system against the wavenumber where both of the fluids have finite depths and the system is heated from below, namely the va por side. The liquid and the vapor layers are taken to be 1mm and 0.1mm deep, respectively. The region above the curve is the stable region. So, this calculation also pertains to a physically realizable experiment. Observe that the curve shows an increase of the cri tical temperature differe nce with wavenumber at low wavenumbers, a maximum and then a de crease. This shows the dual effect of wavenumber. The waviness of the disturban ce finds itself participating in transverse diffusion of heat and momentum as well as in stabilization on account of surface tension. At low wave numbers, surface tension plays a minimal role. Transverse diffusion itself plays two roles. In the domain it offers a dissipative effect and therefore plays a stabilizing role but at the interface it encour ages the movement of the interface as hot locations move toward cold regions and so forth. This is the destabilizing effect. At high wave numbers the stabilizing effect of su rface tension and the critical temperature difference for the realization of stability is decreased.

PAGE 100

87 0.5 1 1.5 2 2.5 3 3.5 x 10-7 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Dimensionless wavenumber, Dimensionless critical vapor temperature gradienta b Figure 5-4. Effect of shallow fluid layers. a) dl = 50e-3m and dg = 50e-3m b) dl = 50e-3m and dg = 5e-3m Note that heating from below can stabili ze an arrangement where a vapor underlies its own liquid at a wavenumber at which the classical Rayleigh-Taylor analysis indicates otherwise. In summary, we conclude in this chapter th at the assumption of infinitely deep fluid layers give results that are very different than for the fluid laye rs with finite depths as one can conclude that a system that is unstable for a given temperature gradient in the vapor based on the infinitely deep layers assumption may actually give stability for fluid layers with finite depths.

PAGE 101

88 0.5 1 1.5 2 2.5 3 3.5 x 10-7 -1 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 Dimensionless wavenumber, T (C) Figure 5-5. Critical temp erature difference vs. dime nsionless wavenumber Another result of this chapter is that th e use of a phenomenological coefficient is questionable given the observation fact that su ch a coefficient itself is a function of the wavenumber.

PAGE 102

89 CHAPTER 6 EFFECT OF VAPOR DYNAMICS ON TH E RAYLEIGH-MARANGONI PROBLEM Introduction In keeping with the general goals of this thesis as stated in Chapter 1 we move on to the discussion of the role of vapor dynamics in Raylei gh Marangoni convection without phase change. In this regard we fi rst discuss the physics of the pure RayleighMarangoni problem. The physical argumen ts made are then supported by the calculations done using nothing more than a li near stability analysis. The results of calculations are followed by a discussion of experiments that were done to demonstrate the calculations. When we have a liquid layer underlying a ga s layer, say air, a nd heat the bilayer from below, there are several mechanisms th rough which the convection will be initiated. To better understand the possibl e convection mechanisms that can occur, let us consider Figure 6-1. a b c Figure 6-1. Different type s of convection mechanisms

PAGE 103

90 When the bottom layer, i.e., the liquid layer, is deep compared to the top layer (case a), it is more likely for the convection to start in the liquid layer due to buoyancy. The convection in the liquid layer will cause the top layer to convect through viscous coupling at the interface. For ce rtain depths in both fluid layers (case b), the interface will be an isotherm and the bilayer system will act as two separate one-layer fluid systems. The convection in each layer will be due to buoyanc y and we say that the layers are coupled thermally. Another mechanism, one that is of importance to us, is when the gas layer is deep compared to the liquid layer (case c). In the latter case, conv ection will start in the gas layer due to buoyancy. In this case, the convection in the air cannot cause the liquid to convect through viscous coupling because of the low ratio of the gas viscosity to the liquid viscosity. However, the convection in the upper layer will cau se perturbations in the temperature at the interf ace, which in turn will cause Marangoni convection along the interface and consequently in the liquid layer. This last case is important to us because it demonstrates the role of an act ive vapor layer on the onset of convective instability. Note that if the gas layer were assu med to be passive, we would not expect to predict the sort of convection depicted in case c. Several experiments have been done as part of this doctoral st udy and the results of these experiments have been compared to the calculations to verify the theory. In the following section, we will briefly describe the experimental setup. This setup was originally designed and constructed by Johns on (1997). Fine details on the experiment such as the data acquisition system and cont rol algorithms are given in the thesis of Johnson (1997).

PAGE 104

91 Experimental Apparatus and Procedure The objectives of the experiments were to study the effect of the air height on the flow pattern of the lower liquid. The e xperimental apparatus was designed so that various air heights, as well as various liquid aspect ratios co uld be used. A schematic of the test section is given in Figure 6-2 and th e picture of the experimental setup is in Figure 6-3. The lower heating element consisted of a heating plate under a hollowed cylinder made of lucite, ” thick, called the lower bath. The top and bottom of the cylinder were capped with thin copper plates, ” thick, and the interior of the cy linder was filled with water. A magnetic stir bar was placed in the water and the entire lower bath sat on top of a magnetic stirrer. The stirred water helped to stabilize any temper ature fluctuations from the heating block, ensuring a uniform temperature up to 0.05oC. Figure 6-2. Schematic of the experimental setup The test section itself consis ted of four pieces: a liquid in sert, an air insert, a zinc selenide IR transparent window and a clamp for st ructural integrity. The oil insert, the air

PAGE 105

92 insert and the clamp are made of lucite. To en sure that the liquid-gas interface was flat, a pinning edge was used in the liquid insert. Th e liquid insert could be made of different radii and heights, to achieve the desired liquid aspect ratio. The radius of the air insert was always the same as the radius of the li quid insert. However, different air heights could be used for the same liquid insert. The lucite clamp fastened the liquid and air inserts down, preventing silicone oil from leak ing and ensuring structural integrity. The outer radius of the oil inserts was 1 13 16”. The outer radius of th e air inserts were 2 ”. However to reiterate, all of the inserts had the same inner radius. The zinc selenide window had a diameter of 2” and was 5 mm thick. The accuracy of each piece was machined to within 0.1 mm. An Inframetrics, model 760, infrared camera was used to visualize the flow patterns. This particular mode l of infrared camera is capabl e of measuring in the 3 to 5 m range or the 8 to 12 m range. However, only the 8 to 12 m range was used. The infrared camera was placed directly above the test section and measured the infrared radiation being emitted by the silicone oil. As silicone oil readily absorbs infrared radiation, only the radiation from the firs t few angstroms of the silicone oil interface could be detected. An effective emissivity could be calibrated, and programmed into the infrared camera to find the temperature of th e interface. The IR camera was never used to measure the temperature rather it was used to detect the variation in interfacial temperature from which the flow structure could be deduced. Zinc selenide is 60% transparent to infrared radiation in the 8 to 12 m range. Additionally an anti-reflective infrared polymer was coated on the zinc sele nide window by II-VI Inc. This coating was

PAGE 106

93 necessary to eliminate any false images generated by reflected, ambient infrared radiation. Figure 6-3. Experimental setup To control the temperature at the top of th e test section, an infrared transparent medium was needed. The medium of choice was air. A lucite box enclosed the test section, the lower bath, the magnetic stirrer, the infrared camera and the heating control elements. A heater, a fan, and a car radiator were used to control the temperature of the ambient air in the lucite box. To increase th e ambient temperature, the heater would turn on. Cool water, pumped through the radiator, continuously removed h eat from the air. Proper mixing was ensured by a blowing fan. The overall temperatur e control across the top of the zinc selenide and the bottom of the top copper plate of the lower bath was never more than 2.4% of the setpoint. When the smoothing of the temperature fluctuations due to the highl y conductive plates (copper and zi nc selenide) are taken into consideration, the actual c ontrol across the two fluid laye rs was probably better. The temperature of the top and bottom plate, and the temperature of the cooling water were

PAGE 107

94 measured. All of the temperatures were r ecorded and controlled using a computerized data acquisition system. The temperatur e control was later improved to 0.05oC, with no substantial effect on the experimental results. The infrared images were sent to a VCR and displayed on a TV. As some of the experiments lasted several days the VCR was wired into the co mputer so that the control program would periodically record and stop record ing the images. Each experiment was conducted in the same systematic manner. First the silicone oil was placed into the liqui d insert. The corre sponding air insert was chosen and placed on top of the liquid inset. Th e lucite clamp was placed on top of the air inset and screwed into the lower bath. The level of the oil-ai r interface could be ch ecked by both a standard bubble level and the infrared camera. As it turn s out, the oil-air interf ace acts as a lens to infrared radiation, concentrati ng infrared radiation if the surface is depressed in the center or diffusing infrared radiation if the interface is elevated in the center. If there is an insufficient amount of silicone o il in the insert, the liquid-air interface will be depressed and the temperature will appear to be higher in the center, even though the temperature is uniform. If there is too much oil in the insert the liquid-air interface will be elevated in the center and the temperature will appear to be lower in the center. Checking the oil level with the infrared camera was more accurate than the bubble level. Once the test section was secured, a temp erature difference was applied across the liquid-gas bilayer. The initial temperat ure difference was less than the critical temperature difference necessary to initiate convection in either layer. When a temperature difference was applied, it was held constant for several time constants. The longest time constant in these experiment s is the horizontal thermal diffusion time

PAGE 108

95 constant. Here the thermal diffusivity of silicone oil is 2 71.110 sec m and the typical diameter is about 25 mm to give a horizontal time constant of 2 3 2 2 72510 d 1.6 1.110 sec m hrs m The temperature difference for these experiments was held constant for four hours. However, steady st ate was usually reached well within the fourhour period. After the temperature differen ce was held constant for four hours, the temperature difference was increased to a ne w set point and held constant for another four hours. This procedure was repeated unt il the fluid began to flow. Once the fluid began to flow, the temperature differe nce and the flow pattern were noted. Experimental Results We have compared the experimental resu lts with two types of calculations. The results of the linear calculati ons are in good qualitative agreem ent with the experiments. Calculations were done to support the physical arguments given above using the physical properties in Table 6-1. The first set of calculatio ns was done using a code written by the author assuming that the two fluid layers ar e bounded vertically by two rigid plates of infinite lateral extent. The physical model and the equations are given in Chapter 1 and 2. Let us start out by verifying the physical arguments we made us ing the first set of calculations. The fluid properties used in the calculations are those that corre spond to fluids used in the experiments, viz. silicone oil and ai r. The calculations were done assuming the liquid and air depths to be 5 mm and 3mm, respectively.

PAGE 109

96 Table 6-1. Physical properties of silicone oil and air at 20 C 3kg 968 m 3kg 1.2 m 3kg 96.810 m sec *5kg 1.910 m sec 1J 1.610 m sec Ck *2J 2.610 m sec Ck 2 7m 1.110 sec 2 *5m 2.110 sec 41 9.610 C *41 32.410 C 3N 20.910 m 4 TN 0.510 m In Figure 6-4, z0 represents the liquid-gas interface. It is evident that at the onset of convection, for this is what the calculations predict, at any point in the xdirection, the z-components of velocity in both fluids have opposite signs, that is, the convection rolls in each phase are in the opposite direction of the other phase. If at any given ‘x’ the liquid moves upward, the air layer at that point in the direction of ‘x’ is moves downward. The flow in the liquid layer is st ronger than the flow in the gas layer as convection initially starts due to buoyancy in the liqui d layer. The calculations therefore represent the scenar io depicted by case (a). Next, the calculations are de picted in Figure 6.5 for a si licone oil-air system where the liquid and air depths are 5 mm and 5mm, respectively. We will see that the calculations correspond to a case si milar to case (b). We obser ve that there is strong flow in both fluid layers due to buoyancy. There is a small roll in the va por near the interface, which arises because of the no-slip condition at the interface. To be tter visualize the flow pattern in the system, let us turn to Figure 6-6. In Figure 6-6, z0 represents the interface between silicone oil and air. It is more obvious from this plot, that the

PAGE 110

97 convection in the gas near th e interface is very weak in comparison. The convection occurs in both domains due to buoyancy. A further adjustment in the gas depth could result in a more accurate representation of case (b) than the above calculations. -5 0 x 10-5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 z-component of velocityDimensionless depth, z Figure 6-4. z-component of velocity in a silic one oil-air system at th e onset of instability where the silicone oil and air depths are 5mm and 3mm, respectively. -2 0 2 4 6 x 10-5 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 z-component of velocityDimensionless depth, z Figure 6-5. z-component of velocity in a silic one oil-air system at th e onset of instability where the silicone oil and air depths are 5mm and 5mm, respectively.

PAGE 111

98 Now let us consider the case c where the air layer is deeper than the silicone oil layer. The convection starts initially due to buoyancy and the flow in the vapor causes convection in the liquid through surface te nsion gradients at the interface. 0 1 2 3 4 5 6 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 Dimensionless width, xDimensionless depth, z Figure 6-6. Flow profiles in a silicone oil-air system at th e onset of convection where the silicone oil and air depths ar e 5mm and 3mm, respectively. Figure 6-7 depicts a silicone oil-air bilaye r system where the silicone oil and air depths are 5mm and 9mm, respectively. In Fi gure 6-7, the convecti on in the air layer is much stronger than the convection in the liqui d layer since the convect ion in the liquid is actually caused by surface tension gradients at the interface, and Marangoni forces are not strong enough to cause flow as strong as the buoyancy forces that initiated the convection in the vapor. The second set of calculations was done using a code supplied by Dauby. In the second set of calculations, a single liquid layer of silicone oil superposed by a conductive, active gas in a closed cylindrical container is modeled. The inte rface between the liquid and the gas is assumed to be flat. All side s of the cylindrical container have no-slip

PAGE 112

99 boundaries. The lower and upper vertical boundari es are rigid, no-slip conducting plates. The lateral sidewalls could be either perfect conductors or perfect insulators. -1 0 1 2 3 4 5 6 x 10-4 -1 -0.5 0 0.5 1 1.5 z-component of velocityDimensionless depth, z Figure 6-7. z-component of velocity in a silic one oil-air system at th e onset of instability where the silicone oil and air depths are 5mm and 9mm, respectively. The calculations are able to predict the cr itical temperature difference (or Rayleigh number or Marangoni number) for different azimuthal modes over a range of the cylinder’s aspect ratio. Figure 6-8 shows the predicted criti cal temperature difference for a range of the air heights for a fixed liquid depth and radius from the calculations as well as the results of experiments performed at se veral air heights. In Figure 6-8 there are three curves for each air height Each of these curves stands for a different azimuthal flow pattern. To better explai n the convection patterns of th ose modes, let us turn to Figure 6-9. m0 is an axisymmetric flow pattern, meaning that the flow will go up in the middle and go down on the sides. In Pict ure 4, one can observe the lighter color indicating the upflow in the middle. m1 on the other hand, is symmetric around an imaginary diameter cutting the cylinder in th e middle. That is, the flow will go up on one side and go down the other, which is very clear in Picture 3.

PAGE 113

100 2 4 6 8 10 12 14 16 18 20 0 1 2 3 4 5 6 7 8 9 10 Air depth, (mm)Critical temperature difference, ( C)1 2 3 4 m=0 m=1 m=2 Figure 6-8. The critical temp erature difference vs. air height for a bilayer system where the silicone oil depth is 8mm a nd the aspect ratio is 2.5. The results of experiments are marked with an ‘x’ on the same plot, and some of the results are also numbered. The pictures of the experiment at the given numbered points are below. 1 2 3 4 Figure 6-9. Pictures of experiments taken by IR camera at the onset of instability We observe that the experiments and cal culations agree perfectly for small and large air depths whereas for intermediate ai r heights the onset of convection is observed before the expected critical temperature diffe rence is reached. When the air depths are very small and very large, the convec tion is due to buoyancy and one can show mathematically that as the Rayleigh convection problem is self-adjoint, the instability is supercritical. For intermediate air depths, Marangoni convection star ts to play a role on the stability of the problem, and Marangoni problem can be s ubcritical for m=0. And in

PAGE 114

101 this case, it turns out to be subcritical. Th at is why we believe that we can observe the convection set in before the critical point is reached. In summary, we conclude that even in the absence of phase-change the fluid flow in the gas layer plays an activ e role in determining the stability of a bilayer system indicating it is a mistake to leave the fluid dynamics of gas layers in convective instability problems. Our conclusions have been backed by experiments. Moreover we have been able to interpret the results of th ese experiments in light of bifurcation theory.

PAGE 115

102 CHAPTER 7 CONCLUDING REMARKS AND FUTURE SCOPE In this chapter, the main discoveries of each chapter are re-evaluated and the future scope for further research is suggested based on the conclusions of this doctoral work. In order to do this, the discussion is divided into three parts. Active Vapor Layers A theoretical model was written to describe the physics of the in terfacial instability arising from evaporation. The numerical cal culations arising form this model showed without any doubt that it is inappropriate to assume that the vapor layer is only a conductive heat transfer medium In fact, it was found that the fluid dynamics of the vapor layer has an important role on stab ilizing the evaporation problem. As past workers have assumed that the vapor can be taken to be infin itely deep and fluid mechanically passive, this assumption has b een shown to not always be true. In a separate experimental study using non-evaporativ e liquids in the presence of air, it was again shown that gas layers are fluid mechani cally active, thus have as strong effect on the onset point as well as the pattern at the onset of instability. These experiments agreed well with companion numerical models. As a result of this work it is suggested that one should numerically investigate the role of vapor layers on the onset flow patterns during phase-change in cylindrical containers. The companion experiments might, however, prove to be harder than corresponding problem in the absence of phase-change. One of the main reasons for the difficulty is the fact that a temperature gradient has to be applied

PAGE 116

103 across the bilayer system and the vapor will c ondense on the cold plate obscuring the IRimaging of the liquid-vapor interface. Th ere are several ways to get around this hindrance. One of them is to heat the syst em from the vapor side. The calculations suggest that the interface will become unstabl e even when heated from above within the range of temperature differences that are suit able for laboratory expe riments, i.e., at low values of temperature gradients. Another way is to not apply a temperature difference across the system and to let the temperatur e gradient occur through evaporation at the interface. This can be done by lowering the va por pressure at the top plate and allowing a small phase-change rate at the interface, which in turn will lower the temperature of the interface. In that case, the liquid will still be po tentially unstable for Rayleigh convection whereas the vapor layer will be Rayleigh stable By using shallow liquid depths, one can delay the buoyancy-driven instability in the li quid and investigate the role of active vapor layer, which can then only convect due to su rface tension gradients or the phase-change. Another necessary calculation to do would be the nonlinear analysis of the problem in its complete form, i.e., including the buoyancyand surface tension gradient-driven convection. Several patterns are expected to emerge with the onset of instability and nonlinear calculations will provide informati on on the stability of these possible patterns and tell us which one of them will prevail after the onset of instability. In the last problem investigated in this wo rk, Chapter 5, the role of an active vapor layer was investigated in relation to a film-bo iling configuration. There is plenty of data in that line of research and the results of th is work will definitely shed more light on the Rayleigh-Taylor instability problem in the presence of phase-change including a finite base state evaporation rate.

PAGE 117

104 Conditions at the Interface of an Already Evaporating Liquid In this doctoral work, all calculations were done assuming a zero phase-change rate in the base state. Natura lly, the thermodynamic equilibr ium assumption at the liquidvapor interface was sufficient to complete th e model for the interf ace conditions. On the other hand, in practical systems this is not always applicable as very high evaporation rates are common in film boiling. In that case, continuum thermostatics is not adequate to describe what is going on at the interface. One can investigate the possibility of a nonequilibrium relation or an accommodation co efficient and a range of possibilities where that coefficient may replace the interf ace conditions used in the theoretical model specified in this study. There is, of cour se, the complication that such an accommodation coefficient would depend upon the wave number. Coupling of Evaporative Instabilities with Other Phase-Change Phenomenon There are plenty of examples where the evaporation at a liquid-vapor interface is accompanied by the solidification at a liquid-soli d interface, such as in coating processes and in the drying of car paints. Evaporative instability and the dept h of the vapor layer can be used to stabilize an otherwise unsta ble solidification front. For that reason, it would be interesting to inve stigate the coupling of evapor ative instability with the solidification problem. Experiments can be performed to investigate the effect of evaporating liquid layers a nd active vapor layers on the morphology of a solidification front. In other words, there is plenty of scope to take the results of this thesis toward other problems with e qually interesting physics and applications.

PAGE 118

105 APPENDIX A THERMODYNAMIC EQUILIBRIUM RELATION AT THE INTERFACE Consider an interface between an evaporating liquid and its vapor. For two phases of a pure species coexisting at equilibrium, LGˆˆ GG where Lˆ G and Gˆ G are the liquid and gas specific Gibbs free energies, respectively. This relation holds regardless of the shap e of the interface. We will make use of this relation to get information about the effect of curvature on the pressure. So, we write LG curvedcurvedˆˆ GG and LG flatflatˆˆ GG Combining all these equations we can write LGLG curvedcurvedflatflatˆˆˆˆ GGGG0 or LLGG curvedflatcurvedflatˆˆˆˆ GGGG0 The changes in ˆ G due to changes in T and P are evaluated using ˆˆ ˆ dGSdTVdP where ˆ S and 1 ˆ V are the specific entropy and volume, respectively. Assuming that the specific entropy and vol ume do not change much on account of small changes in T and P, we get

PAGE 119

106 0 dP V dT S dP V dT SG G G L L L or 0 dP V dP V dT S SG G L L L G The above equation can be rea rranged and rewritten assuming GL GL11 dPdP to GLGG G1 SSdTdP0 Assuming an ideal gas, we then have 0 dP P RT dP V dT S SG G L L L G This becomes 0 dP P RT dP V dT T H HG G L L L G But, HGHL is equal to the latent heat of evaporation, Therefore, 0 dP P RT dP V dT TG G L L or 0 dP P R dP T V dT TG G L L 2 We can integrate both sides to get a rela tion between the temperature and pressure at the interface. The integrat ion limits are between the flat and curved states. Hence 0 dP P R dP T V dT TG C G F L C L F C FP P G G P P L L T T 2

PAGE 120

107 This equation can be worked on further to get the first order terms, namely the perturbed equation. It is 2 1 2 0 L G 0 1 G 1 G 0 0 1 1 2 0x Z T V dz dP Z P P R dz dT Z T T

PAGE 121

108 APPENDIX B A SIMPLE DERIVATION TO OBTAIN TH E CRITICAL CONDITION FOR ONSET OF CONVECTION IN A FLUID LAYER Consider an infinite horizont al layer of fluid of a finite depth, d, where a constant temperature gradient, 0 0dT T dzd is maintained by heating from below. Here, hotcoldTTT Making use of the continuity e quation, the dimensionless perturbed Navier-Stokes equations are br ought into the following form: 22 222 z11 22dd1 vRa T dzdzPr B.1 2 2 1z1 2d Tv dz B.2When we set 0 to determine the critical Rayleigh number and eliminate T1 between the two equations, we get 3 2 22 z1z1 2d vRa v dz B.3The temperatures at the top and bottom plat es are kept constant Using the domain energy equation we turn the cons tant temperature conditions at 0z0 and 0z1 into 2 2 2 z1 2d v0 dz B.4The fluid is confined between the top and the bottom plates, that is z1v0 at 0z0 and 0z1 Two more boundary conditions are ne eded. For a simpler algebra, let

PAGE 122

109 us require that both surfaces are stre ss free and the boundary conditions at 0z0 and 0z1 are given by 2 z1 2dv 0 dz B.5Using the no mass transfer and stress fr ee boundary conditions it is easy to convert the boundary condition obtained from the temperature boundary conditions into 4 z1 4 0dv 0 dz B.6From equations A.3 through A.6, we determine that 6 z1 6 0dv 0 dz By differentiating the equation A.3, we find that 8 z1 8 0dv 0 dz By further differentiations of A.3, we conclude that all even derivatives of z1v vanish on the boundaries. Consequently, 2m z1 2m 0dv 0 dz at 0z0 and 0z1 for m = 1, 2, … At this point, it is not hard to guess that the solution is in the form z10vA sin nz (n = 1, 2, …) where A is an unknown constant, and n is an in teger. After substituting this solution into A.3, we get the following relation: 3 222 2n Ra For any given positive value of 2 the lowest value of Ra occurs when n = 1, accordingly, 3 22 2Ra

PAGE 123

110 Our aim is to determine the critical wa venumber that corresponds to the lowest value of Ra. In order to determine th e critical Rayleigh number along with its corresponding wavenumber, let us take the first derivativ e of Ra with respect to 2 and set it equal to zero. It is 23 2222 224d Ra 30 d from which we find that 2 2 c2 The corresponding value of Ra is 4 c27 Ra657.5 4 Below is the plot of the Raylei gh number against the wavenumber. 0 1 2 3 4 5 6 7 8 0 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Dimensionless wavenumber, Ra Figure B–1. Critical Rayleigh num ber vs. dimensionless wavenumber

PAGE 124

111 APPENDIX C A SIMPLE DERIVATION FOR THE RAYLEIGH TAYLOR INSTABILITY When a vapor underlies a liquid, the resulti ng configuration is pot entially unstable. This problem, known as the Rayleigh-Taylor problem, is well known and a nice account of this instability is given by Chandrasekhar (1961). Here I will give a derivation of the problem for infinitely deep fluid layers and find the critical wavenum ber at the onset of instability. Linearized stability is performed in th e standard way leading to the following dimensionless equations of motion 2 2 2l z1 2d v0 dz w and 2 2 2g z1 2d v0 dz w When we subject the vapo r and liquid domain equati ons to the conditions at 0z and 0z respectively, we get the following expressions for the velocity and temperature fields. 00zz lll z10vA B z wwee 00zz ggg z10vA B z wwee At the interface, 0z0 the mass balance turns into

PAGE 125

112 z1v0 and z1v0 The no-slip conditions becomes z1z1 00dvdv dzdz The tangential and normal stress balances assume the following dimensionless form: 2g2l 2g2l z1z1 z1z1 22 00dvdv v v0 dzdzww and lg2 3gg3ll 2222 z1z1z1z1 1 33 0000d dvdvdvdv 33Z0 dzdzdzdz g wwww where is the ratio of the vapor visc osity to the liquid viscosity. Using the interfacial conditions we arrive at 1 2 lg cdg w

PAGE 126

113 APPENDIX D THE BOUSSINESQ APPROXIMATION The Rayleigh problem is a stability probl em because of the nonlinearity introduced to the problem through the convective heat eq uation term in the domain energy equation. After we linearize the domain equations we obs erve that the temperature and velocity are still dependent on each other through the energy equation. However, that relation is not the only reason for the Rayleigh convection. The assumption that th e fluid density is a linear function of temperature makes the fl uid dynamics problem dependent on the heat transfer problem. Clearly, the analysis of the Rayleigh problem may be very complicated. For that reason, the assumption of Boussinesq approximation is employed. It is based on the assumption that 0TT1 T1 T which is basically 1 where is the change in the fluid density for the temperature difference applied. The density can be expressed by the first te rm of a Taylor series expans ion about the base state as 000 TTTT T Let us now look at the perturbed equati on of motion and the energy equation, and perform a dimensional analysis. The perturbed energy equation is 2 1 10001 0T vT T t and we learn from a dimensional analysis that the prevailing velocities would be of order d, then plugging that velocity scale into the motion equation

PAGE 127

114 2 1 01101 0v P g v t we determine that the term where the perturbed density multiplies gravitational acceleration has to be of order 10-4 or 10-5 for most fluids for a depth of 1 cm. 0TT1 T itself is of the order of 10-4 or 10-5 before it even multiplies the gravitational acceleration and a temperature difference. A simple dimensional analysis of the motion equation tells us that the density, can be replaced by its base state value, 0 in all the terms except in the term that multiplies the gravitational acceleration. During the course of my study, Prof. Narayanan and I had an interesting observation on Boussinesq approximation. A Remark on the Boussinesq Approximation Consider a fluid in a rigid container. The energy balan ce for the fluid is given by 2d1 ˆ Uvqv dt2 where PIS S is the viscous stress tensor a nd for Newtonian fluids it is S2D Now let us take the volume integral of th e equation above and rearrange it. We get 2d1 ˆ Uvqvn dt2VVS At steady-state all terms excep t the heat conduction term are equal to zero. The last integral is zero because of the rigid walls and the no slip at the boundary which require vn0 and vt0 So, we have q0V

PAGE 128

115 Now let us look at another form of the energy equation for a fluid, which can be brought into the same form as the one given in the beginning of this Appendix. It is d ˆˆ UUvunqPv2D:D dtVSVVV At steady-state and rigid boundaries, the inte grals on the left hand side are equal to zero. We have shown that q0V and then the equation reduces to Pv2D:DVV The continuity equation in the domain in our analysis is v0 Thus, we are only left with one term and the equation further reduces to 2D:D0V When there is convection in our system, we know that 2D:D0V This poses no problem in the linear stabil ity analysis because this term disappears in the first order perturbation of the equation and is automatically equal to zero, however, it seems that nonlinear analysis there is a violation when one uses Boussinesq approximation.

PAGE 129

116 APPENDIX E THE UNIT NORMAL OF A SURFACE In this Appendix, the unit normal vector in 2D system where the normal vector is on a surface on the z-axis. Le t a surface be denoted by (,)0 fzZxt Then, f is positive on one side of f = 0, negative on the other, and the normal pointing into the region where f is positive is given by n f f Here, f f f ik x z Then, the equation for the normal stress balance is given by 1 2 2n 1Z ik x Z x The derivation of the unit tangent vector is straightforward using the definition nt0 from which we get 1 2 2t 1Z ik x Z x

PAGE 130

117 APPENDIX F DERIVATION OF THE SURFACE SPEED Let a surface be denoted by (,)0 frt where (,,) rrxyz and the unit normal of the surface is given by n f f Let the surface move a small distance s along its normal in time t. Then, ( n,) f rstt is given by (,) ( n,)(,) n(,) ... tfrt frsttfrtsfrtt whence ( n,)0(,) f rsttfrt requires (,) n(,) t f rt sfrtt The normal speed of the surf ace, u, is then given by (,) u n(,) f rt s t tfrt Now, using the definition of the unit normal given earlier we get u f t f In our problems, the definition of u becomes

PAGE 131

118 1 2 2u 1 Z t Z x

PAGE 132

119 APPENDIX G THE INTERFACIAL ENERGY BALANCE There are three difference equations of intere st to us in thermal problems and they are n u v n u v* n T u v v t n H 2 n T u v v * s 2 ***2****1 ˆ vvuqTvn2H nu 2 1 ˆ vvuqTvn 2 U U where uu n Notice that n points out of the non-asterisk phase into asterisk phase and 2H is taken to be negative at a crest, looking along n The first of these is form invariant under a change of observer; the second and third are not. The second and third equations are taken to hold in a laboratory frame and our aim is to use the second equation to rewrite the mechanical terms in the third equation, thereby producing a form invariant differen ce equation that accounts for the latent heat released as a phase-change front moves into the asterisk fluid. To do this, work only on the left hand-si de, take the dot pr oduct of the second equation 2 with u, subtract the result from the third equation to obtain

PAGE 133

120 2 ******2***1 ˆ vu vuvvuqTvun2H nu 2 1 ˆ vu vuvvuqTvun 2 U U and use TIS P and ˆˆ P HU to get 2 2 2 ***2****11 ˆ vuuvuqSvun 22 11 ˆ vuuvuqSvun 22 H H The resulting difference equation is then n u v S q u v u v 2 1 ˆ n u v S q u v u v 2 1 ˆ* * 2 * 2 H H and observe that it is form invariant.

PAGE 134

121 APPENDIX H THE EXPANSION OF A DOMAIN VARIA BLE AND ITS DERIVATIVES ALONG THE MAPPING In this appendix, the pert urbation equations used in the theoretical work are explained. This is done for a one-dimensional domain and the reader is referred to Johns and Narayanan (2002) for details Let ‘u’ denote the solution of a problem in an irregular domain D where D is not specified and must be determined as part of the solution. Imagine that D lies in the vicinity of a reference domain 0D and can be expressed in te rms of the reference domain via a parameter Let u be a function of a spatial c oordinate ‘z’. Then ‘u’ must be a function of ‘ ’ directly because it lies on D and also because it is a function of ‘z’. By an irregular domain it is mean t a domain that is not convenient. Now to solve for u and obtain D simultaneously we can solve a seri es of companion problems defined on the nearby regular or convenient domain, 0D, which we called the reference domain. What needs to be done is to discover how to determ ine u in terms of the solutions to problems defined on 0D. The points of 0D will be denoted by the coordinate 0z and those of D by the coordinate z. Imagine a family of domains D growing out of the reference domain 0D. u must be determined on each of these. The point z of the domain D is then determined in terms of the point 0z of the reference domain 0D by the mapping 0zz, f

PAGE 135

122 The equation above is a little more genera l than necessary to explain the mapping for our purposes. In the problem s of this thesis, only one part of the domain is irregular, namely the interface. Now, to get going let us expand the function f in powers of as 2 00 2 00 2z,0z,0 1 z,z,0... 2ff ff where 00z,0zf and the derivatives of f are evaluated holding 0z fixed. Then in terms of the notation 0 10z,0 zz,f 2 0 20 2z,0 zz,f etc. the mapping can be written as 2 010201 zzzz,0zz,0... 2 Also, the boundary of the reference domain must be carried into the boundary of the present domain by the same mapping. Th e function, Z, which describes the boundary of the new domain, inherits its expansion in powers of from the mapping given in and can be written as 2 010201 ZZZZ,0ZZ,0... 2 It is the iZ’s that need to be determined to specify the domain D in terms of the domain 0D.

PAGE 136

123 Accordingly, uz, can be expanded in powers of along the mapping as 2 00 2 0 2duzz,0duzz,0 1 uz,uzz,0... d2d where du d denotes the derivative of th e function u depending on z and taken along the mapping. To obtain a formula for 0duzz,0 d differentiate u along the mapping taking z to depend on holding 0z fixed. Using the chain rule, this gives 0duz,uz,uz,z, dzf Now, set to zero in the above equation to get 000 10duzz,0uz,0uz,0 zz, dz Then, introduce the definition of 1u via 0 10uz,0 uz and observe that 000 0uz,0uz zz to get 000 1010 0duzz,0uz uzzz, dz All the other orders of the derivatives of u can be determined the same way. Finally, if a domain variable needs to be specified at the boundary it is written as

PAGE 137

124 000 1010 0duzZ,0uZ uZZZ, dz When additional derivatives are obtained and plugged into the expansion of u, it becomes 2 22 000 1 0112112 2 0000uuu u 1 uz,uuzu2zzz... z2zzz The careful reader will notice that the mapping 1z does not appear in the domain equations given in the thesis. To show why, let us work out an example and so let u satisfy an equation u 0 z in the new domain. Our reference domain, however, is in the coordinate system, 0z. Thus, 0 0z uu zzz Hence, we must differentiate the right hand side of the expansion of u with respect to 0z, holding fixed and then multiply it by 2 0 12z zz 1 1... zz2z This gives us to the first order 2 00 1 1 2 000uu u u z... zzzz See how the derivatives of 1z and 2z are lost during the step where we multiplied the derivative of u with respect to 0z by 0z z Now going back to our example and plugging this in, we get

PAGE 138

125 2 00 1 1 2 000uu u u z...0 zzzz That is to say that in the domain 0 0u 0 z and accordingly, 2 0 2 0u 0 z which gives 1 0u 0 z The mapping is lost from domain equations If surface variables were considered the mapping would not be lost. What we show in this appendix is mostly sufficient for us to get on with the work of this thesis.

PAGE 139

126 APPENDIX I NUMERICAL METHOD The Chebyshev spectral tau method was us ed to solve the resulting eigenvalue problem. Numerical methods using this method are often considerably faster with greater accuracy than other standard techniques such as finite differencing. The tau method easily incorporates complicated boundary conditions. A great explanation of the Chebyshev spectral tau method is given by Johnson with nice examples. In this section, the Chebyshev spectral tau method is described br iefly, the reader is referred to Canuto et al. and Gottlieb and Orszag. Spectral method is discretiza tion scheme developed from the method of weighted residuals (Finlayson and Scriven). He re follow the details of the method: Consider the following problem: u LuAu t with boundary conditions represented by Bu0 where L and A are linear operators, B is a linear boundary operator, and is the eigenvalue. The solution vector u can be expressed in terms of infinite series of trial functions. The trial functions chosen are the Chebyshev polynomials, then u is written as 0u,ttT()nn n x ax

PAGE 140

127 The solution vector u is then approximated by the truncated series up to some finite value of N and is given by N N 0u,tu,ttT()nn n x xax Now, Nu,t x must explicitly satisfy the boundary conditions. Since this is not an easy task, k more terms are added to the polynomial, so the expression for Nu,t x becomes Nk N 0u,ttT()nn n x ax Here, N is the dimension in the dom ain and k is the number of boundary conditions. Now there are N+k unknowns and N equations. The k additional equations obtained from the k boundary equations. Nk 0tB T()0nn nax The Chebyshev polynomials are orthogonal on th e interval (-1,1) with respect to the weight function 1 2 2()1 wxx These polynomials satisfy the following: 1 10 T(x),T()()T(x)T() 0 1 0 2mnmnmn xwxxdxmn mn Now choose a test function T(x)m and take the inner pro duct with the problem to be solved. It gives NNNNNN 000000d(t) T(),T()(t)T(),L T()(t)T(),A T() dtn mnnmnnmn mnmnmna x xaxxaxx

PAGE 141

128 which reduces to NNNN 0000d(t) (t)T(),L T()(t)T(),A T() dtn nmnnmn mnmna axxaxx The equation above can be solved know using the definition of Chebyshev polynomials.

PAGE 142

129 LIST OF REFERENCES Adham-Khodaparast, K., Kawaji, M. & Antar, B.N., The Rayleigh-Taylor and KelvinHelmholtz stability of a viscous liquid-vapor interface with heat and mass transfer, Phys. Fluids, 7 359-364 (1995). Badratinova, L.G., Colinet, P., Hennenberg, M. & Legros, J.C., On Rayleigh-Taylor instability in heated li quid-gas and liquid-vapor systems, Russ. J. Eng. Thermophys., 6 1-31 (1996). Bnard, H., Les tourbillons Cellulaires da ns une Nappe Liquide, Rev. Gen. Sciences Pure Appl., 11 1261 (1900). Block, M., Surface tension as the cause of Bnard cells and surface deformation in a liquid film, Nature, 178 650 (1956). Chandrasekhar, S., Hydrodynamic and Hydromagnetic Stability, ch. 10 (Oxford University Press, Oxford, 1961). Grindrod, P., The Theory and Applications of Reaction-Diffusion Equations Patterns and Waves, Second Edition, (Oxford University Press, Oxford, 1996). Ho, S.P., Linear Rayleigh-Taylor stability of viscous fluids with ma ss and heat transfer, J. Fluid Mech., 101 111-128 (1980). Hsieh, D.Y., Effects of heat and mass tran sfer on Rayleigh-Taylor instability, Transactions of the ASME, 94 156-162 (1972). Hsieh, D.Y., Interfacial st ability with mass and heat transfer, Phys. Fluids, 21 745-748 (1978). Huang, A., Joseph, D.D., Instability of th e equilibrium of a liquid below its vapour between horizontal heated plates, Journal of Fluid Mechanics, 242 235 (1992). Johns, L.E., Narayanan, R., Interfacial Inst ability, (Springer-Verlag, New York, 2002). Johnson, D., Chebyshev polynomials in the Sp ectral Tau method and applications to eigenvalue problems, NASA Contractor Report 198451 (1996). Johnson, D., Narayanan, R., A tutorial on th e Rayleigh-Marangoni-B enard problem with multiple layers and side wall effects, Chaos, 9 124 (1999).

PAGE 143

130 Margerit, J., Colinet, P., Lebon, G., Iorio, C.S ., Legros, J.C., Inter facial nonequilibrium and Benard-Marangoni instabil ity of a liquid-vapor system, Physical Review E, 68 041601 (2003). Mihaljan, J.M., A rigorous exposition of the Boussinesq approximations applicable to a thin layer of fluid, Astrophysical Journal, 136(3) 1126-1133 (1962). Nield, D.A., Surface tension and buoyancy e ffects in cellular c onvection, J. Fluid Mech., 19 1635 (1964). Oron, A., Nonlinear dynamics of thin evaporat ing liquid films subject to internal heat generation, in Fluid Dynamics at Interfaces edited by W. Shyy and R. Narayanan (Cambridge University Press, Cambridge, 1999), 3-15. Oron, A., Davis, S.H., Bankoff, S.G., Longscale evolution of thin liquid films, Reviews of Modern Physics, 69 931 (1997). Ozen, O., Johnson, D. & Narayanan, R., Obs ervations on interfacial convection in multiple layers without a nd with evaporation, in Interfacial Fluid Dynamics and Transport Processes Lecture Notes in Physics Vol. 628, edited by R. Narayanan and D. Schwabe (Springer-V erlag, Berlin, 2003), 59-77. Pearson, J.R.A., On convection cells indu ced by surface tension, J. Fluid Mech., 4 489-500 (1958). Rayleigh, Lord, On convection currents in a ho rizontal layer of fluid, when the higher temperature is on the under side, Phil. Mag., 32(6) 529 (1916) Rednikov, A., Colinet, P., Velarde, M.G., and Legros, J.C., Rayleigh-Marangoni oscillatory instability in a horizontal liquid layer h eated from above: Coupling between internal and surface waves, Phys. Rev. E, 57 2872 (1998). Rubinstein BY, Bankoff SG, Davis SH, I nstability of subcooled boiling film on a vertical wall, International J. of Heat and Mass Transfer, 45 (25) 4937-4948 (2002). Sarma, G.S.R., Interaction of surface-tens ion and buoyancy mechanisms in horizontal liquid layers, J. Ther mophys. Heat Transfer, 1 129 (1987). Scriven, L.E & Sternling, C.V., Interfacial turbulence, hydrodynamic instability and the Marangoni effect, AIChE J., 5(4) 514 (1959) Ward, C.A., Stanga, D., Inter facial conditions during evap oration or condensation of water, Physical Review E, 64 051509 (2001).

PAGE 144

131 BIOGRAPHICAL SKETCH Ozgur Ozen was born in Turkey. He gr aduated from Bogazici University in Istanbul, Turkey, receiving a B.S. degree in Chemical Engineering in 1999. He then attended the University of Florida for gra duate studies under the supervision of Prof. Ranga Narayanan. In 2004, he graduated from the University of Florida with a Ph.D. in chemical engineering.


Permanent Link: http://ufdc.ufl.edu/UFE0006000/00001

Material Information

Title: Effect of Vapor Dynamics on Interfacial Instabilities
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0006000:00001

Permanent Link: http://ufdc.ufl.edu/UFE0006000/00001

Material Information

Title: Effect of Vapor Dynamics on Interfacial Instabilities
Physical Description: Mixed Material
Copyright Date: 2008

Record Information

Source Institution: University of Florida
Holding Location: University of Florida
Rights Management: All rights reserved by the source institution and holding location.
System ID: UFE0006000:00001


This item has the following downloads:


Full Text












EFFECT OF VAPOR DYNAMICS ON INTERFACIAL INSTABILITIES


By

OZGUR OZEN

















A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


2004


































Copyright 2004

by

Ozgur Ozen

































This doctoral thesis is dedicated to the teachers all around the world who dedicated their
lives to future generations.















ACKNOWLEDGMENTS

First of all, I would like to thank Professor Ranga Narayanan for his support and

advice. He has been both a mentor and a friend. His enthusiasm for his work and

dedication to teaching have inspired me to accept every challenge.

Many thanks and love go to my wife, Jen, for her constant support and generous

love throughout my graduate education. I also have the highest appreciation for my

parents and my sister for their love and nurturing throughout my educational career.

Despite the large distance, they have always been there for me and have always

encouraged me to overcome every obstacle.

The members of my PhD committee, Dr. Jason Butler, Prof. Oscar Crisalle and

Prof. Loc Vu-Quoc, also deserve my gratitude. I have really enjoyed taking classes from

Prof Vu-Quoc and being a teaching assistant for Prof. Crisalle. Their teaching

philosophies have deeply influenced me.

Many thanks go to my friends Mehmet Ali Elicin, Eric Theisen, Kerem Uguz and

Berk Usta for always being there when I needed help and friendship.

Special thanks go to Dr. Charles Martucci, whose advice changed the direction of

my life in order to pursue a career in academics. I cannot repay him for his mentoring

and support.

My final acknowledgments are to my parents-in-law, Betty and Dave Hackworth,

who have looked out for Jen and me in every possible way. Without their help this thesis

would have taken longer to write.
















TABLE OF CONTENTS

page

A CKN OW LED GM EN TS ....................................................... 4

LIST OF TABLES .............. .... .... .... .................... viii

LIST OF FIGURES ...................................... ix

ABSTRACT.................................... xi

CHAPTER

1 INTRODUCTION ................... ...................................... ......... .......

Physics of Various Instabilities Connected to Vapor-Liquid Phase Change................
Evaporative Instability......................... ............. ... ....1
Rayleigh Convection or Buoyancy-Driven Convection........................6
Pure Marangoni Convection or Interfacial Tension Gradient-Driven
C onvection ................................... .................. ...... ...8
R ayleigh-T aylor Instability ....................................................... 12
Plan of the Thesis....................................... .............................. .........13

2 THE M ATHEM ATICAL M ODEL.................................................................... 15

Nonlinear Equations ..................... ................. 16
D im ensionless N onlinear Equations..................................................................... 19
Som e Observations .............. ................ .. ........ .......................21
The Linear M odel ................. ..............................22
Perturbation Equations .......................................................23
The Base State Solution and the Perturbed Equations ......................................24
Endnote: Challenges in phase-change problems: Interface conditions ....................28

3 INSTABILITY DUE TO EVAPORATION THE LINEAR MODEL....................30

The Pure Phase-Change Problem ....................................................... 30
Stabilizing Effect of the Vapor Flow................................... ....... 38
Endnotes ..................................... ........... .. ... ... ....... .......... 42
The Exchange of Instability when Heated from the Vapor Side........................42
The Assumption of Flat Surface in the Stability Analysis of the
Phase-Change Problem ................. ...............................42


v









4 WEAKLY NONLINEAR ANALYSIS ........................................................ 44

The Pure Phase-Change Problem ........................................................ 44
Perturbation Expansions ........................ .................47
Overview ................................................48
Perturbed Equations.................................... ............... 50
R results of C calculations ............................................................57
Endnote: The Determination of a ................................... ............... 62
A different approach to the weak nonlinear analysis...................... ............... 64

5 THE RAYLEIGH-TAYLOR PROBLEM WITH PHASE-CHANGE...... .....73

Introduction ................... ...................................... ......................... 73
Physical and M them atical M odel.................................................................... 74
Analytical Results for Infinitely Deep Fluids....................................................78
Numerical Results for Finite Depth Fluids.......................... ...............84

6 EFFECT OF VAPOR DYNAMICS ON THE RAYLEIGH-MARANGONI
PROBLEM ........................................................89

Introduction ...................................... ............. ...... .........................89
Experimental Apparatus and Procedure .......................................91
Experim mental Results .................. ........................... ......... .. ...... ................. 95

7 CONCLUDING REMARKS AND FUTURE SCOPE.............. ............ 102

A ctive V apor L ayers .............. ..... ........... ....... ..... ........... .. .................... 102
Conditions at the Interface of an Already Evaporating Liquid............................. 104
Coupling of Evaporative Instabilities with Other Phase-Change Phenomenon....... 104

APPENDIX

A THERMODYNAMIC EQUILIBRIUM RELATION AT THE INTERFACE........ 105

B A SIMPLE DERIVATION TO OBTAIN THE CRITICAL CONDITION FOR
ONSET OF CONVECTION IN A FLUID LAYER.................................108

C A SIMPLE DERIVATION FOR THE RAYLEIGH TAYLOR INSTABILITY ....111

D THE BOUSSINESQ APPROXIMATION .................. .... ............113

E THE UNIT NORMAL OF A SURFACE............................. 116

F DERIVATION OF THE SURFACE SPEED...........................................................117

G THE INTERFACIAL ENERGY BALANCE.......................................... ..................119









H THE EXPANSION OF A DOMAIN VARIABLE AND ITS DERIVATIVES
A L O N G TH E M A PPIN G .................................................................................... 12 1

I N U M ER IC A L M E TH O D ................................................................................... 126

L IST O F R E F E R E N C E S ............................................................................................ 129

BIOGRAPHICAL SKETCH ............................................................... ...............131
















LIST OF TABLES

Table page

3-1 Physical properties of water and water vapor at 1000C .......... ..... .... ........ 31

3-2 Definitions of R R2 and R3 .................. ............................. ......31

3-3 Comparison of critical temperature differences ........................................... 32

3-4 Definitions of R4, R5 and R6................................................................ .....................36

6-1 Physical properties of silicone oil and air at 200C ...................................... 96
















LIST OF FIGURES


Figure page

1-1 Physical model of the phase-change problem....... ...............................2...

1-2a H eat input from the liquid side....................................................3

1-2b Heat input from the vapor side................................ ........ .. .. .................4

1-3 Physics with fluid dynamics..............................................5

1-4 Rayleigh convection ................... ........... .. ........ ........... .... ....... .6

1-5 M arangoni convection..............................................8

3-1 R' vs. wavenumber. ...................... ............. ........33

3-2 R 2 and R3 vs. w avenum ber. ........................................................... 36

3-3 R4 through R6 vs. wavenumber. ............................... ............... 38

3-4 The velocity vectors in both phases at the onset of instability. z = 0 represents the
interface position in the base state................................... ....... 40

3-5 The critical temperature difference versus the dimensionless wavenumber for
water-water vapor system where d = 3mm and d* =1.5mm...............................41

3-6 Evaporation number vs. dimensionless wavenumber for two water-water vapor
system s ............................................................43

4-1 Heat transfer across the bottom plate. ....................................... 59

4-2 The scaled pressure perturbations in the vapor domain for the case where d = Imm,
3 = 0.1 and a = 1. z = 0 represents the interface position in the base state..........60

4-3 The scaled pressure perturbations in the vapor domain for the case where d = 1mm,
S= 0.3 and o = 1. z = 0 represents the interface position in the base state..........61

4-4 The scaled pressure perturbations in the more viscous vapor for the case where
d = Imm, 3 = 0.3 and o) = 1. z = 0 represents the interface position in the base
state ..................................................62









5-1 The physical model .................................. .............. ........ 75

5-2 The dependence of the phenomenological coefficient on the wavenumber ............83

5-3 Com prison of m odels ................................................. ............... 85

5-4 Effect of shallow fluid layers ........................................... ............... 87

5-5 Critical temperature difference vs. dimensionless wavenumber.............................. 88

6-1 Different types of convection mechanisms ........................................89

6-2 Schematic of the experimental setup............................................. .. ......... 91

6-3 Experim mental setup ................................................................. .. ......... 93

6-4 z-component of velocity in a silicone oil-air system at the onset of instability
where the silicone oil and air depths are 5mm and 3mm, respectively.................97

6-5 z-component of velocity in a silicone oil-air system at the onset of instability
where the silicone oil and air depths are 5mm and 5mm, respectively.................97

6-6 Flow profiles in a silicone oil-air system at the onset of convection where the
silicone oil and air depths are 5mm and 3mm, respectively. ................................98

6-7 z-component of velocity in a silicone oil-air system at the onset of instability
where the silicone oil and air depths are 5mm and 9mm, respectively.................99

6-8 The critical temperature difference vs. air height for a bilayer system where the
silicone oil depth is 8mm and the aspect ratio is 2.5................ ................100

6-9 Pictures of experiments taken by IR camera at the onset of instability ................. 100

B-l Critical Rayleigh number vs. dimensionless wavenumber ..................................110
















Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy

EFFECT OF VAPOR DYNAMICS ON INTERFACIAL INSTABILITIES

By

Ozgur Ozen

December 2004

Chair: Ranga Narayanan
Major Department: Chemical Engineering

Instabilities at fluid-fluid interfaces are of interest to researchers on account of their

role in many industrial and natural processes. These instabilities are often accompanied

by phase-change between liquids and their vapors, and in some instances are caused by it.

Examples of processes involving such instabilities are film boiling and drying processes.

An everyday occurrence is seen in the pattern formation caused by the drying of coffee

drops on a kitchen counter.

It turns out that the instabilities are accompanied by fluid flow in both phases but

many past studies of the phenomena have been done by considering only the liquid fluid

dynamics while the role of the vapor phase has been either left out of the study or

relegated to heat transfer only while its fluid dynamics has been ignored.

In this study, evaporative instabilities are analyzed theoretically and experimentally

at the interface between a liquid and its vapor where both the liquid and the vapor are

fluid dynamically active in order to investigate and understand the role of vapor

dynamics on evaporative and interfacial instabilities. In the theoretical model, the Navier









Stokes and energy equations in the liquid and vapor domains were using linear and

nonlinear perturbation methods with the assumption that both fluids are incompressible

using the Boussinesq approximation. The linear theory tells us the conditions at the onset

of instability and the critical point at which the interface becomes unstable. The weak

nonlinear method tells us the behavior of the system slightly beyond the onset of the

instability. The linear and weak nonlinear theories are of importance for the experimental

reasons as they provide us with the information on how and when the instability will

occur.

In a separate study of non-evaporative systems, experiments were performed in a

liquid-gas bilayer system in a circular container with the intention of studying instabilities

in the absence of evaporation. The depth of the gas layer was changed to observe the

effect of the fluid dynamics in the gas layer on the onset flow pattern and conditions.

These experiments were then compared with a corresponding theoretical model and it

was concluded that the gas dynamics played an ever-increasing role as the gas layer

height was increased. Past workers had assumed the vapor to be a conductive, fluid

dynamically passive medium. However, the experiments show us that deep enough

vapor layers can cause the bilayer system to convect earlier than expected by thermally

coupling with the liquid at the interface.

In summary, from theoretical studies it was found out that active vapor layers play

an important role on the evaporative instabilities even in the absence of natural

convection. Moreover, it is discovered that the fluid dynamically active vapor plays a

very strong stabilizing role on the evaporation problem. In addition to that, from the

experimental studies in the absence of phase-change it was revealed that the vapor layers









can and do affect the onset flow pattern and conditions. Both theoretical and

experimental studies also indicate that the commonly used infinitely deep vapor layer

assumption is not always true and can lead to errors.














CHAPTER 1
INTRODUCTION

Evaporation and condensation involve the transformation of a substance from the

liquid to the vapor state or the other way around by the addition or depletion of energy. It

turns out that these processes can lead to an interfacial instability. This instability is

typically manifested by the sudden undulation of an erstwhile flat surface between a

liquid and a vapor as a control variable, such as the temperature drop across the bilayer,

exceeds a critical value. There is a need to understand the phenomenon of the instability

in the light of the application of evaporation to many technologies and this is the focus of

our work. Evaporation arises in a variety of industrial and natural processes, such as in

heat pipe technology, the drying of oxide films, the drying of lakebeds and in dry-eye

syndrome.

However, evaporation comes accompanied by fluid flow and natural convection

due to buoyancy and interfacial tension gradients. In order to have an understanding of

the various phenomena associated with evaporation and the first order effects that

promote instability during evaporation we advance some "picture arguments," first for

the pure phase-change problem and later for connected instabilities such as convection

and gravitational instability.

Physics of Various Instabilities Connected to Vapor-Liquid Phase Change

Evaporative Instability

Evaporative instability is best understood by considering an evaporating liquid in

contact with its own vapor. To get an idea of the physical model the reader's attention is









drawn to Figure 1-1, which depicts a liquid-vapor bilayer where a temperature difference

is imposed across it. The input variables are the fluid depths, the temperature difference

and the pressure of the vapor at the upper plate. As a result of fixing these variables the

evaporation rate, the fluid behavior and the interface morphology can be determined. It

may be noted that evaporation can be obtained in the absence of a temperature difference

by manipulating the pressure at the upper plate. However, in this work we shall imagine

that a temperature difference is always imposed.

The interface is originally thought to be flat and it is the stability of this state that is

in question. The liquid and vapor depths may be fixed by suitably adjusting the liquid

feed and vapor removal rates. In fact, given a temperature drop the upper plate pressure

may be adjusted in such a way that the evaporation or condensation rate is zero. We will

consider much of the analysis in this thesis for a case where initially there is no mass

transfer across the interface as it simplifies the calculations and is sufficient to reveal the

physics behind this problem. We now proceed in a stepwise manner to consider

evaporative instability and return to Figure 1-1.

r COLD PLATE r

LG Vapor



LL Liquid

JA HOT PLATE

Figure 1-1. Physical model of the phase-change problem

First, imagine that the base flow rate is zero and the fluid dynamics is left out. By

imagining that the fluid mechanics is omitted, we are also pretending that pressure










perturbations do not arise. This will obviously change when the perturbed fluid motion is

taken into account. Interesting conclusions can be drawn by simply looking at the

temperature profiles in the base state. The schematics in Figures 1-2a and 1-2b depict an

evaporation problem at a zero evaporation rate. Figure 1-2a shows the heating

arrangement where the heat needed for evaporation is supplied from the liquid side while

Figure 1-2b shows heat input from the vapor side. The arrows in the figures show the

direction of the front motion if the base state were to have some evaporation though this

is not central to the physical arguments presented here.

COLD PLATE

Vapor




Front motion Liquid


HOT PLATE
Figure 1-2a. Heat input from the liquid side

In both figures the dotted wave represents a perturbed interface while the dotted

lines represent the temperature profiles in the perturbed state at the interface. Note that

the temperature at the perturbed interface is the same as the unperturbed state because the

temperature is assumed to be at equilibrium with the pressure at the interface, an

assumption that is justified at low evaporation rates. When the fluid dynamics is ignored

the pressure is not perturbed and consequently also the temperature is unperturbed at the

interface (see Appendix A for a description of the thermodynamic equilibrium relation at

the interface and its development upon the imposition of an infinitesimal perturbation).

When heated from the liquid side, a crest is further away from the heat source and also

closer to the heat sink. Therefore the temperature gradient on the vapor side gets sharper










whereas the temperature gradient on the liquid side becomes less sharp (i.e., less heat is

transferred to that point and more heat is taken away). The rate of evaporation from the

crest is thus decreased which makes the interface unstable. On the other hand, when

heated from the vapor side, the trough gets closer to the heat sink and is further away

from the heat source. Consequently, less heat is transferred to the front and more heat is

taken away from it. The speed of the trough is thus decreased and as a result the front is

stabilized. Note that the picture argument presented here suggests that when the interface

is unstable, it is unstable for disturbances of all periodicities or wavenumbers and when it

is stable, it is stable for all wavenumbers.

HOT PLATE

Vapor




Liquid Front motion

COLD PLATE
Figure 1-2b. Heat input from the vapor side

Now include the fluid dynamics in both phases in the thought experiment, but do

not consider any other convective effect. Figure 1-3 is similar to Figure 1-2.a except that

fluid dynamics is now taken into account in both phases. Observe that a deflected surface

does not have the same temperature as an erstwhile flat interface for a trough is expected

to be at a higher temperature than a crest as it is closer to the heat source. Remember also

that we started out with a base state where there was thermodynamic equilibrium at the

interface and no phase-change was taking place. Now that the trough is at a temperature

higher than the saturation temperature there will be evaporation taking place there and

consequently, when perturbed, the vapor will condense into its own liquid at the crest










where the temperature is colder than the saturation temperature. There will be upward

flow at the troughs and downward flow at the crests as indicated in Figure 1-3. Again,

instability occurs at every wavenumber.

COLD PLATE



(- N






HOT PLATE
Figure 1-3. Physics with fluid dynamics

Now surface tension acts on pressure differences in fluid dynamics problems via

the surface curvature and in many phase-change instability problems it is the surface

tension that comes in to stabilize the situation. Its effect is felt most strongly at small

wavelengths and weakly on lazy perturbations. But, surface tension comes attended by

the complication that it depends upon temperature and we will take this into account as

we continue to study the problem in stages. There is also another complication and this

time it is due to gravity. If the vapor were at the bottom as in pool boiling then gravity

would play the destabilizing role of a heavy fluid on top of a light fluid. The

destabilization would take place in that case when gravitational potential energy

overwhelms surface potential energy. This is the familiar Rayleigh-Taylor instability

problem and we give a brief discussion of this instability later in this chapter. However,

even if the vapor were on the top and the liquid were heated, gravity still plays a role via

buoyancy and to understand the nature of natural convection due to buoyancy and

interfacial tension gradient we turn to more picture arguments.









Before we move on to a discussion of these other effects, it might be noted that

evaporative instability can take place even if the interface were assumed to be flat. In

this case, the instability would be manifested by fluid motion. This case will not be of

interest to us.

Rayleigh Convection or Buoyancy-Driven Convection

Buoyancy-driven convection, often termed natural convection or Rayleigh

convection, occurs when a fluid is subject to a temperature gradient in a gravitational

field and when there is a variation of density with respect to temperature. Imagine a layer

of liquid bounded vertically by two horizontal surfaces, with the lower surface at a

temperature greater than the upper one as in Figure 1-4.

COLD


Figure 1-4. Rayleigh convection

As density typically decreases with an increase in temperature, the fluid near the

top surface is heavier than the fluid at the bottom plate, creating a top-heavy arrangement,

one that is gravitationally unstable. For sufficiently small temperature differences, the

fluid simply conducts heat from the lower surface to the upper one, creating a linear

temperature field across the fluid. Now imagine a mechanical perturbation being

imparted to an element of fluid so that it is projected downward. As the density of the

"fluid packet" is greater than its environment it will proceed downward. Due to mass

conservation, fluid from below will be displaced and will move upward. This motion









will continue unabated unless the fluid's kinematic viscosity and thermal diffusivity are

large enough that the mechanical perturbations and thermal perturbations die out quickly

and the original mechanically quiescent state is re-established. From this thought

experiment it is immediately clear that unlike many classical problems of fluid

mechanics, the problem of the onset of convection can not be considered for fluids in the

limit of vanishing viscosity or thermal diffusivity. In other words, disturbances to the

fluid are dampened by the kinematic viscosity and thermal diffusivity of the fluid and

enhanced by the buoyancy forces. Under steady conditions there is a balance that is

established between these effects, represented by the dimensionless Rayleigh number, a

group that arises from the scaling of the modeling equations. It is given by


Ra= agAd3 1.1
VK


Here, ac is the negative thermal expansion coefficient and is usually positive, g is

the magnitude of gravity, AT is the vertical temperature difference across the fluid layer

and is taken positive when the fluid is heated from below, d is the depth of the fluid, v is

its kinematic viscosity, and K is its thermal diffusivity. If the temperature difference is

increased beyond what will be called the critical temperature difference, then the

gravitational instability overcomes the viscous and thermal damping effects and the fluid

is set into motion, causing buoyancy-driven convection. (See Appendix B for an

elementary derivation of the critical condition for the onset of buoyancy driven

convection.) If the upper surface adjoins a vapor phase then convection can initiate in

either phase. If it begins in the lower phase, it will propagate to the upper phase by way

of viscous dragging. If convection is simultaneous in both phases they can become









thermally coupled with the interface becoming an isotherm. If it initiates in the upper

phase, it cannot propagate to the lower phase except by way of generating interfacial

tension gradients, which in turn will generate weak flow in the lower liquid. For that

reason, we turn to interfacial tension gradient-driven or Marangoni convection in the

following section.

Pure Marangoni Convection or Interfacial Tension Gradient-Driven Convection

Interfacial tension gradient-driven convection or Marangoni convection is unlike

buoyancy-driven convection and can occur in a fluid even in the absence of a

gravitational field. Imagine a layer of fluid which is bounded below by a rigid plate and

whose upper surface is in contact with a passive gas as depicted in Figure 1-5.





HOTu







Figure 1-5. Marangoni convection

A passive gas is not one without density or viscosity. It is merely a gas wherein

one tacitly assumes that thermal or mechanical perturbations cannot occur. Pretend that

above the passive gas, there is another rigid plate. For the sake of consistency with the

earlier argument, let the lower plate be at a temperature higher than the upper plate's

temperature. Now, imagine that the interface between the lower liquid and the passive

gas is momentarily disturbed. The regions of the interface that are pushed up experience

a cooler temperature. Likewise, the regions of the interface that are pushed down









increase in temperature. Typically, interfacial tension decreases with an increase in

temperature. Therefore, the regions of the interface, which are pushed up increase in

interfacial tension, which pulls on the interface, while the regions of the interface, which

are pushed down, decrease in surface tension. When the fluid is pulled along the

interface, warmer fluid from the bulk replaces the fluid at the interface near 'x' enhancing

the interfacial tension gradient-induced flow. If the temperature difference across the

liquid is sufficiently small, then the thermal diffusivity of the fluid will conduct away the

heat faster than the disturbance can amplify. The dynamic viscosity will also resist the

flow causing the interface to become flat again and the interfacial tension to become

constant. As in the case of buoyancy-driven convection, there exists a critical

temperature difference when the surface tension gradient-driven flow is not dampened by

the thermal diffusivity or viscosity, and the fluid is set into motion. Scaling of the

equations that model the interfacial momentum balance leads to a dimensionless group

that exhibits the balance between dissipative effects of the viscosity and thermal

diffusivity and the promoting effects of interfacial tension gradient. This dimensionless

group is called the Marangoni number and is given by


Ma = )/TATd 1.2
IUK


Here y, is the change in the interfacial tension with respect to the temperature, and

[t is the dynamic viscosity of the fluid. It is noteworthy that while the above argument

assumes the interface to be deflectable it is not necessary for the interface to undulate in

order for Marangoni convection to take place. To see this a little more clearly, look at

Figure 1-5, but pretend that the interface is flat. Points 'x' and 'y' will therefore be on









the same horizontal level. Imagine an element of fluid projected downward from the

interface at point 'y'. Fluid nearby must replace the lost fluid at 'y' and this fluid comes

from the region near 'x'. Consequently, hot fluid from below moves to 'x'. If the

interfacial tension decreases with temperature, as is the case in most substances, then

fluid flows from warm surface regions to cold surface regions. Therefore, fluid moves

from 'x' to 'y' and the convection continues unless large thermal diffusivity and viscosity

dampen out the disturbances. This means, theoretically speaking, that Marangoni

convection can be analyzed assuming a flat free surface, i.e., a surface with an infinite

surface tension. This fact was observed by Pearson (1958) who gave the first theoretical

analysis of Marangoni convection using a methodology similar to that used by Rayleigh

for the case of buoyancy driven convection.

The physical argument presented above can easily be extended to a bilayer system

where the gas layer is fluid dynamically active. In this case, the gas layer will also

convect due to surface tension gradients and will bring in colder fluid from the cold plate

to 'x' opposing the destabilizing nature of the convection in the liquid; however, the

effect of the convection in the gas is less because of its high kinematic viscosity.

The extent of either Rayleigh convection or Marangoni convection is primarily a

function of the fluid depths. By examining equations 1.1 and 1.2, we notice that

Rayleigh convection is proportional to the cube of the fluid depth and that Marangoni

convection is directly proportional to the fluid depth. From these scaling arguments, we

can conclude that for deep fluids, buoyancy-driven convection is more prevalent, whereas

for shallow depths, interfacial tension gradient-driven convection is more prevalent, and

at intermediate fluid depths, both Rayleigh and Marangoni convection can occur. If the









interface is allowed to deflect one might anticipate that in the case of pure Marangoni

convection, warm fluid flows towards a trough while in the case of Rayleigh convection

cold fluid flows away from a trough leading us to believe that Marangoni motion can

counter Rayleigh driven motion. This counter effect was demonstrated computationally

by Sarma (1987) and could not be seen in the study by Nield (1964) wherein flat surfaces

were assumed.

When a layer of fluid is heated from above, it creates stable density stratification.

Therefore the buoyancy force does not cause convection. Marangoni convection, though,

may still occur in fluids being heated from above. For example, consider a single layer of

fluid superimposed by a passive gas and heated from above. One might imagine that

purely steady Marangoni convection would not occur. For instance, suppose a random

perturbation causes some part of the surface to become warmer than the rest of the

surface. The interfacial tension will decrease in this region and the higher interfacial

tension elsewhere will pull fluid away from this hot spot. By continuity, fluid lying

below the hot spot will move up to replace the displaced fluid. As the lower fluid is

cooler than the surface, this region now cools off and the interfacial tension increases,

thereby either stabilizing the region or causing an overshoot in temperature of that spot

and generating oscillatory motion. Only calculations can tell precisely which case will

occur. In real systems the upper fluid is never truly passive. Given the same scenario,

fluid movement along the interface will also drag warmer fluid from above and cooler

fluid from below. Depending upon the fluid properties and depth ratios the instability

will or will not be reinforced.









Oscillatory convection can occur if the fluid is heated from above and buoyancy

effects are included. This was shown by Rednikov et al. (1998) for the case when the

upper fluid is passive. The explanation is as follows. If a small volume of fluid is

displaced within the bulk of the fluid, the density stratification acts as a restoring force.

The Marangoni force acts similarly and the two effects reinforce each other causing an

overshoot and sustained oscillations. Indeed, as was demonstrated in their paper, this

only occurs in certain fluids, at certain depths, where the buoyancy and interfacial tension

forces are comparable.

There is another phenomenon associated with surface tension gradient-driven

convection, often called the long wavelength Marangoni instability. This instability

typically occurs when either the surface tension or the depth of the fluid layer is very

small. The initiation of this instability is similar to the description given above.

However, in the long wavelength scenario, the convection cells are much larger than the

regular, or short wavelength, Marangoni convection. As the convection propagates, it

causes large-scale deformation in the interface, which can actually cause the interface to

rupture; that is, the interface deforms to such an extent that it comes in contact with the

lower plate. This phenomenon occurs in the drying of films and in coating processes

(Schatz).

Rayleigh-Taylor Instability

So far we have talked about the liquid and gas systems where the liquid in

consideration underlies the gas phase. There are however situations, such as film boiling,

where the opposite arrangement is encountered. When the vapor phase underlies a liquid

phase, this is a top-heavy arrangement called Rayleigh-Taylor instability and it is

potentially unstable. The instability occurs when the gravitational potential energy









exceeds the stabilizing effect of surface potential energy. In the classical treatment of the

problem, it is assumed that there is no heat and mass transfer, and it is well known that

the system is unstable to disturbances with a wavenumber smaller than the critical

wavenumber given by







where y is the surface tension, p and p* are the liquid and vapor densities, respectively.

A derivation for the classical Rayleigh-Taylor critical wavenumber from the book by

Chandrasekhar (1960) for infinitely deep fluid layers is given in Appendix D. This

instability is also discussed in connection to the phase-change problem in Chapter 5.

What makes the Rayleigh-Taylor problem of special interest to us is the unavoidable

existence of both the heat and mass transfer in film boiling.

Plan of the Thesis

The rest of the thesis contains the modeling equations, their scaling and the method

of solution to determine the onset conditions of the instabilities. This is followed by a

discussion of the results of calculations for the onset of instability in the evaporation

problem, the Rayleigh-Taylor problem and the pure convection problem. In each of these

discussions, the role of vapor dynamics is stressed: a role that has not been emphasized

up to this point in this chapter. The vapor is a source of heat transfer in all of these

problems. To assume that the vapor plays a passive role is then ordinarily incorrect.

When the vapor moves it does so on account of several mechanisms. It may move

because pressure perturbations are imparted to it or because the thermal profiles generate

phase-change at least upon perturbation, or it can move on account of viscous dragging.






14


Whatever the mechanism for its motion, it must change the mode of heat transport from

conduction to convection. This change in the transport must feed back to the interface

and either assist the undulations or hinder them. Our goal is to understand the role of the

feedback, the modeling and the physics of instability mechanisms.

This thesis is different from the traditional type of dissertations because of the

absence of a chapter on the literature search. I have decided to refer to the literature in

the context of the work presented.














CHAPTER 2
THE MATHEMATICAL MODEL

This chapter includes all of the equations used to analyze the problems described in

Chapter 1. The equations are divided into the field equations that are meant to describe

the behavior of the state domain variables and the interfacial equations that describe the

relation of the domain variables with surface variables such as surface speed and surface

curvature which, in turn, depend upon surface position. The domain equations and

interfacial conditions are written in a general manner without regard for distinction

between the various problems. We shall draw the attention of the reader to the

differences as we go along and thereby indicate how each problem belongs to a special

case of the general equations.

What may not have been apparent from the discussion in the earlier chapter is that

the instabilities are related to a nonlinear property of the modeling equations. The

nonlinearities arise on account of more than one reason. One of the reasons to expect

nonlinear characteristics in the modeling equations is that the interface position is

coupled to the transport behavior and the two depend on each other. Another reason is

that the domain variables are coupled. For example, the interface temperature

distribution affects the pressure distribution, which in turn affects the flow profiles. The

flow profiles then affect the temperature profiles, either reinforcing them or causing the

disturbances to subside. We now turn to the nonlinear equations and shall point out their

characteristics as we go along.









Nonlinear Equations

The various problems we will deal with are all fluid flow problems; hence, the

equations that model the fluid physics in unsealed form are given by the Navier Stokes

equations in each phase, i.e., by


p -+vVV =VP+p(T) g+, V2v 2.1

and


PK +V.VV VP*+p(T*) g+,U V2V 2.2

Each of the problems is connected with temperature gradients. These affect the

energy transport. Therefore the energy equations in each phase are needed. They are

given by

dT
-+v VT = KV2T 2.3
Bt
and


++ VT = Kr V2T* 2.4
9t
Mass conservation in each phase requires the continuity equations. They are

V-= 0 2.5
and

V *= 0 2.6
In the equations above, v, P and T are velocity, pressure and temperature fields,

respectively, and the asterisk denotes the vapor phase. In the equations, p, [t, and K are

the density, viscosity and thermal diffusivity, respectively.

A major assumption that has tacitly been made is that the densities of the fluids are

treated to be constant in all of the terms save those associated with gravitational









acceleration. This assumption is fair provided the temperature drop is not very large or at

least as long as the total density change divided by its mean value in each phase is a small

fraction. This is the celebrated Boussinesq approximation. It can be justified by a

perturbation series in a scaled thermal expansion coefficient. This justification, given by

Mihaljan (1962), has proven to be invaluable in studies in thermal convection greatly

simplifying the calculations. More on this approximation is given in Appendix E.

We continue with the modeling equations and assume that the bottom and top plate

temperatures are kept constant and therefore T (-d) = Tbottom and T' (d*) = T hold.

The top and bottom walls are impermeable, therefore, v- n -d = 0 and

v* n = 0 hold. The no-slip condition applies along the top and bottom plates, and


gives rise to t -d = 0 and =d. = 0.

At the interface, the mass balance equation is given by

p(v-u) = -n = (v* -u) -n 2.7

OZ -
-- +k
where ni = is the unit outward normal (Appendix F) and




az
u n =- is the interface speed (Appendix G). When there is no phase-




change at the interface, we set the interfacial mass balance equal to zero, and we obtain

two interface conditions from the interfacial mass balance equation:

p(V-u)i- = 0









and

p (v8 -u )- =0

At the interface, the temperatures and the tangential components of velocities of

both fluids are equal to each other, so T = T* and v t = v* I t hold where

OZ-
I +-k
dx k
t = is the unit tangent vector.




The interfacial tension and its gradient with respect to the temperature at the

interface come into the picture through the force balance. By taking the dot product of

the stress balance with the unit normal and unit tangent vectors separately, we get the

normal and tangential stress balance equations. The force balance is given by

p V (v; -i)- j i+;y2Hh+Vy i =(p V ( -U) -T ii 2.8

where T = -PI + S is the total stress, y is the interfacial tension and

02Z
nx2
2H = is the surface mean curvature.




In addition we need the interfacial energy balance, viz.


H+ (_- iiv-+ +q-S.-v- ii
2.9
= pH +(v -u) (v -u)+q -S -(v -u) -n

Here H is the enthalpy per unit mass, S is the extra stress, and q is the heat flux. Now,

if there is phase-change across the interface, we are short of one equation as the









interfacial mass balance (Equation 2.7) gives us only one relation. To get this additional

equation we observe that when the evaporation rate is very small we may assume that

there is thermodynamic equilibrium at the phase-change boundary, that is


P*-Pbase= p'hln 2.10
Base

where Pbase and Tbase are the base state pressure and temperature of the fluid, and h is the

latent heat of evaporation per unit mass.

Dimensionless Nonlinear Equations

The nonlinear equations in dimensionless form, if necessary, are given below. The

/. d 'Uv
following scales are used: length, d; velocity, v = ; time, _; pressure, P = and
d v d

T-T
temperature, T = -old where AT = Thot Tld The instability in these problems is
AT

caused by the temperature gradients, so it is reasonable to scale the velocities with respect

to the thermal diffusivity since it is the dissipation of the temperature perturbations that is

important. The temperature is scaled with respect to the total temperature difference

across the bilayer system and the reference temperature is the cold temperature, which in

many systems we cover will be the top plate. This will prove to be of convenience when

we move to higher order expansions of the problem. Throughout the course of our study

we will keep the top and bottom plate temperatures constant and by scaling the

temperature with respect to our choice of reference temperature and the total temperature

difference, we will significantly simplify the top and bottom temperature boundary

conditions. Hereafter, we will drop the tilde on the scaled temperature. Using the

expansion for the density given in Appendix E, the momentum equations in each phase

become










-+-V = -VP -G + Ra T +V2 2.11
Pr At

and


-+v -V v VP VGGi + Ra T* + V2v 2.12
v Pr A 'U v

V gd3
where Pr = and G = are the Prandtl and Gravity numbers, respectively.
K VK

The energy equations in each domain are given by


V2T = T + -'VT 2.13

and


VT +V V T 2.14

The mass balance at the interface, z = Z(x, t), now becomes



p
-)n= ^ *-a* 2.15

The normal stress balance in dimensionless form is given by


Ca v (V -) T : n + 2H = Ca p *T :i fi 2.16
Pr p Pr P/

where Ca = K is the Capillary number.
yd

The tangential stress balance becomes

S :lit -S:t + MaVT =0 2.17

In addition, the interfacial energy balance becomes










C l+K'41 -(V -U2 ii)-.i -EKvT.-iik VT rn

( /* / ^2.18



2 C AT VK
Here, KP = E= p and Vp =h d2 E is the Evaporation number.


At the interface, the temperatures of both fluids are equal to each other, and the

thermodynamic equilibrium at the phase-change boundary holds, that is


HKE (P*-PBASE) =ln T 2.19
Tbase

P VK
where IKE = Vc
p h d2

Some Observations

We can make several observations on the above equations. First, as noted, the

interfacial mass balance is simplified when there is no phase-change and the

thermodynamic equilibrium relationship is then ignored. It is also possible to consider

phase-change with a zero temperature drop across the bilayer. However, in this case a

non-zero evaporation base state is required. This case is not considered in this study.

Second, when phase-change takes place and the thermodynamic equilibrium is taken into

account, the effect of curvature on latent heat also needs to be taken into account.

However, the effect is weak and this fact is proven in Appendix A. Third, by imposing a

restriction of a flat undulated interface and infinite interfacial tension, the equations

simplify considerably. The normal component of the stress balance is then used only to

determine the value of the curvature times the interfacial tension but not to close the

problem. Fourth, if the Marangoni effect were excluded, then the tangential stress









balance at the interface simplifies by rendering the interfacial tension gradient to zero but

no singularity arises on account of this. Finally, the interfacial balances are all frame

invariant. This is seen in the force balance if the mass balance is taken into account and

it is seen in the energy balance in the form that it is presented. The energy balance is

derived for the reader's benefit in Appendix H.

The Linear Model

To see the equations that determine the onset of instability requires us to linearize

the equations about a base state. First, a word about the base state: The base state is

always a solution to the nonlinear equations and often it might seem defeating to look for

a base state if it means solving a full set of intricate coupled nonlinear equations. But in

practice for a large class of problems the base state is seen almost by inspection or by

guessing it. For example, in the problems considered in this thesis an obvious base state

is the quiescent state with a flat interface. It is not in global thermodynamic equilibrium

as a temperature gradient is imposed but it is a state whose stability is called into

question.

To study the stability of the interface, arbitrary disturbances of infinitesimal

amplitude are applied and for a given set of input variables, the growth and decay time

constants of these disturbances or equivalently the critical temperature difference across

both layers that result in marginal stability are determined. The reader might wonder

why the disturbances are taken to be small and what this restriction implies. If a state is

unstable to infinitesimal disturbances it must be unstable to all disturbances. The same is

untrue of finite amplitude disturbances. Thus, stability to small disturbances usually

means very little but instability implies a lot. Using small disturbances also allows the

dropping of terms with quadratic interactions and higher degree interactions. Therefore,









the plan is to linearize the above equations about a known base state and to find the onset

of the interfacial instability from the perturbed model.

Perturbation Equations

Hereon, the variables of the base state are denoted by the subscript 0, and the

variables of the perturbed state by the subscript 1. Thus, the temperature, when

perturbed, is described by


T= T +s T, + dT z++O(2)
dz )

where E is a small perturbation parameter representing the deviation from the base state

and z, is the mapping of the perturbed configuration onto the reference configuration. Its

meaning is explained in the Appendix I and, at the interface, it is simply the perturbation

of the surface deflection to first order, Z1, a variable to be determined during the course

of the calculation. We can further expand T, and other subscript 'one' variables using a

normal mode expansion. For example,

T, = T, (z) e'te .

Here, a is the inverse time constant while wo is a wavenumber associated with the given

perturbation. A wavenumber arises because the system is infinite in lateral extent.

However, the wavelength, which is proportional to reciprocal of the wavenumber, may be

related to the spacing between the sidewalls in containers with periodic boundary

conditions.

The same expansion as above is used for both components of velocity and pressure

as well as temperature in both phases. It also holds for the interface deflection variable









Zi. The first priority is to get the base state, as this is the state whose stability is in

question and about which the expansions must be done.

The Base State Solution and the Perturbed Equations

In the base state there is no flow in either the liquid or the gas. For the phase-

change problem, assume that there is no phase-change. In other words, assume that the

upper plate pressure is fixed so that the phase-change rate in the base state is zero. This

of course does not mean that the phase-change rate will be zero upon perturbation.

Consequently, both fluids are in a conductive state, thus V, = v0 = 0 Hereon, we will

drop the overbar of the scaled temperature. The temperature profiles in the liquid and its

vapor in dimensionless form become

k* k
T=- z+
k +k* k +k*

and

k k9
T z+
k +k*z k 8+k*


where = d-
d

After we perturb the domain and boundary equations in the manner given in the

previous section, we arrive at the following equations in the domain in dimensionless

form:

d2 1
w2 v 1 P2 d =U- vx1 2.20
dz Pr


d2 2 W2vz1 + Ra T, cr v z 2.21
dz Oz Pr












dz


dT0
dz Vz =1 T,
dz


dvX
dz


for the lower phase,





d 2


and

Sd2 2 1 V
d2 2 Vx1 -i w P1= vPry
dz U Pr vP


2 V*z a Ra Ta T -* 1 V
a v 'U dz Pr v

d2 dT* K
dz2 i T d vz z I T*
Idz ^ A dz A:


dv1
w vx, + = 0
dz


for the upper phase.

At the lower plate, z = -1, the liquid is subject to vI = 0, vx,


while the conditions at the top plate become v*I = 0, vI = 0 and T,1

At the interface the mass balance turns into


P P
vz --vzi 1-j I Z-
P P

which, in the absence of phase-change, becomes

Vz1 = 7 Z,


O and T, = 0


0.





2.28


vTi = C Z,.

The continuity of temperature and the no-slip conditions become


2.22


2.23


2.24



2.25



2.26


2.27










T, + dTo, = T + 1O Z, 2.29
dz dz
and

vx, = v,. 2.30
The normal and tangential stress balances assume the following dimensionless

form:


Ca P,-P1)-2Ca dvzz dvz (Bo+w2) Z,=0 2.31
aP y dz u dz

and


dvx +w v dvx +iw vz iw Ma T, + dTo Z =,1 0 2.32
,u dz dz 1 dz )

where Bo = p -)gd2 -, and again Ca = K and Ma = AT d are the Bond,
Y yd /UK

Capillary and Marangoni numbers, respectively. The interfacial energy balance is


vz E k dT dT1 = Z, 2.33
k dz dz

C AT
where E stands for the Evaporation number and is given by E = C, is the liquid


heat capacity. Again, one has to be careful in order to avoid a mistake here. The

interfacial energy balance in the absence of phase-change simply becomes

k* dT1 dT,
k dz dz

The equilibrium condition at the interface, when there is phase-change across the

interface, is

dT*
T1 + oZ
P -n Zdz 1 2.34
KE TbPE 1 se
Base










where HPE g d is a dimensionless parameter from the linearized Clapeyron equation,
h

and Tase is the scaled base state temperature of the vapor at the flat interface and it is

equal to TO.

There are several observations to be made in the perturbed equations. First, the

growth constant a appears in the domain as well as interfacial equations. Dropping them

in the domain equations lead to an interfacial mode, however dropping them may not

always be justifiable but when the growth constant is dropped from the domain its

calculation becomes straightforward. It must then also be real as it is given by the ratio

of two determinants. Second, the growth constant is ordinarily an eigenvalue for a

sensible control parameter such as E. Third, it is conceivable that solutions exist such

that Zi is zero but they may not be the first to become unstable as a control parameter,

such as E, changes. Fourth, when we pay closer attention to some of the dimensionless

numbers we recognize that their values are very small. A typical value for Capillary

number is of the order of 10-8 and for HKE is of the order of 10 '0. This indicates that

for systems where there is a high surface tension between the two fluids it may not be a

bad idea to set Ca equal to zero, which consequently will imply that the surface deflection

Z, is equal to zero. That is to say that we can assume the interface to be flat at all times

and the instability will still set in because the temperature and velocity will still affect

each other through the domain equations. The same goes also for the HKE term.

In the next chapter the perturbation problem is solved for the special case of pure

evaporation in the absence of natural convection, and from the calculations a physical

interpretation of the problem is sought.









Endnote: Challenges in phase-change problems: Interface conditions

There is a crucial difference between the classical Rayleigh-Marangoni-B6nard

problem and the evaporation problem. The difference lies in the interfacial conditions.

In the Rayleigh-Marangoni-Benard, the interfacial mass balance gives us 2 interface

conditions as a result of the no mass transfer condition at the interface. However, in the

evaporation problem, we only get one condition from the interfacial mass balance, and

consequently, we are 1 interface condition short. The interface conditions are completed

with our choice of equilibrium conditions for the temperature and pressure at the

interface. The choice of conditions at the interface between the liquid and the vapor

poses a challenge in phase-change problems. Ordinarily, the continuity of temperature is

used along with some other thermodynamic equilibrium condition (Huang and Joseph,

1992) or a relation derived from the kinetic theory (Oron et al., 1997) is used instead.

Ward and Stanga (2001) have investigated the assumption of continuity of temperature

and they have verified the existence of a temperature discontinuity at the interface.

Margerit et al. (2003) studied the role of such interfacial nonequilibrium effects on

Benard-Marangoni instability. However, Ward and Stanga report that the magnitude of

the temperature discontinuity increases with the evaporation flux and note that when

equilibrium exists between the liquid and vapor phases, the temperature is the same in

each phase. Consequently, when the evaporation rate is very small we may assume that

there is thermodynamic equilibrium at the phase-change boundary and the temperature of

both fluids are equal to each other at the interface, that is

P = f(T) and T=T*.






29


It is no problem to consider non-zero evaporation in the base state. The

equilibrium relation is replaced by one that involves a phenomenological coefficient,

which may, in fact, depend upon the wavelength of the disturbance.














CHAPTER 3
INSTABILITY DUE TO EVAPORATION THE LINEAR MODEL

In this chapter we shall concentrate on the solution of the linearized equations

given in the last chapter to find the effect of active vapor layers on the evaporative

instability without and with convection. The calculations are presented first for the

phase-change problem without the complications of Rayleigh and Marangoni

instabilities. The effects of Rayleigh and Marangoni convection are added later and each

step towards the complete phase-change problem is discussed. The results for different

cases are presented in the form of plots and physical explanations are given.

The Pure Phase-Change Problem

To understand the physics of the pure phase-change problem, we first display the

results of calculations using the full set of phase-change problem equations as mentioned

in Chapter 2. The only terms left out of the equations are the terms that have gravity and

surface tension-gradient in them. The first goal of the calculations is to determine the

conditions at the onset of instability. In order to do this, we set a, the dimensionless

inverse time constant, to zero. This implies that we seek only the conditions for steady

bifurcation from the base state. Another set of calculations, where we have calculated the

values of a, was also performed for the case of fluids heated from the liquid side where

it was discovered the leading a to be always real. Upon setting a to zero, the unsealed

eigenvalue in the linearized problem is the temperature drop. The Chebyshev spectral tau

method (Johnson, 1996) is used to solve the resulting eigenvalue problem that determines

the critical temperature difference. The interested readers are referred to Appendix J for











more details on the method. The calculations are done assuming the physical properties


of water and water vapor at the saturation temperature under 1 atm, i.e., at 100 OC as

given by Rubinstein et. al. (2002). The properties used in the calculations are given in

Table 3.1.


Table 3-1. Physical properties of water and water vapor at 1000C

p = 960 p = 0.6
m m

= 2.9.10-4 kg =1.3105 kg
msec msec

k = 6.8 101 k- = 2.5 10 2 -
m sec OC m sec OC
7m2 2
K =1.7.-10 -K = 2.0.105 M
sec sec

= 6.0.104- a* = 6.0.10-3 1
C C

S= 2.3.106 J 5.8-10 2 N
kg m

y, = 2.0 *104 N
moC



To make a proper comparison we consider different physical problems and present

graphs of the ratio of the critical temperature differences for the different cases against

the wavenumber. The definitions of these ratios are shown in Table 3.2 where the

subscripts pc, Ma, g, and Ra stand for phase-change, Marangoni convection, gravity and

buoyancy-driven convection, respectively.

Table 3-2. Definitions ofR1, R2 and R3

AT
pATc-Ma
R -



Tpc-Ma
ATpc-Ma-Ra
ATpc-Ma









For example, ATpeua denotes the critical temperature difference for the case where

there is phase-change and Marangoni convection, and R' therefore is the ratio of AT Mc-a

to the critical temperature difference for the case of pure Marangoni convection, ATma -

First, we present the results for the case where the critical temperature difference

for a pure phase-change problem is compared to the critical temperature difference for a

case where the Marangoni effect is added to the phase-change problem. This comparison

is given in Table 3-3 and the readers can see for themselves that the addition of the

Marangoni effect does little to change the threshold of evaporative instability.

Table 3-3. Comparison of critical temperature differences
Wavenumber, co ATpc (C) ATpc-Ma (C)
0.1 0.03 0.03
0.2 0.47 0.47
0.3 1.78 1.78
1.1 6.70 6.69
1.2 6.96 6.95
1.3 7.25 7.24
2.1 10.64 10.63
2.2 11.24 11.23
2.3 11.89 11.88
3.1 19.21 19.19
3.2 20.44 20.42
3.3 21.75 21.73


Marangoni convection is driven by the temperature variations along the interface.

In order to understand the effect of phase-change on temperature gradients at a deflected

interface, let us first consider such an interface in the absence of gravity and observe that

the deflection must make the temperature vary along the interface. As discussed earlier,

the temperature at a trough is expected to be higher than that at a crest as it is closer to the

heat source. When perturbed, the liquid starts to evaporate into its own vapor at a trough

since the temperature at the trough is now higher than the base state temperature.











Likewise, condensation takes place at the crests. When the liquid evaporates at the

trough, the otherwise hot trough temperature goes down since part of the energy input is

now used for phase-change whereas at the crest the heat is released, and the otherwise

low crest temperature goes up. In other words this phase-change action reduces the

temperature perturbations along the interface in phase-change problems, which, in turn,

reduces the surface tension gradient-driven convection. Marangoni convection therefore

ultimately never plays a significant role on the stability of the phase-change problem.

We can also predict the reduction of the temperature perturbations at the interface by

considering the linearized thermodynamic equilibrium relation at the interface. They are

seen from the H terms in the equation 2.34 in Chapter 2, which are very small for most


fluids, thereby making the temperature perturbations very small.


7

6-

5-

4-
R
3-
R
2-

1 -

0
0.3 0.35 0.4 0.45 0.5
Dimensionless wavenumber, o


Figure 3-1. R' vs. wavenumber.

In Figure 3-1, we observe that at large wavenumbers the pure Marangoni problem

is considerably more unstable than the phase-change problem with Marangoni









convection. At low wave numbers the reverse is true. Before we go on and explain these

observations let us make clear the restrictions behind the calculations that generated

Figure 3-1.

In producing Figure 3-1 we have performed calculations by taking two similar

systems where a liquid underlies its own vapor at zero gravity, and in each system heat is

supplied from the liquid side. In one case, we pretend that there is no phase-change

across the interface between the liquid and its vapor. In this case, the flow can only be

caused by surface-tension gradients along the interface. In the other case, we allow

phase-change across the interface. To make a fair comparison, as there is initially no

flow in the pure Marangoni problem, we adjust the pressure on the vapor side in the

second case such that the phase-change rate is zero. As a result, phase-change is only

possible through perturbations at the interface. Accordingly, we have two systems with

no-flow in their base states. They are two similar systems with only one difference; one

has the possibility of phase-change upon perturbation while the other has no such

possibility at all. We then plot the ratio, R', versus wavenumber, wo, in Figure 3.1. At

small wavenumbers, not at or very near the zero wavenumber, the pure Marangoni

problem is more stable than the phase-change problem with Marangoni convection for a

water and water vapor system. Small wavenumbers imply that there are reduced

transverse thermal gradients thereby stabilizing the Marangoni effect, yet allowing the

phase change instability to continue unabated. At large wavenumbers, we observe the

strong action of thermal and momentum diffusion.

In many physical situations and laboratory experiments, gravity plays an important

role. Gravity comes equipped with the burden of buoyancy-driven convection.









However, the effects of Rayleigh convection can be avoided in applications that involve

very thin liquid and vapor layers even though gravity continues to play a role. Gravity

simply pulls the perturbed interface back to its original position. In other words, gravity

plays a stabilizing effect at all wavenumbers. Even so, its effect on the stability of the

interface is more obvious at small wavenumbers, as small wavenumbers imply weaker

surface curvature, and that explains why the critical temperature difference is higher at

small wavenumbers when gravity is present. This is depicted by the curve labeled R2 in

Figure 3-2.

Working with very thin fluid layers is difficult in a laboratory environment. Thus,

eventually we have to account for the effect of buoyancy on the problem. For that

reason, Rayleigh convection is now introduced into the problem. It is well known from

the Rayleigh problem that, in the absence of phase-change, buoyancy-driven convection

can easily be stabilized at very small and very large wavenumbers. Consequently, in

Figure 3-2 we observe that when R2 and R3 are plotted against the wavenumber, both

plots understandably merge at small wavenumbers. At large enough wavenumbers, in the

absence of buoyancy-driven convection, the stability offered by gravity is swept away by

the stabilizing effect of the phase-change mechanism and the curve, depicted by R2,

merges with the unity line. This occurs because the surface curvature and, in turn, the

phase-change becomes more important compared to gravity at large wavenumbers. In

contrast, R3 diverges to become less stable because of the destabilizing character of

buoyancy-driven convection. Yet, at even larger wavenumbers R3 ultimately merges

with the unity line, not shown in Figure 3-2.



















R




0.5-



1 2 3 4 5 6 7 8 9
Dimensionless wavenumber, mo


Figure 3-2. R2 and R3 vs. wavenumber.

In order to demonstrate the crucial role of the active vapor layer in the stability of a

bilayer phase-change problem, we vary the depth of the vapor layer for a constant liquid

depth and calculate the critical temperature difference for a phase-change problem with

Marangoni and Rayleigh convection. The results are plotted in Figure 3-3, where the

ratio of the critical temperature difference for systems with different values of 8 is

graphed against wavenumber. The ratios are defined in Table 3-4.

Table 3-4. Definitions ofR4, R5 and R6
ATcritical of a system with 6 = 0.2
ATcritical of a system with 6 = 0.1
ATcritical of a system with 6 = 0.3
R5
ATcritical of a system with 6 = 0.1
6 ATmtical of a system with 6 = 0.4
ATcritical of a system with 6 = 0.1
Figure 3-3 shows us that increasing the vapor depth stabilizes the interface. The

depth of the vapor layer plays a dual role when we account for the Rayleigh convection in

our problem. The deeper the vapor layer, the more unstable it becomes due to Rayleigh

convection which is scaled as the cube of the depth; however, the important stabilizing









role of the increased vapor depth on the interfacial stability due to phase-change is still

present. The stabilization comes about because of the increase in the vapor flow;

meanwhile the destabilizing effect of heat transfer is also weakened. To understand the

stability mechanism in Figure 3-3, let us first look at the problem from a heat transfer

point of view.

The heat that is transported to the interface is used for phase-change and kinetic

energy with the balance being conducted away from the interface through the vapor. In

our model, the interface temperature in the base state is always equal to the saturation

temperature. Its value depends on the pressure at the top plate and the top plate pressure

is set to make the evaporation rate in the base state equal to zero. Hence, a change in

fluid depths must influence the pressure required at the top plate to keep the phase-

change rate equal to zero and thence the saturation temperature at the interface.

Therefore, the proximity of the top plate to the interface affects the temperature gradients

at the interface in both phases. The further away the top plate, or in other words, the

greater the ratio of the vapor depth to the liquid depth, the smaller the temperature

gradients in the base state. We learned from our preliminary analysis of the phase-

change problem in the absence of fluid dynamics that the heat transfer mechanism

destabilizes the interface. Thus, from a heat transfer point of view, the deeper the vapor

depth the less the destabilizing effect of the heat transfer. From a fluid dynamic point of

view, the deeper the vapor layer the less the top plate impedes the stabilizing flow in the

vapor. Hence, the interface becomes more stable for these two reasons. It is this

stabilizing role of the vapor flow that makes the passive vapor assumption questionable.














20 -
R6
15-
R
R'
10 -






1 2 3 4 5 6 7 8 9
Dimensionless wavenumber, o


Figure 3-3. R4 through R6 vs. wavenumber.

Stabilizing Effect of the Vapor Flow

In order to understand how the vapor flow stabilizes the system, it is best to look at

the velocity profiles at the onset of instability of a pure phase-change problem in the

absence of gravity and Marangoni effects.

In Figure 3-4 the velocity profiles are shown for a case where the vapor depth to the

liquid depth ratio is taken to be equal to 0.4. The position, z0 = 0, represents the liquid-

vapor interface in the base state. To determine the direction of the flow at a trough (or at

a crest) we calculate the value of the vertical component of the velocity in either phase at

the interface relative to the computed value of Z,. This is done using the eigenvector


that corresponds to the critical temperature difference. Then, we find the sign of v z (or
Z,


v'1 ) where vz1 (or vj ) is the liquid (or vapor) velocity at the interface and Z, is the
Zi

surface deflection. This ratio turns out to be negative, implying that when the deflection









is negative, the flow is in the positive direction. In other words, there is up flow at the

trough and down-flow at the crest. The fluid flows upset the temperature gradients in

both phases at the interface and more so in the vapor because of the stronger flow at the

interface in the vapor. Once the phase-change takes place upon perturbation, the flow in

the liquid will bring the warmer liquid from below attempting to make the trough hotter

whereas the stronger flow in the vapor tries to convect this heat away from the trough.

Observe the difference in the magnitudes of the velocities and the nature of the profiles

between the two phases. The horizontal component of the velocity in the vapor is very

strong compared to the liquid velocity. The vapor brings the warm fluid from the trough

to the colder crest, making its temperature go back up again. In other words, the liquid

flow convects heat to the interface and destabilizes the interface while the heat

conduction from the bottom plate to the interface also encourages instability. The vapor

flow offers stability though the heat conduction in the vapor encourages instability. The

velocity profiles of the perturbed flow tell us that the liquid flow plays a destabilizing

role whereas the vapor flow plays a stabilizing role.

To give further credibility to the statement that the vapor flow is stabilizing

whereas the liquid flow is destabilizing, we can perform a little test. By deliberately

changing the viscosity of the vapor we can strengthen or weaken the vapor flow. Thus, if

we increase the vapor viscosity by a small amount, the vapor flow should be impeded and

the system will become unstable. Indeed, an increase in the vapor viscosity destabilizes

the interface at all wavenumbers. Accordingly, an increase in the liquid viscosity

stabilizes the interface because the destabilizing liquid flow is impeded. It is useful to

remember that in both phases the flow is caused by the temperature perturbations at the







40


interface. This is obvious from the energy equation at the interface wherein one sees that

perturbed temperature gradients must result in perturbed flows at the interface. This flow

at the interface generates pressure perturbations, which in turn feedback to the

temperature through the thermodynamic equilibrium at the interface. To summarize, the

stability of a phase-change problem is determined by the competitive nature of heat

conduction and heat convection, where the relative values of physical properties of both

the liquid and the vapor affect the magnitude of the heat transfer mechanisms. The

phase-change problem in the absence of gravity becomes unstable at all wavenumbers

given that a certain critical temperature difference is applied across the fluid layers.





0.2 -

0-


-0.2 -

-0.4 -

-0.6 -

-0.8-

-1
0 1 2 3 4 5 6



Figure 3-4. The velocity vectors in both phases at the onset of instability. z = 0
represents the interface position in the base state.

Thus, if we were to conduct experiments in space, we would have to cut-off the

small wavenumbers by limiting the size of the container so that the system is not unstable

from the beginning of the experiment. However, when gravity is present, there exists a

critical wavenumber at which the interface becomes unstable, as depicted in Figure 3-5.



















AT (C)


5-




0
1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6
Dimensionless wavenumber, o

Figure 3-5. The critical temperature difference versus the dimensionless wavenumber for
water-water vapor system where d = 3mm and d* = 1.5mm.

This figure shows the critical temperature drop in 'C vs. the dimensionless

wavenumber, a). We immediately observe that the critical AT is within an

experimentally realizable range.

In this chapter, we have attempted to explain the physics of the phase-change

problem at the onset of instability in a comprehensive manner using a straightforward

linear stability calculation. We conclude that the consideration of the vapor fluid

dynamics is vital for the study of stability of the phase-change problem and leaving it out

is an invalid assumption. The linear calculations tell us what the flow and temperature

profiles are at the onset of instability. This information is enough for us to understand the

physics of the instability; however, one needs to support the theory by conducting

experiments. In order to know what to expect during the course of an experiment we

need to turn to the nonlinear analysis. This analysis will give us the information on the

type of the instability; in other words, it will tell us whether the instability is subcritical or









supercritical. In Chapter 2, we defined the input variables to be the fluid depths, the

pressure at the top plate and the overall temperature drop between the hot and cold plates.

The critical conditions for the onset of the instability are determined along with the

critical wavenumber. Advancing the temperature drop beyond the critical value will not

only reveal the magnitudes of the flow and temperature profiles, but will also give the

liquid feed rate and the vapor removal rate, as well as the enhancement in heat transport.

In order to find the answer to the behavior of the two-phase problem beyond the onset of

instability, we turn to the nonlinear problem in the next chapter.

Endnotes

The Exchange of Instability when Heated from the Vapor Side

Another set of calculations was performed for the evaporation problem when

heated from the vapor side. It is not necessary for the onset mode to be steady.

Sometimes the onset modes can also be oscillatory and one should always check the

values of the inverse time constants with the leading real part. When the bilayer system

is heated from the vapor side and only gravity effect, but not buoyancy-driven and

surface tension gradient-driven convection, is taken into account it turns out that onset of

instability is oscillatory.

The Assumption of Flat Surface in the Stability Analysis of the Phase-Change
Problem

The phase-change problem can also be studied under the assumption that the

interface does not deflect with the onset of instability, i.e., the onset of instability is

represented by the transition from the conductive state to the convective state. If the

interface were artificially constrained to be flat then surface curvature would be zero and

interfacial tension would be infinity making the product a determinate value only upon










introduction of the solved variables into the normal component of the force balance at the

interface. Under the flat interface restriction the pure phase-change problem exhibits a

minimum as seen in Figure 3-6. Weak transverse gradients stabilize on the low end, as

one needs these gradients in order to have local evaporative regions and adjoining

condensation regions. In the absence of transverse gradients we would not have local

evaporation and condensation. Both are needed to have an unstable base state. At the

high end of wave number stabilization is once again due to viscous and thermal diffusive

effects.


0.015





0.01 -


E


0.005-
b




0 0.5 1 1.5
Dimensionless wavenumber, o
Figure 3-6. Evaporation number vs. dimensionless wavenumber for two water-water
vapor systems. Curve (a) represents the results for calculations where the
interface was assumed to be flat whereas curve (b) represents calculations
with a deflected interface.














CHAPTER 4
WEAKLY NONLINEAR ANALYSIS

The Pure Phase-Change Problem

In Chapter 3, linear stability methods were used to analyze the phase-change

between a liquid and its own vapor where the vapor, in addition to the liquid, was

assumed to be fluid dynamically active. We concluded there that the active vapor layer

plays a significant role in determining the stability of the phase-change problem, and

offered explanations for the physics behind this stability. In this chapter, we will

examine what happens once we pass the onset of instability. There are several reasons to

study the post onset region. One reason is to construct the steady solution beyond the

onset of the instability and to find the magnitude of the state variables at the given

conditions as this is of experimental importance. In a laboratory experiment the onset of

instability can be detected by different means, such as IR-imaging, flow tracers or

measuring the change in heat transport with respect to a change in the control variable

once the critical state is passed. The use of tracers to observe the onset of the instability

may cause problems since it could affect the surface tension at the interface between the

liquid and the vapor. Moreover, tracer accumulation at the bottom of the liquid or at the

liquid-vapor interface may present another experimental difficulty. IR imaging, though

proven to be efficient in flow visualization of bilayer convection by Johnson &

Narayanan (1996), may be difficult because of the likelihood of the vapor condensing on

the cold top plate when the bilayer is heated from the liquid side; thereby obstructing the

imaging.









The change in the heat transfer at the bottom plate is one way to detect the onset of

instability in many convective instability problems (Chandrasekhar, 1961), provided that

the change is measurable. Therefore, another reason to study the post-onset region via

nonlinear calculations is to find the magnitude of the change in the heat transfer.

However, there are other advantages of doing the nonlinear calculations. Knowing the

magnitude of the perturbed domain variables will also tell us how the flows in both

domains change when we change any one of the input variables. The linear calculations

only reveal the velocities of both fluid phases relative to each other for a specific set of

input variables, whereas the nonlinear calculations give us the chance to compare

different bilayer systems with each other. This way, we can better understand the effect

of each input variable or physical property on the instability.

The physical model is the same as described in Figure 1.1. To get an idea on the

nonlinear features, we will investigate the phase-change problem in the absence of

gravity and Marangoni convection in this chapter.

The problem is studied in two spatial dimensions, x and z. The nonlinear equations

in dimensionless form for the pure phase-change problem are repeated below for the sake

of convenience. The momentum equations in each phase are


Pr- -+v -Vv -VP +V2v 4.1
Pr At


and


V +v Vv VP *+V2 4.2
v Pr At U


Recall that the energy equations in each domain are given by










V2T = -+v -VT 4.3


and


V2 T = T +V* NV T 4.4


Also remember that the continuity equations in each phase hold assuming that the

fluids are incompressible in the sense of the Boussinesq approximation.

We will complete the mathematical model by repeating the interfacial conditions at

the phase-change boundary. The mass balance at the interface, z = Z(x, t), is given by



p
)n= -u)-ri 4.5


At the interface, the tangential components of velocities of both fluids are equal to

each other, so v t = v* I holds.

The interfacial tension comes into the picture through the force balance. This

balance equation has normal and tangential components and is given by


Ca -V (V-P)-T1 i+2Hi=Ca -V* v -u 4.6
Pr v )p Pr P


In addition, we need the interfacial energy balance, viz.


-+K. ( -1V 2-u ( V -)--E VT.fin- VT -n


-VP .(V ) *( ) =0 4.7
-VC S(v-) -v u =









At the interface, the temperatures of both fluids are equal to each other and the

thermodynamic equilibrium at the phase-change boundary is


HKE (P-PBASE) ln T* 4.8
BASE

In the nonlinear analysis one goes beyond the critical state by advancing a control

variable such as the temperature drop across the plates. The scaled temperature drop,

however, appears as E, the Evaporation number, therefore it is E that will be advanced

from its critical value. Once E is increased all of the state variables will respond and it is

their response that interests us. Whether the response is linearly proportional to the

change or not is also of interest and will be determined. To this end we introduce

perturbation expansions about the base state for each response state variable.

Perturbation Expansions

The domain variables and their derivatives are expanded along the mapping as

follows:

I u 1 u 8 u 2u1
u=u+(A-Ac) u,+z au (A Ac) u +2 2z U1+z 0 +z
I zoj 2 1zo 2z0 zo

+-(A-Ac) u U+3z -+3z2 +3z Iz 2 +3z -- +z +z +...
16 1o az 0 zo

and

0u 0u 0fu I 8 uo 1 iB 02U Iu + 'u 0 u 8u0
S+A A +z 0 + ( A -A +2z, + +
1z azz z 0 2 az u 0 z 0 0
1 A 3 a 0U 3 2 2 2U, I 3U 0 Z2 3U u 3 0 2U 0
+-( A) +3z 2 +3z 2- +3zz 2 +3z +z + +..
6 azzo 0z 1z +z1 z 8z )

In the above expansions, z,, z, and z3 specify the mappings from the current state

in the (x,z)-domain to the reference quiescent state in the (xo,zo)-domain at first, second









and third orders, respectively. Higher derivatives with respect to zo and derivatives in

other directions follow in a like manner. It has been shown by Johns and Narayanan

(2002) that the derivatives of the spatial mappings z,, z2, z3 etc. disappear due to

algebraic cancellations and never need be of any concern. In the above expansions, A is

the dimensionless control variable, and AC is its value at which the system becomes

unstable, namely the critical value. In our problem A will be replaced by the parameter


E. We go on and set the terms equal to zero and investigate the nature of the steady-
8t

state solutions to the problem once we cross the threshold value. To do this, we will go

beyond the critical point by an amount e, i.e., A = Ac + E. The control variable is now

defined as E = E + E We will substitute the expansions of the domain variables into the

nonlinear equations and find the values of a, along with the amplitude of the

perturbations that excited the instability, using the method of dominant balance at every

order of the calculation as given by Grindrod (1996). The equations are very long and

the process of balancing the orders is tedious; thus, in order to assist the reader, we give

an overview of what we are going to do in the calculation

Overview

Suppose that the steady nonlinear problem is in the form of Nu = f, where N is a

nonlinear operator containing the control variable k, u is the solution vector and b is an

inhomogeneity vector that drives the problem. The vector u contains the velocity

components, the pressure and the temperature, in addition to the surface position, Z,

which is the last element of the vector u. The base problem must satisfy the nonlinear

equation and its solution is trivially known. In our application, the base state is a pure









conductive state, so in each fluid the velocity is zero and the temperature field is linear.

The first order problem is required to be an eigenvalue problem and this means that a

cannot be determined at this order. The first order problem is written as:

Lu, = 0.

Here, the operator L derived from the O(" ) equations by linearization of N and contains

the control variable k. The eigenvector ul is, of course, known only up to an arbitrary

multiplicative constant, A, which is the amplitude of the disturbances that cause the

instability. At first order we can determine ul and the critical value of k, but not a and

A. It is our aim to determine what A, as well as a, is through the dominant balance

method, and therefore we proceed to the next order. The second order problem is in the

form

Lu2 =fi(u1, 1)

The solvability condition is applied to this system of equations. But, in our

problem the operator L is not self-adjoint. Hence, to make the algebra simpler, the

elements of the vector u2 are solved in terms of its last element, i.e., Z2 and Z,2, where Z2

is the second order correction to the surface position Z. Solvability is then applied to the

differential equation for Z2, which is actually the normal stress balance at the liquid-vapor

interface as the differential operator for Z at every order is self-adjoint. The normal stress

balances in the first and second orders look as follows;


d2Z0
dxO









d aZ =Z (Z2



We find in this problem that, because Zi is a multiple of cos(wox), the solvability

condition is automatically satisfied at the second order; thus we move on to the third

order problem. In our problem, at third order, we find that a balance can be struck

between the various terms for a specific value of a. We then introduce the specified

value of a into the equations and apply the solvability condition on the normal stress

balance, which is now in the form of

d2- aZ3= b3(Z3,ZZ ).


By doing so, we get an equation that determines the value of A2. Hence, in the

dominant balance method, we find the value of a and A at the same order. The value of

a tells us what type of bifurcation we should expect to see once the instability sets in.

For example, if a is equal to one the bifurcation is trans-critical; and if a is equal to 1/2 it

is a pitchfork. To determine whether the pitchfork is forward or backward we focus on

the sign of A2. If it is positive the pitchfork is forward, else it is backward. In the latter

case, E must be expanded as E = E, e The work showing how to determine a is

shown in detail as an endnote to this chapter.

Perturbed Equations

Using the method outlined above, we expand the nonlinear domain equations given

earlier up to the third order. Thus far we have pretended that our physical system is

infinite in lateral extent, but the pure phase-change problem is unlike the standard Benard

problem, which produces a nonzero critical wavenumber. For that reason, in this chapter









we will choose an arbitrary wavenumber, and will proceed with the calculations using

that wavenumber. As a result, we have tacitly assumed that the fluids are confined within

a laterally finite system with slippery and insulating walls. An analysis for long

wavelength modes assuming a passive vapor has been considered by Oron et. al. (1997).

Returning to the expansions of the governing equations, we first write the

momentum equations in both phases as


(E-Ec) -VPV +V0,)+ (E-EC2 -V V PIPV2) 2 10
4.9
+ (E E )3-V P3 V3 V0v 2 2V0 1)=0


and


(E- -Ec) _* V 1 + ( E-EC 2 OP2 +V -_ -V 2.
) 2 V P v* Pr 1
4.10
+(E-EC)3a VOP3 +V2V 3 0(V 2V* VOV 0
6 P v Pr


The energy equations in each phase are


VTo +(E -Ec)a (VT -, 1 VoTo)+ 2 (E-EC ,2 2 V OTo -2v VT )
4.11
+ (E-EC )3a T3 VOTo -3 ( OT +, VoT,))=0


and


VTo +(E-Ec) VoTi -- I -VoTo

\ K K
+l(E-EC 2 T2 2- OTo* -2- V -VoT* 4.12
2 r r

+(E-Ec) T3 V3 -VOTo -3- (v2 -VOTi +V -VoT =0
6 I KI









The continuity equations become

(E -Ej .(E EC)2a 0V'V 1(E- E)3a V .3 =0
2 6


4.13


4.14


(E-E (E EC)2 V0 '.v1 (E- E)3" VO .
2 6


The unit normal and the unit tangential vectors, upon perturbation up to the third

order, become


ni= k +(E-Ec)a


8Z + l1 8EZ)( aZ azj,2
'-I- +-(E-EC) 2. -\2 1
O 0 2 xo 0 8xo 0


4.15


1 C8Z
+ (E-EC)3
6 Ox





t +(E-Ec) + I+(E-EC)2a
ax0 2 /


1
+(E
6


8ZJ 3aZ, aZj
O3 -3 O2


Z 1 I


C aZ, azaZ aZ
) xo axo Laxo


+ Z2 k
(7X0


3

lxo


At the lower plate, z


1, the boundary conditions v n = 0, -t = 0, and T = 1


are expanded as


(E-Ec) vz +-(E-Ec ) v2 + (E-Ec)3 vz3
2 6


(E-Ec) vx + (E-Ec) vx2 + (E-Ec) vx3
2 6


4.16


4.17


4.18










TO +(E-EcaT, 1 (E -E,)2aT2 + (E-E,)3a T3 = 1 4.19
2 6


while the conditions at the top plate v* n = 0, v* t = 0, and T* = 0 become

(E -E) v +1(E-Ec)v 12 + -(E-Ec v3 = 0 4.20
2 6


/ \1 2x \3 1
(E -Ec) v, +-(E-E)2a v +(E- E vx3 = 0 4.21
2 6


and


To+(E E T + (E E)2a T + (E- E)3a T = 0 4.22
2 6


For reasons of brevity, we do not provide all of the expansions of all interfacial

conditions save three that are crucial to the analysis of the problem. One of those three

conditions is the interfacial energy balance and, after simplifications, it assumes the form











-E r aT+E E ,)( k* vI-l
azo k azo c/ zo k azo

vz2+2 Z2vx
Izo c xo
a T a2T aT, aZ
2+2 Z -2
S1 E)2a B z zo 1 ax0 0
2O2 +2 Cq2 i -2
k [8z0 8 8xo 8xo

-2VPc ^zi P azi
cIzo p 8zo


aZV 8v 8v 82V
vz3-3v, 2 +3 zZ+3 c Z2 zZ
8 xa0 x zo a zo a
-3v 2 3v2, v-6 v Z, 12

+3Kv, v z +v vzi2 Vz 2


+ (E Ec)3a
6


aT, a2T, a82T,
3 3 2 Z, +3 1Z
azo 0z az0
a7T, 1z,2 T az,
axg ax0 axg ax0
-3 3

aT a2T
3 + 3 2 Z, + 3
k* 8z, 8z
k _Ta* Z2 T_
-3 3 T2
SO x xo a 3xo


33T
2 +3 ZI2

a2T,
6 Z1

0
Z2 +3
Z,

-6
axo azoc


aZ,
c5xo


8z

* 7Z,
xo cxo


6vz 2 1-2xi + +2 a21Z,
c8zo 0 80x 8zo xo Czq I

-12vx-, + 6 vz2 + 2-Z


Xl a CA 0 ^ l cZ, :
[8zo, 8xo, 8zo 8x, 8z,2

-12v, + 6 v + 2 z Z
8zo 80xo 8zo 8z2 co J


4.23









Another one is the thermodynamic equilibrium relation between the vapor pressure

and the vapor temperature;


HKE (P Psat) -In
T sat


+(E- EC))aj HKEPl


Tst1 T
Tsat


+dTo,
dzo


Zj


1
+-(E
2


1 (E E,)3a
6


HKE P2

1
T
Tsat


+2 B Z1
Oz0


aTB dT
+2 'Z,+ 0o
Ozo dzo


HKE P +3 +
np;+ 3J
T sa
+ T T+32


Tsat

2 T1 2a T1
+3. IT' <


-2T,J I
TsfaTf


2 Z, +
3 3
+
3z0


,OP
azo

3
Oz


) + Tat )IT


Z2 +3 1Z2
zo


a2T
Z +3
aZ2z


0*

'o


aT* dT
+2 Z, +
Ozo dzo


SdT
dzo


We also give the third order normal stress balance equation because we will use it

to find the values of a and A. The third order normal stress balance is given by


2
dT'
d To* 2
+ oz
dz 1


4.24


2' dTo

dT*
+ Z,
dzo










Ca(P -Po)

+(E-Efc)





+ (E-E )2
2


+ (E -EC)3a
6


av, aa2Z



P2 + 2 Z, P2 2 Z,
9zo dzo
Ca
Nf 2 0 '2 z P* "V :2* ,2V I
-2 +2 z +2 +2PzZ

Ca P 2 Z 2Z ,



P3 +3 Z, +3 Z2 + 3 Z2

-P -3 Z, 3 Z 3 2I
BzO 1 Z O z2 1Z2 0

Cva3 a2v2 v 2VZ aa3VzZ2 vza Oz1 )
Ca -2 az3 az2Z1+3 zjZ2+3 +z 1 a2a12
Oz oz2 oz2 Oz3 I Oz Ox

+2 3 +3 2 Z, +3 a Z +3 Z2 +12
S ZO9 o Bz2 Z2 BZ3 1 BzO Ox


4.25


Ca avi az), .
+3 vz2+2 Zl -2vxl (v, -v0)
Pr azo xo)

Ca av az v z2 OZ,
+3 vz +2 Z 2vx, v2 -2 Z, + 2vix v
Pr azo axo azo I xo

+2Z3 2Z aZ)2
+ 09
ax2 axo2 lax


Remember that in all of the expansion formulas, E will be replaced by Ec+e. The

first order problem is an eigenvalue problem and it was solved in an earlier chapter. The

solution to the first order problem is used at every subsequent order along with the

solvability condition. We learn from the first order equations that the eigenvector ul in

terms of Z1 is









U, = A i, (Z,)cos xo,.

Recall that o) is an input and it is determined by the width of the compartment in

the lateral direction. We cannot determine the value of a either from the first or the

second order equations. As noted in the earlier section, the solvability condition is

automatically satisfied at the second order. The solution to the second order, U2, will go

as

u2 =B u2a (2)coS(oX+A2u2b(Z12)+A2u2c(ZI2)cos2mxo.


It is at the third order that we learn that a is equal to 2. At this order, there exists


a value for A other than zero. A appears as A2, and the value of A2 is obtained using the

solvability condition on the third order normal stress balance. The third order solution,

U3, depends upon the form of ul, u2p as well as Z3, i.e.,

u3 =C u3a (Z3) COS O)X + A3u3b (ZI3, ZZ2)COS 0)X
+A3U3c(Z3,ZZ2)cos3mxo +Au3d(Z,)cosxO +B Au3e(Z,Z2) (cos 20xo +1)


Results of Calculations

The Chebyshev spectral tau method is used to solve the eigenvalue problem that

determines the critical temperature difference and the value of the unknown

multiplicative constant A. Once we know what A is, we can calculate the actual change

in the heat transfer rate at the bottom plate when we advance AT a little beyond the

critical point. Our system is finite in the lateral extent with periodic conditions. The heat

transfer up to the second order is calculated as










2 k rdT, dT 1 dT2 dT2
q-n= -k +E2- Cosx0+-E +- cos2mxo dxo 4.26
0 0 dz, dz 0 2 dz+ dzo


Clearly, coswex,0 and cos 2wx,, when integrated over the width, give zero. We are

more interested in the percentage of change of the heat transfer, so we have plotted


0 4- against in Figure 4-1 for a system where the liquid and vapor depths are
qo *no EC

4mm and 0.4mm, and the wavenumber is co = 0.3. In nonlinear calculations the

wavenumber of the disturbance is fixed. Accordingly, one can find the width of the box

using the wavelength of the disturbance. The wavenumber, co = 0.3, corresponds to a

cell width of about 84mm.

We learn that once the interface is unstable, more heat is transferred through the

bottom plate to the system. Thus, if we wanted to know when the system becomes

unstable, we only have to measure the amount of heat that is needed to keep the bottom

plate at a constant temperature. The heat necessary changes with the onset of instability.

For a typical problem we can observe a 2% increase for a supercritical change of about

10%.



















-S 1






0


-0.5
0.98 1 1.02 1.04 1.06 1.08 1.1
E/E


Figure 4-1. Heat transfer across the bottom plate.




The nonlinear analysis results do more than giving the heat transfer enhancement.

They can also be of use in helping us to better understand the physics of evaporative

instability. To see how, let us first observe from the linear calculations that the deeper

vapor layers make the interface more stable. Then apply the nonlinear calculations to

two cases with the same base state interface pressure and temperature where in one of the

systems 3 = 0.1 and in the other 3 = 0.3, and compare the magnitudes of the pressure


perturbations in the vapor at the interface. Because P = P= + E2P* at the interface and

because we wish to compare the problems with each other for the same percentage of


departure from the critical state, we will express P as



PP E -E P1,
Q ,E
^C











In Figures 4-2 and 4-3, we have plotted the magnitudes of pressure perturbations in


the vapor, EY2 P,, for the cases where 3 = 0.1 and 3 = 0.3. We learn from the


calculations that even though the nature of the profiles of the vapor pressure perturbations


are the same for both cases, the magnitudes are not. The vapor pressure perturbation at


the interface for the case with the shallow vapor layer (3 = 0.1) is about three times more


than the pressure perturbation at the interface for the case with 3 = 0.3, resulting in a


more unstable configuration. In Chapter 3, during our discussion of the physics of the


pure phase-change problem we had deliberately changed the value of the vapor viscosity


to see its effect on the instability. The system, where the vapor viscosity was increased,


proved to be more unstable.



0.1

0.09

0.08

0.07

0.06

Z 0.05 -

0.04

0.03

0.02

0.01 -

0
-7.174 -7.172 -7.17 -7.168 -7.166 -7.164 -7.162
E /2 P 10-3


Figure 4-2. The scaled pressure perturbations in the vapor domain for the case where
d = 1mm, 3 = 0.1 and ) = 1. z = 0 represents the interface position in the
base state.















0.25


0.2


Z 0.15


0.1


0.05


-2.23 -2.22 -2.21 -2.2
E 1P x10-3


Figure 4-3. The scaled pressure perturbations in the vapor domain for the case where
d = Imm, 3 = 0.3 and c) = 1. z = 0 represents the interface position in the
base state.

As shown in Figure 4-4, our nonlinear calculations return a larger pressure

perturbation value at the interface for the system with the more viscous vapor. Thus, the

nonlinear calculations support the argument that the vapor flow dynamics affect the

stability of the phase-change problem.














0.25


0.2


Z 0.15


0.1


0.05


-3.28 -3.27 -3.26 -3.25 -3.24 -3.23
E1 P. x103

Figure 4-4. The scaled pressure perturbations in the more viscous vapor for the case
where d = Imm, 3 = 0.3 and o) = 1. z = 0 represents the interface position in
the base state.

Endnote: The Determination of a

There is a lot of analytical work that goes into determining the value of a For

that reason, we will demonstrate this cumbersome procedure on a simpler problem.

Assume the following domain equation

d2
d + Af(u) =0 4.27
dx


df d2f
where f(0) = 0, f' (0) (0)> 0 and f" (0)= 2 (0) > 0 are given, and the domain
dx dx

equation is subjected to Dirichlet condition at x = 0 and x = 1. To determine the

behavior of the system around the critical point, in this case Ac, we will go a little beyond


the critical point and replace A by Ac + When we expand the domain variable u


around a base state uo as a power series of "a, the domain equation can be written as









dx22 au1 + I2a,2 + 3aU3
dx 2 6

+A f (O) Eu, +-E 2aU +-E 3a 31 +-1f(0) u +-E 2au, + 3a U3 ..
2 6 2 2 6

+e f' (0)EU +- 2aU2 +2 3a U3 +-f"(0) Eau, +-E2a u2 + 3a +... = 0
2 6 2 2 6

The E" order problem is then

d2
U+ f (0)u, =0 4.28
dx2

subject to u, (0) = 0 and u, (0) = 0. When we solve for u1, we find that u, = A sin(rx)

2
and A, = Now we will solve for u2 and u3 and so on. To do that we will move to
f' (0) 2

the next order of equations, which is unknown since a is yet unknown and needs to be

determined. There exist 3 possibilities:

1+a < 2a is one of the possibilities, and this implies that

f' (O)u1 =0.

Now, this can only be true if only if u, = 0, so obviously this is not a plausible

choice. Another possibility is that 1 + a = 2a. When we set a = 1, we get the following

equation as the next order equation:

1 2 1 1 (0)f"(0)u2 '(0)u1=0 4.29
+ Pyf'(O)u2+ 1yf"
2 dx 2 2

If we multiply equation 4.28 by u2 and equation 4.29 by ul, then integrate both

equations over the domain, and subtract them from each other we get










1f "(0) jfU3 dx+f'(O)Jf u' dx= 0
0 0

Using the relation above, we can determine a value for A which is non-zero, and it

is


A 3 f (f))2
4;r f"(0)

So far, one of the three possibilities gives us a non-zero value for A. The next

possibility is that 1+ a > 2a. Now this gives us

1Af (o0) u dx = 0




This is not possible since ju' dx = Thus, we conclude that a = 1. Note that
f 3;r

the value of A is determined at the same order where we determined the value of a .

A different approach to the weak nonlinear analysis

There is more information one can actually get from the weak nonlinear analysis. To do

that, we will have to modify our physical model. This time we will again set the

evaporation rate in the base state equal to zero by adjusting the vapor pressure, however,

now we will leave the feed and removal valves open, and as we increase the temperature

difference applied across the bilayer system we will adjust the vapor pressure accordingly

such that the evaporation rate is equal to zero in the base state even if we cross the critical

point. Intuitively, one expects the net evaporation rate to change once we cross the

critical point. Beyond the critical point the net evaporation rate, U, can be expanded as

U = Ubase + Uextra 4.30










where Uextra is


Uetra = Uflat +3 4.31


In equation 4.31, Uflat is the change in the net evaporation rate if the instability had not

occurred, and 3 is the correction to it because of the onset of instability. In this analysis,

we will adjust the vapor pressure as we increase the temperature difference applied,

accordingly, Ulat is always zero even after we cross the critical point. Thus, the only

new unknown in the problem is 3. The sign of 3 tells us what will happen to the net

evaporation rate beyond the critical point. If it turns out to be positive, then with the

onset of instability the system will suck in liquid through the bottom valve and remove

the excess vapor through the top valve. One does need to solve the complete weak

nonlinear problem to get this information. Moreover, one can even determine the sign 3

by analytical calculations. From our earlier analysis we know that the solvability

condition was automatically satisfied at the second order. That means that we can go

ahead and solve the second order problem. By doing that we will find 3 in terms of the

amplitude of the deflection. The momentum equations in each phase at the second order

are

2
-V P +V2 = v -V V v 4.32
Pr


and


AV 2 v 2 V o 4.33
0V p" +V'-Vl V0V 4.33
'u v* Pr


The energy equations in each phase are







66


V2T2 2 VOTo = 2v VoT,


V2 K OC *
0 2 2 0 V0O


The continuity equations become


Vo V2


Vo v =0


At the lower plate, z


1, the boundary conditions v- i


U, v-t=0, and


T = 1 are expanded as


Vz2 = U2


vx2=0


T, = 0


while the conditions at the top plate V Ui = U*, V


0, and T = 0 become


V 2 = U


Vx2 =


4.34


2 K V* T


4.35


4.36


4.37









For reasons of brevity, we do not provide all of the expansions of all interfacial

conditions save three that are crucial to the analysis of the problem. One of those three

conditions is the interfacial energy balance and, after simplifications, it assumes the form

vza2 E 2 -2 zI Z, + 2vx Z
azo k zo H azo Cx0
r1 ^ j 4.38
az2T T, 8Z, kx2 C9 2Tv2 ,, ST ZV, 'v
-E 2 Z 22 1 1 C+2
z k 0 xo xo 0zo 0 zo


Another one is the thermodynamic equilibrium relation between the vapor

pressure and the vapor temperature;

0P(PI 1 + dT'
HT P2+2 Z2 0Z2
K2 azo') Tbase dzo
2 2 4.39
1 T LT1 dT1
2 L Z, -T + 1 Z
base 0 Tba dzbs

The normal stress balance assumes the following form


Ca (2 -P2 z2 +2 z2 2


=Ca 2 Z, -2 Z,+4 Z -4 Z 2 vz vzI v- v
Szo azo az 1 p 1z Pr

When one solves the first order problem, one determines that the domain

variables are actually trigonometric functions in the x0 -direction. They are, in the liquid,


Vx1 = x1 (Z0)sin Ox0

Vz1 = i1 (Zo)COS Xo

T, = T (zo)cosCxo









P, = P (zo)cos x0

The vapor domain variables assume similar forms. We plug them into the second order

domain equation and find the form of the second order domain variables. What follows

is the demonstration of how we find the forms of the second order domain variables. The

second order energy equation is

VT2 -V VOT = 2v, VoT,

After substituting the first order domain variables we get

V2T dT aTdT
V0 22 2vx, + vz1
dzo OXo OZo

which translates to

S dT1 2 C2
V T2, -z V dz =2 -o, (zo)iT (zo)sin2 O +X z dZo o) COS2XO


Using trigonometric manipulations we get

V2T2 v 2 dT0 K,. dT1 d ( 1
dzo dzo dzo

+ (zo) (zo) -t (zo) (zo) cos2wxo

From the manipulations above, one learns that

T2 = (z )+ (+zT )cos2wxo

Following the same method, one obtains

vx2 = x2 (ZO)+ vx2 (Z ) sin 2wxo

Vz2 = Vz2 (ZO) + Vz2 (Z ) COs 2xo











P2 = (ZZo)+P2 (Zo) COs 2cxo

We are trying to find the sign of 3. At the second order, 3 is actually equal to the x

independent part of the correction to the net evaporation U, which is U2 Thus, we only


need to solve for the x-independent part of the particular solution for U2. The domain

equations in the liquid for x-independent part of the second order problem are

d2
dzo


dP2 d2 z2 2 d, i
+- =-v
dzo dzg Pr zdz


dT 2 dTo = d, di,z
--- v 2 = v i -- + Ti --
dz0 dz0 z dzo dzo


and


6z2 =o
dzo

In the vapor domain, they are

d2 *
dzo


SdiP d2 z2 2 v dvv,
S+ -- --V
pt dz0 dz2 Pr v dzo


dT_ K, dT_ K dT,' dv
v V z2 -- Vz 1
dzo K dzo K dzo dzo


and



dzo







70


At the bottom plate, the boundary conditions are


vz2 = U2

T =0


At the top plate,


V* 0
x2=


z2 2

1 =0
2


At the interface, we have


Vz2 = -Vz2



Vx2 = Vx2


dTz 2
dz 0 2 '2


dTz
dz0 )


dzx2 dxi 2
dzo t dzo


1 j2* dT Od-
2T dz 2 KE HP2

1 dT I 1 1
Tbase dzo0 Tbase


4.45


+ d Z, + IKE d, 21
dzo dzo


4.41



4.42



4.43


dTi1
dzo


dTl 1
dzo


4.44


(2









dTi, k* di di, / v. z d6i
-E --- = E zi 1Z+2V, 1-- iz 4.46
dzo k dzo dzo K v dzo


S di diP z 2 Ca p p
Ca P -P -2 z2+2 +_=+ 1-- 447
2 dz, p dz, 0x; Pr p p)




These equations allow an analytical solution for the problem. First, we learn that

everywhere in each domain




Uv2(zo)=

UZ2 (ZO)= f


v*z2 P U


We also learn that





p* 1 ,2
p Prv

For reasons of brevity, I will omit giving the analytical solution for T2 (z,) and T2 (z,) .

In conclusion, we find the correction to the net evaporation, U,, and it turns out

to be positive at each wavenumber. This information is useful if we were to do

experiments since now we have an additional indicator, other than measuring the heat

transfer at the bottom wall, to determine whether the instability has set in or not.

Opening the valves and allowing the bilayer system such in liquid through the bottom






72


wall with the onset of instability is a preferable method since it is a visual method, unlike

measuring the change in the heat transfer.














CHAPTER 5
THE RAYLEIGH-TAYLOR PROBLEM WITH PHASE-CHANGE

Introduction

As we discussed earlier in Chapter 1, a vapor underlying a liquid results in a

configuration that is potentially unstable. This problem, known as the Rayleigh-Taylor

problem, is well known and an account of its underlying physics is given by

Chandrasekhar (1961). It may be recalled from Chandrasekhar that the interface is

unstable for wavenumbers smaller than a dimensionless critical wavenumber given by


WRT = (I g p) gdj-2


Here, p' and pg are the liquid and vapor densities, respectively, y is the interfacial

tension, d is the horizontal length scale and g is the coefficient of gravitational

acceleration. The derivation of this result is also given in Appendix D. The classical

Rayleigh-Taylor instability problem is ordinarily studied under the assumption that there

is no mass transfer across the interface. The problem was reformulated by Hsieh (1972)

who considered the Rayleigh-Taylor problem with mass and heat transfer in inviscid

fluids of finite vertical depths. He solved both the velocity and the temperature fields in

both fluids and derived an analytical dispersion relation i.e., a relation between a and 0c .

Later, Hsieh (1978) studied the same problem again for inviscid fluids but introduced a

phenomenological coefficient, a, which allowed him to decouple the velocity and

temperature fields for both infinitely deep fluids as well as for fluid layers with finite









depths. Using the same idea of a phenomenological coefficient, Ho (1980) analytically

obtained a result for viscous fluids with finite depths. In the same paper, Ho also derived

the critical dispersion relation for viscous fluids with equal kinematic viscosities for

infinitely deep fluid layers. As an application to pool boiling, Badratinova et al. (1996)

have found an analytical expression for the critical wavenumber for a system where a

vapor with finite depth underlies its own infinitely deep liquid, solving for both the

velocity and the temperature fields. However, we have seen in Chapter 3 that the

assumption of large depth fluids may lead to conclusions that do not quite represent the

physics. We might also like to know what the phenomenological coefficient used by Ho

(1980) and other workers such as Adham-Khodaparast et al. (1995) really means. These

workers found a critical dispersion relation for infinitely deep fluid layers with equal


kinematic viscosities to be given by w2 +4- w w = 0. Here, v is the cinematic


viscosity of both the liquid and the gas phases. Note that upon setting the

phenomenological coefficient equal to zero in the above relation, one recovers the

classical Rayleigh-Taylor instability result.

In this chapter we therefore investigate the Rayleigh-Taylor problem with

interfacial mass and heat transfer and present the results of our calculations for viscous

fluids where both fluid layers have finite depths. In addition, we derive an analytical

expression for the phenomenological coefficient, a, albeit for infinitely deep viscous

fluid layers and show that this coefficient is actually a function of the wavenumber.

Physical and Mathematical Model

Unlike the earlier chapters this time the physical model is a system where a vapor

of depth, dg, underlies its own liquid of depth, d', and is heated from below (Figure 5-1).









COL4XD PUATE





Vapor

HWPLATE M
Figure 5-1. The physical model

As before the base state is one where there is a flat interface between the liquid and

its underlying vapor, and the stability of this base state is in question. Also as before the

liquid and vapor depths are assumed constant by suitably adjusting the vapor removal and

liquid feed rates. The phase-change rate at the interface can be controlled by adjusting

the vapor pressure at the vapor-side wall. The phase-change rate is not constant as it can

change upon perturbation. Consequently, in this problem as in the pure evaporation

problem the input variables are the liquid and vapor depths, the temperature difference

between the bottom and top plates, and the pressure at the lower plate. Furthermore, the

pressure at the lower plate can be adjusted such that the phase-change rate at the interface

is equal to zero. In the analysis of the problem, we assume as usual that the liquid and its

vapor are in thermodynamic equilibrium because we will only consider either zero or

very small phase-change rates at the interface in the base state.

The equations that model the physics are given by the Navier-Stokes, energy and

continuity equations in each phase assuming that the fluids are incompressible. It is

assumed that the bottom and top plate temperatures are kept constant. When we assume

that the liquid and its vapor are in equilibrium and that no phase change is taking place

across the interface, we need not feed and remove fluid through the bottom and top

plates. The no-slip condition applies along the top and bottom plates. Hereon, vx, vz,









and T are taken to be the x and z-components of velocity, and temperature fields,

respectively, while the asterisk denotes the vapor phase.

The equations are brought into dimensionless form using the following scales:

U/ K d
length, d = ; velocity, v = ; time, where ul and K are the liquid viscosity
y d v

T
and the liquid thermal diffusivity, respectively. The temperature is scaled as T =
dTog d
dzo

dT'
In the temperature scale, 0 is base state temperature gradient in the vapor at the
dzo

interface. Linearized stability is performed in the standard way by expanding the state

variables in the form of u = z)e'o e"o leading to the following dimensionless

equations of motion

d2 WK d2 w2 av1 0 5.1
dz dz Prz


and


K w2 -w2J)d2 w2 V 1 = 0 5.2
dz2 dz2 Priv zv


where a and Pr are the inverse time constant and liquid Prandtl number, respectively, v

is the ratio of the vapor kinematic viscosity to the liquid kinematic viscosity and w is the

dimensionless wavenumber of the perturbation.

The perturbed energy equations in dimensionless form become


S -w2 -- T1 kv 5.3
2 dz C









and

wd2 2 g 54
-2 -- Tzi= vz 5.4
dz / ) /


In the equations above, k is the ratio of the vapor thermal conductivity to the liquid

thermal conductivity, and K here is the ratio of the vapor thermal diffusivity to the liquid

thermal diffusivity. At the lower plate, the no-flow, no-slip and constant temperature

dvp
imply that the vapor is subject to vg = 0 dz = 0 and Ti, = 0 while similar conditions
dzo


dv1
at the top plate become v1 = 0, dz = 0 and T, = 0 .
dzo

At the interface the mass balance is given by

vy, -p V = C (1- p) Z, 5.5


where Z, is the perturbed surface deflection, and p is the ratio of vapor density to liquid

density.

The continuity of temperature and the no-slip conditions at the interface become


T,1 +k Z, = T1 +Z1 and dvz dv1z
dzo dzo

The normal and tangential stress balances assume the following dimensionless

form:


d vz 3w2 + zl z1 3w2 + z1 2 w 2)Z, = 0 5.6
dz Pr v dzo dz' Pr dz RT










S +w2 v1 Kil +1w2 V1 = 0 5.7


where p is the ratio of the vapor viscosity to the liquid viscosity.

The interfacial energy balance is

vI dT1 k dTi>
vz -a Z, +E k = 0 5.8
dzo dz0


C' AT
where again E stands for the Evaporation number and is given by E = Here, C'
h P

and h are the liquid heat capacity and the latent heat of evaporation per unit mass,

respectively. The equilibrium condition at the interface becomes

d3V+ 1> dv1>
HKE 2 + 7 z1 2 1PE Z, = W21 (T" + Z,) 5.9
E dzj Pr v dz0


V9K1 gd
where again KE 2 and PE = are two dimensionless parameters from the
h d h

linearized Clapeyron equation, and 7 is the inverse of the scaled saturation temperature

at the flat interface.

Analytical Results for Infinitely Deep Fluids

Ordinarily the fluid layers are assumed to be of finite depths but to provide the

reader with an expression for the phenomenological coefficient, a by way of an

analytical procedure and avoid a numerical calculation we assume that they are infinitely

deep.

The Rayleigh-Taylor problem is essentially an interfacial instability problem, so to

find the leading inverse time constant, o, one can omit the 'o's coming from the domain

equations and domain variables, and keep the 'o's that appear only along with the surface









deflection. Doing this makes it easier to find an analytical dispersion relation for a .

When we subject the vapor and liquid domain equations to the conditions at z, = -0o and

zo = oc, respectively, we get the following expressions for the velocity and temperature

fields.

v, = A e-wzo +B1 z, e"-wzo


v1 = A ewzo +Bg Zo eWZo



T = C1 eWZO 2- A + zo e-wzo +-Bl z ew
12w 2w) 2



Tig = Cg ewzo +-- A ) zo ewzo +Bgz2 ewzo
2w ic 2wK w 2 0


The steps taken to get the analytical result are given below. Using the no slip

condition and the tangential stress balance one gets

-A'w + B = Ag'w+Bg

and

Alw-B' =1 (Aw+Bg).


We deduce from these two equations that

Alw = B1


and


Aw = -Bg.









Now we can write two of the unknowns, B' and B1, in terms of A' and of A'. We

then plug in the expressions for v1, and vg, into the normal stress balance and we get

2A'w + 2uAg'w = (w, -w2 )Z,.


The mass balance gives us

A' = pA +(1- p))cZ,.


We can use the relation obtained from the interfacial mass balance to eliminate A'

from the normal stress balance and we get Ag in terms of the first order perturbation of

the surface deflection, Z1, which is

Ag --w2 2w(1- p)
2w(P+p)


When we plug the expressions for T,1 and T,1 into the continuity of temperature

and the thermodynamic equilibrium relation we get

C9 + Z, = C' + k Zi


and

17(C + Z,) = -2IKE wA 17PE .


Now we are in a position to write C' and Cg in terms of Z1. We still have not used

the interfacial energy balance. When we plug everything into the energy balance at the

interface we get

(....)Z, =0.


For this equality to hold either Zi or the term in the parenthesis must be equal to

zero. Zi equal to zero is not an option since then there would be no solution to the









problem solved. Setting the terms in the parenthesis equal to zero we get the following

expression as the dispersion relation:

F+G
a=- 5.10
H

where


2 2 W2 '(~ ) 3k (p+ 1
w w-2) RT E +(1+ k) I
2w(/+p) 77 4w




G = wE HPE (1+ k) + 2k


and



3k (1- p) r KE (1+ k) 3k(
H= p+ 3E(1-p)+ (IUP)p+E I
4w (P+p) 7 4w



When we set a = 0, we get an expression for the critical wavenumber in terms of

the vapor temperature gradient at the interface and the physical properties of both fluids.

We mentioned earlier that Ho (1980) and Adham-Khodaparast et al. (1995) derived

a dispersion relation for the phenomenological coefficient, a, by decoupling the

temperature and velocity fields. Once a is introduced only the velocity fields in both

fluids were obtained. To obtain the result of Adham-Khodaparast et al., one only needs

to drop all the energy equations from our model and replace the interfacial energy

balance by the following balance equation in dimensional form:

vli -UZ, = aZ,









Simply put, the phenomenological coefficient is introduced in the interfacial energy

balance where it replaces the heat conduction terms. The dispersion relation of Adham-

Khodaparast et al. using dimensionless equations becomes


w +2av +v ')d w-w =0


When we assume that the kinematic viscosities of both fluids are equal to each

other, we get the familiar result of Adham-Khodaparast et al. in dimensionless form. It is

2 av d 2
w +4- w-wRT = 0


Going back to analytical dispersion relation, when we rearrange the equation for

the critical wavenumber, we get the expression for the phenomenological coefficient, a,

at the onset of instability. It is


a = d2 (G

S[nKE2(1+k) 3k p+
77 4w



We see that a, at the onset of instability, itself is a function of the wavenumber as

well as the physical properties of the fluids.

In Figure 5-2, the dimensionless phenomenological coefficient, a, is plotted

against the dimensionless wavenumber to see the shape of the curve and investigate

whether there is any range of w in which ac remains roughly unchanged. The

calculations are done assuming the physical properties of water and water vapor at its











saturation temperature under 1 atm, i.e., at 1000C. The physical properties used in the


calculations are given in Table 3-1.


Figure 5-2 shows the strong dependence of the phenomenological coefficient on the


wavenumber in the plotted wavenumber range. The classical Rayleigh-Taylor


dimensionless critical wavenumber for the same system is 3.4e We observe that at


exactly the classical Rayleigh-Taylor critical wavenumber the value of the


phenomenological coefficient is equal to zero. When we look at the relation derived by


Adham-Khodaparast et al., it is easier to see why a = 0 when w = wRT.


x 106





a
5-

4-
0 4
0


2


0 \
0 05 1 15 2 25 3 35
Dimensionless wavenumber, x 10-


Figure 5-2. The dependence of the phenomenological coefficient on the wavenumber

The graph therefore tells us that the use of a phenomenological coefficient to


decouple the temperature and velocity fields is not a good assumption as the


phenomenological coefficient itself depends strongly upon the wavenumber. In fact, the


coefficient even turns negative for large wave numbers.









Numerical Results for Finite Depth Fluids

To see whether the assumption that one or both fluids are infinitely deep is a valid

assumption, we have plotted the critical temperature gradient in the vapor against the

wavenumber in Figure 5-3 for different systems. In one of the cases, the fluids are

assumed to be infinitely deep whereas in the others both fluid layers have finite depths.

Positive temperature gradients mean that the bilayer system is heated from the liquid side,

namely from the top.

We know from the classical Rayleigh-Taylor problem that the system is unstable

for wavenumbers smaller than the critical wavenumber, however, when there is phase-

change taking place at the liquid-vapor interface, the interface can be stabilized by

heating from the liquid side even at small wavenumbers. In the plots, the region under

the critical curve is the unstable region. As we increase the liquid and vapor depths, the

critical values of the vapor temperature gradients approach the values for the infinitely

deep layers. We observe that for large wavenumbers the curves overlap whereas for

small wavenumbers the critical temperature gradients in the vapor at the interface diverge

from the values for the infinitely deep layers. From Figure 5-3 we can tell that the

assumption of infinitely deep layers is also not a good one for all wavenumbers. For a

water-water vapor system, the finite layers, where d' = d9 = 0.1 m. deep, give the same

results as the infinitely deep layers at all wavenumbers, however, a 0.1-meter depth for a

vapor layer is a relatively large fluid depth under many practical circumstances.













09-
08 -
07
06-

a05
S04-

03 c
02 -

5 01 -b
a
05 1 15 2 25 3 35
Dimensionless wavenumber, x 107


Figure 5-3. Comparison of models a) Infinitely deep fluid layers
b) d'= 0.05 m and d = 0.05 m c) d'= 0.01 m and d = 0.01 m

In many applications of the Rayleigh-Taylor instability, the vapor layer is

shallower than the liquid layer. For that reason, in Figure 5-4 we plotted the

dimensionless critical temperature gradient against the wavenumber for two different

cases to see the effect of the vapor layer depth on the instability. Observe that when

realistic fluid depth ratios are considered the effect of the fluid depths is considerably

more important. The larger the ratio of the liquid depth to the vapor depth is, the easier it

is to stabilize the bilayer system by heating from the liquid side. That is, the minimum

temperature gradient required to stabilize the system is less. For shallow vapor layers,

the flow in the vapor is impeded by the bottom wall and the flow in the liquid starts to

play a role on the stability. The liquid flow brings warmer liquid to the troughs

encouraging evaporation, thus preventing them from growing deeper. More thorough

explanations of the stabilization mechanism through the flow generated by phase-change

are given in chapter 3. Thus far, in all the cases we considered, the instability was

delayed by heating from the liquid side. Again, when we turn to industrial applications,









we recognize that the system is mostly heated from the vapor side, namely heated from

below.

Figure 5-5 gives a plot of the critical temperature difference across the bilayer

system against the wavenumber where both of the fluids have finite depths and the

system is heated from below, namely the vapor side. The liquid and the vapor layers are

taken to be Imm and 0.1mm deep, respectively. The region above the curve is the stable

region. So, this calculation also pertains to a physically realizable experiment. Observe

that the curve shows an increase of the critical temperature difference with wavenumber

at low wavenumbers, a maximum and then a decrease. This shows the dual effect of

wavenumber. The waviness of the disturbance finds itself participating in transverse

diffusion of heat and momentum as well as in stabilization on account of surface tension.

At low wave numbers, surface tension plays a minimal role. Transverse diffusion itself

plays two roles. In the domain it offers a dissipative effect and therefore plays a

stabilizing role but at the interface it encourages the movement of the interface as hot

locations move toward cold regions and so forth. This is the destabilizing effect. At high

wave numbers the stabilizing effect of surface tension and the critical temperature

difference for the realization of stability is decreased.














09-

08 -

07-

06

a05


03

S02

S01

05 1 15 2 25 3 35
Dimensionless wavenumber, c x 10


Figure 5-4. Effect of shallow fluid layers.
a) d'= 50e-3m and d = 50e-3m b) d'= 50e-3m and d9= 5e-3m

Note that heating from below can stabilize an arrangement where a vapor underlies


its own liquid at a wavenumber at which the classical Rayleigh-Taylor analysis indicates


otherwise.


In summary, we conclude in this chapter that the assumption of infinitely deep fluid


layers give results that are very different than for the fluid layers with finite depths as one


can conclude that a system that is unstable for a given temperature gradient in the vapor


based on the infinitely deep layers assumption may actually give stability for fluid layers


with finite depths.