UFDC Home  myUFDC Home  Help 



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 VuQuoc, also deserve my gratitude. I have really enjoyed taking classes from Prof VuQuoc 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 parentsinlaw, 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 VaporLiquid Phase Change................ Evaporative Instability......................... ............. ... ....1 Rayleigh Convection or BuoyancyDriven Convection........................6 Pure Marangoni Convection or Interfacial Tension GradientDriven C onvection ................................... .................. ...... ...8 R ayleighT 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 phasechange problems: Interface conditions ....................28 3 INSTABILITY DUE TO EVAPORATION THE LINEAR MODEL....................30 The Pure PhaseChange 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 PhaseChange Problem ................. ...............................42 v 4 WEAKLY NONLINEAR ANALYSIS ........................................................ 44 The Pure PhaseChange 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 RAYLEIGHTAYLOR PROBLEM WITH PHASECHANGE...... .....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 RAYLEIGHMARANGONI 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 PhaseChange 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 31 Physical properties of water and water vapor at 1000C .......... ..... .... ........ 31 32 Definitions of R R2 and R3 .................. ............................. ......31 33 Comparison of critical temperature differences ........................................... 32 34 Definitions of R4, R5 and R6................................................................ .....................36 61 Physical properties of silicone oil and air at 200C ...................................... 96 LIST OF FIGURES Figure page 11 Physical model of the phasechange problem....... ...............................2... 12a H eat input from the liquid side....................................................3 12b Heat input from the vapor side................................ ........ .. .. .................4 13 Physics with fluid dynamics..............................................5 14 Rayleigh convection ................... ........... .. ........ ........... .... ....... .6 15 M arangoni convection..............................................8 31 R' vs. wavenumber. ...................... ............. ........33 32 R 2 and R3 vs. w avenum ber. ........................................................... 36 33 R4 through R6 vs. wavenumber. ............................... ............... 38 34 The velocity vectors in both phases at the onset of instability. z = 0 represents the interface position in the base state................................... ....... 40 35 The critical temperature difference versus the dimensionless wavenumber for waterwater vapor system where d = 3mm and d* =1.5mm...............................41 36 Evaporation number vs. dimensionless wavenumber for two waterwater vapor system s ............................................................43 41 Heat transfer across the bottom plate. ....................................... 59 42 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 43 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 44 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 51 The physical model .................................. .............. ........ 75 52 The dependence of the phenomenological coefficient on the wavenumber ............83 53 Com prison of m odels ................................................. ............... 85 54 Effect of shallow fluid layers ........................................... ............... 87 55 Critical temperature difference vs. dimensionless wavenumber.............................. 88 61 Different types of convection mechanisms ........................................89 62 Schematic of the experimental setup............................................. .. ......... 91 63 Experim mental setup ................................................................. .. ......... 93 64 zcomponent of velocity in a silicone oilair system at the onset of instability where the silicone oil and air depths are 5mm and 3mm, respectively.................97 65 zcomponent of velocity in a silicone oilair system at the onset of instability where the silicone oil and air depths are 5mm and 5mm, respectively.................97 66 Flow profiles in a silicone oilair system at the onset of convection where the silicone oil and air depths are 5mm and 3mm, respectively. ................................98 67 zcomponent of velocity in a silicone oilair system at the onset of instability where the silicone oil and air depths are 5mm and 9mm, respectively.................99 68 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 69 Pictures of experiments taken by IR camera at the onset of instability ................. 100 Bl 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 fluidfluid interfaces are of interest to researchers on account of their role in many industrial and natural processes. These instabilities are often accompanied by phasechange 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 nonevaporative systems, experiments were performed in a liquidgas 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 everincreasing 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 phasechange 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 dryeye 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 phasechange problem and later for connected instabilities such as convection and gravitational instability. Physics of Various Instabilities Connected to VaporLiquid 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 11, which depicts a liquidvapor 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 11. r COLD PLATE r LG Vapor LL Liquid JA HOT PLATE Figure 11. Physical model of the phasechange 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 12a and 12b depict an evaporation problem at a zero evaporation rate. Figure 12a shows the heating arrangement where the heat needed for evaporation is supplied from the liquid side while Figure 12b 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 12a. 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 12b. 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 13 is similar to Figure 12.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 phasechange 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 13. Again, instability occurs at every wavenumber. COLD PLATE ( N HOT PLATE Figure 13. Physics with fluid dynamics Now surface tension acts on pressure differences in fluid dynamics problems via the surface curvature and in many phasechange 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 RayleighTaylor 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 BuoyancyDriven Convection Buoyancydriven 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 14. COLD Figure 14. 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 topheavy 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 reestablished. 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 buoyancydriven 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 gradientdriven or Marangoni convection in the following section. Pure Marangoni Convection or Interfacial Tension GradientDriven Convection Interfacial tension gradientdriven convection or Marangoni convection is unlike buoyancydriven 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 15. HOTu Figure 15. 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 gradientinduced 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 buoyancydriven convection, there exists a critical temperature difference when the surface tension gradientdriven 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 15, 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, buoyancydriven convection is more prevalent, whereas for shallow depths, interfacial tension gradientdriven 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 gradientdriven 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 largescale 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). RayleighTaylor 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 topheavy arrangement called RayleighTaylor 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 RayleighTaylor 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 phasechange problem in Chapter 5. What makes the RayleighTaylor 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 RayleighTaylor 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 phasechange 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 noslip 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(vu) = 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(Vu)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+ +qS.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 phasechange 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 phasechange 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 TT 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 phasechange 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 phasechange and the thermodynamic equilibrium relationship is then ignored. It is also possible to consider phasechange with a zero temperature drop across the bilayer. However, in this case a nonzero evaporation base state is required. This case is not considered in this study. Second, when phasechange 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 phasechange. In other words, assume that the upper plate pressure is fixed so that the phasechange rate in the base state is zero. This of course does not mean that the phasechange 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 1j I Z P P which, in the absence of phasechange, becomes Vz1 = 7 Z, O and T, = 0 0. 2.28 vTi = C Z,. The continuity of temperature and the noslip 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 phasechange simply becomes k* dT1 dT, k dz dz The equilibrium condition at the interface, when there is phasechange 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 108 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 phasechange problems: Interface conditions There is a crucial difference between the classical RayleighMarangoniB6nard problem and the evaporation problem. The difference lies in the interfacial conditions. In the RayleighMarangoniBenard, 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 phasechange 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 BenardMarangoni 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 phasechange 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 nonzero 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 phasechange 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 phasechange problem is discussed. The results for different cases are presented in the form of plots and physical explanations are given. The Pure PhaseChange Problem To understand the physics of the pure phasechange problem, we first display the results of calculations using the full set of phasechange problem equations as mentioned in Chapter 2. The only terms left out of the equations are the terms that have gravity and surface tensiongradient 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 31. Physical properties of water and water vapor at 1000C p = 960 p = 0.6 m m = 2.9.104 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.103 1 C C S= 2.3.106 J 5.810 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 phasechange, Marangoni convection, gravity and buoyancydriven convection, respectively. Table 32. Definitions ofR1, R2 and R3 AT pATcMa R  TpcMa ATpcMaRa ATpcMa For example, ATpeua denotes the critical temperature difference for the case where there is phasechange and Marangoni convection, and R' therefore is the ratio of AT Mca 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 phasechange problem is compared to the critical temperature difference for a case where the Marangoni effect is added to the phasechange problem. This comparison is given in Table 33 and the readers can see for themselves that the addition of the Marangoni effect does little to change the threshold of evaporative instability. Table 33. Comparison of critical temperature differences Wavenumber, co ATpc (C) ATpcMa (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 phasechange 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 phasechange whereas at the crest the heat is released, and the otherwise low crest temperature goes up. In other words this phasechange action reduces the temperature perturbations along the interface in phasechange problems, which, in turn, reduces the surface tension gradientdriven convection. Marangoni convection therefore ultimately never plays a significant role on the stability of the phasechange 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 31. R' vs. wavenumber. In Figure 31, we observe that at large wavenumbers the pure Marangoni problem is considerably more unstable than the phasechange 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 31. In producing Figure 31 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 phasechange across the interface between the liquid and its vapor. In this case, the flow can only be caused by surfacetension gradients along the interface. In the other case, we allow phasechange 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 phasechange rate is zero. As a result, phasechange is only possible through perturbations at the interface. Accordingly, we have two systems with noflow in their base states. They are two similar systems with only one difference; one has the possibility of phasechange 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 phasechange 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 buoyancydriven 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 32. 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 phasechange, buoyancydriven convection can easily be stabilized at very small and very large wavenumbers. Consequently, in Figure 32 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 buoyancydriven convection, the stability offered by gravity is swept away by the stabilizing effect of the phasechange mechanism and the curve, depicted by R2, merges with the unity line. This occurs because the surface curvature and, in turn, the phasechange becomes more important compared to gravity at large wavenumbers. In contrast, R3 diverges to become less stable because of the destabilizing character of buoyancydriven convection. Yet, at even larger wavenumbers R3 ultimately merges with the unity line, not shown in Figure 32. R 0.5 1 2 3 4 5 6 7 8 9 Dimensionless wavenumber, mo Figure 32. R2 and R3 vs. wavenumber. In order to demonstrate the crucial role of the active vapor layer in the stability of a bilayer phasechange problem, we vary the depth of the vapor layer for a constant liquid depth and calculate the critical temperature difference for a phasechange problem with Marangoni and Rayleigh convection. The results are plotted in Figure 33, 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 34. Table 34. 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 33 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 phasechange 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 33, 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 phasechange 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 33. 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 phasechange problem in the absence of gravity and Marangoni effects. In Figure 34 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 downflow 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 phasechange 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 phasechange 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 phasechange 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 34. 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 cutoff 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 35. AT (C) 5 0 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5 6 Dimensionless wavenumber, o Figure 35. The critical temperature difference versus the dimensionless wavenumber for waterwater 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 phasechange 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 phasechange 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 twophase 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 buoyancydriven and surface tension gradientdriven 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 PhaseChange Problem The phasechange 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 phasechange problem exhibits a minimum as seen in Figure 36. 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 36. Evaporation number vs. dimensionless wavenumber for two waterwater 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 PhaseChange Problem In Chapter 3, linear stability methods were used to analyze the phasechange 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 phasechange 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 IRimaging, 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 liquidvapor 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 postonset 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 phasechange 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 phasechange 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 phasechange 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 (VP)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 2u ( 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 phasechange boundary is HKE (PPBASE) 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+(AAc) u,+z au (A Ac) u +2 2z U1+z 0 +z I zoj 2 1zo 2z0 zo +(AAc) 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 selfadjoint. 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 liquidvapor interface as the differential operator for Z at every order is selfadjoint. 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 transcritical; 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 phasechange 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 (EEc) VPV +V0,)+ (EEC2 V V PIPV2) 2 10 4.9 + (E E )3V P3 V3 V0v 2 2V0 1)=0 and (E Ec) _* V 1 + ( EEC 2 OP2 +V _ V 2. ) 2 V P v* Pr 1 4.10 +(EEC)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 (EEC ,2 2 V OTo 2v VT ) 4.11 + (EEC )3a T3 VOTo 3 ( OT +, VoT,))=0 and VTo +(EEc) VoTi  I VoTo \ K K +l(EEC 2 T2 2 OTo* 2 V VoT* 4.12 2 r r +(EEc) 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 (EE (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 +(EEc)a 8Z + l1 8EZ)( aZ azj,2 'I +(EEC) 2. \2 1 O 0 2 xo 0 8xo 0 4.15 1 C8Z + (EEC)3 6 Ox t +(EEc) + I+(EEC)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 (EEc) vz +(EEc ) v2 + (EEc)3 vz3 2 6 (EEc) vx + (EEc) vx2 + (EEc) vx3 2 6 4.16 4.17 4.18 TO +(EEcaT, 1 (E E,)2aT2 + (EE,)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(EEc)v 12 + (EEc v3 = 0 4.20 2 6 / \1 2x \3 1 (E Ec) v, +(EE)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* vIl 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 vz33v, 2 +3 zZ+3 c Z2 zZ 8 xa0 x zo a zo a 3v 2 3v2, v6 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 12xi + +2 a21Z, c8zo 0 80x 8zo xo Czq I 12vx, + 6 vz2 + 2Z 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) +(EEfc) + (EE )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 qn= 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 41 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 41. 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 42 and 43, 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 phasechange 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 103 Figure 42. 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 x103 Figure 43. 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 44, 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 phasechange 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 44. 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 nonzero, and it is A 3 f (f))2 4;r f"(0) So far, one of the three possibilities gives us a nonzero 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, vt=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 xindependent part of the particular solution for U2. The domain equations in the liquid for xindependent 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 RAYLEIGHTAYLOR PROBLEM WITH PHASECHANGE 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 RayleighTaylor 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) gdj2 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 RayleighTaylor 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 RayleighTaylor 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 AdhamKhodaparast 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 RayleighTaylor instability result. In this chapter we therefore investigate the RayleighTaylor 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 51). COL4XD PUATE Vapor HWPLATE M Figure 51. 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 phasechange rate at the interface can be controlled by adjusting the vapor pressure at the vaporside wall. The phasechange 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 phasechange 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 phasechange rates at the interface in the base state. The equations that model the physics are given by the NavierStokes, 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 noslip condition applies along the top and bottom plates. Hereon, vx, vz, and T are taken to be the x and zcomponents 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 noflow, noslip 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 noslip 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 RayleighTaylor 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 ewzo +B1 z, e"wzo v1 = A ewzo +Bg Zo eWZo T = C1 eWZO 2 A + zo ewzo +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 AlwB' =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 w2) 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(1p)+ (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 AdhamKhodaparast 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 AdhamKhodaparast 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 ww =0 When we assume that the kinematic viscosities of both fluids are equal to each other, we get the familiar result of AdhamKhodaparast et al. in dimensionless form. It is 2 av d 2 w +4 wwRT = 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 52, 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 31. Figure 52 shows the strong dependence of the phenomenological coefficient on the wavenumber in the plotted wavenumber range. The classical RayleighTaylor dimensionless critical wavenumber for the same system is 3.4e We observe that at exactly the classical RayleighTaylor critical wavenumber the value of the phenomenological coefficient is equal to zero. When we look at the relation derived by AdhamKhodaparast 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 52. 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 53 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 RayleighTaylor problem that the system is unstable for wavenumbers smaller than the critical wavenumber, however, when there is phase change taking place at the liquidvapor 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 53 we can tell that the assumption of infinitely deep layers is also not a good one for all wavenumbers. For a waterwater 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.1meter 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 53. 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 RayleighTaylor instability, the vapor layer is shallower than the liquid layer. For that reason, in Figure 54 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 phasechange 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 55 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 54. Effect of shallow fluid layers. a) d'= 50e3m and d = 50e3m b) d'= 50e3m and d9= 5e3m Note that heating from below can stabilize an arrangement where a vapor underlies its own liquid at a wavenumber at which the classical RayleighTaylor 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. 