UFDC Home  myUFDC Home  Help 
Title Page  
Acknowledgement  
Table of Contents  
List of Tables  
List of Figures  
List of abbreviations and...  
Abstract  
Introduction  
Literature review  
Modeling equations  
Linearized equations  
Spectral solution method  
Experimental design  
Results and discussion  
Conclusions and possible future...  
Appendices  
References  
Biographical sketch 



Table of Contents  
Title Page
Page 1 Page 2 Acknowledgement Page 3 Page 4 Table of Contents Page 5 Page 6 Page 7 List of Tables Page 8 List of Figures Page 9 Page 10 Page 11 Page 12 List of abbreviations and symbols Page 13 Page 14 Page 15 Abstract Page 16 Page 17 Introduction Page 18 Page 19 Page 20 Page 21 Page 22 Page 23 Page 24 Page 25 Page 26 Page 27 Page 28 Page 29 Page 30 Page 31 Page 32 Page 33 Page 34 Page 35 Page 36 Page 37 Page 38 Page 39 Page 40 Page 41 Page 42 Page 43 Page 44 Page 45 Page 46 Page 47 Page 48 Page 49 Literature review Page 50 Page 51 Page 52 Page 53 Page 54 Page 55 Page 56 Page 57 Page 58 Page 59 Page 60 Modeling equations Page 61 Page 62 Page 63 Page 64 Page 65 Page 66 Page 67 Page 68 Page 69 Linearized equations Page 70 Page 71 Page 72 Page 73 Page 74 Page 75 Page 76 Page 77 Page 78 Page 79 Spectral solution method Page 80 Page 81 Page 82 Page 83 Page 84 Page 85 Page 86 Page 87 Page 88 Page 89 Page 90 Page 91 Page 92 Page 93 Page 94 Page 95 Page 96 Page 97 Experimental design Page 98 Page 99 Page 100 Page 101 Page 102 Page 103 Page 104 Page 105 Page 106 Page 107 Page 108 Page 109 Page 110 Results and discussion Page 111 Page 112 Page 113 Page 114 Page 115 Page 116 Page 117 Page 118 Page 119 Page 120 Page 121 Page 122 Page 123 Page 124 Page 125 Page 126 Page 127 Page 128 Page 129 Page 130 Page 131 Page 132 Page 133 Page 134 Page 135 Page 136 Page 137 Page 138 Page 139 Page 140 Page 141 Page 142 Page 143 Page 144 Page 145 Page 146 Page 147 Page 148 Page 149 Page 150 Page 151 Page 152 Page 153 Page 154 Page 155 Page 156 Page 157 Page 158 Conclusions and possible future studies Page 159 Page 160 Page 161 Page 162 Appendices Page 163 Page 164 Page 165 Page 166 Page 167 Page 168 Page 169 Page 170 Page 171 Page 172 Page 173 Page 174 Page 175 Page 176 Page 177 Page 178 Page 179 Page 180 Page 181 Page 182 Page 183 Page 184 Page 185 Page 186 Page 187 Page 188 Page 189 Page 190 Page 191 References Page 192 Page 193 Biographical sketch Page 194 

Full Text  
CONVECTIVE INSTABILITY IN ANNULAR SYSTEMS By DARREN MCDUFF 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 2006 Copyright 2006 By Darren McDuff ACKNOWLEDGMENTS I begin by thanking my advisor, Dr. Ranga Narayanan, whose great ideas, support, and motivation have guided me through this work. I thank my mom, dad, and sister. They have always been the best of friends to me and been there for me during any difficult times in my life. On the subject of family, I should point out here that I first considered the idea of graduate school as a result of my dad suggesting it. I am also grateful to those on my advisory board. These include Dr. Martin Volz, Dr. James Klausner, Dr. Dmitry Kopelevich, and Dr. Mark Orazem. Their different viewpoints and ideas about my research have helped me to improve my work and avoid potential mistakes. Dr. Volz is my mentor in the N.A.S.A. Graduate Student Researchers Program, and traveled from Marshall Space Flight Center in Huntsville, Alabama, to be here for my presentations. Dr. Klausner comes from the Department of Mechanical and Aerospace Engineering at this university. Dr. Kopelevich and Dr. Orazem are from my own department. I must express my gratitude to Ken Reed. Ken is a talented machinist with T.M.R. Engineering. Highquality components made by T.M.R. were essential to the design of my experimental apparatus, and it was great to work with someone as helpful, knowledgeable, and professional as Ken throughout the project. Additionally, I thank Gerard Labrosse, a professor from L.I.M.S.I. in Orsay, France. During visits with our research group here at the University of Florida, he has spent a great deal of time sharing his knowledge on computational fluid dynamics (in particular, spectral methods) and providing helpful answers to questions from my research group. I would also like to thank a couple of undergraduates who have worked with me on this project, particularly in the design of my experimental apparatus and the process control program for it. Joe Cianfrone helped me to understand the computer software, Lab VIEWTM, and with that software, we designed a very effective process control program for my experiment. Joe also helped me to design a large electrical circuit to send data from my experiment into a process control computer, and send process control signals from that computer to components of my experiment. Scott Suozzo helped me with the design and redesign of several parts of my experimental apparatus, which helped the experiment to run well enough to produce results that match theoretical predictions. Of course, I thank my research group and the other graduate students in this department whose advice and ideas have helped me get past difficult parts of my project and helped me in my ongoing education in this field. In particular, I thank Dr. Weidong Guo, of my research group, for spending a good deal of time discussing computational issues with me during the last couple of years. Lastly, I am deeply grateful to N.A.S.A. for supporting this work through a fellowship from the N.A.S.A. Graduate Student Researchers Program (grant number NGT852941). TABLE OF CONTENTS A C K N O W L E D G M EN T S .................................................................................... 3 L IST O F T A B L E S .................................. ...................................................... 8 L IST O F F IG U R E S ......................................................................................... 9 LIST OF ABBREVIATIONS AND SYMBOLS.................................. ............... 13 A B S T R A C T ....................................................................................................... 1 6 CHAPTER 1 IN TR O D U C TIO N ................................................................... .......... 18 1.1 Physics of Convection ...................................... ................ ........... 18 1.1.1 Rayleigh Convection ................................. .............. .......... 19 1.1.2 M arangoni C onvection............................................ ... .............. 21 1.2 R ayleigh N um b er...................................................... ... .............. 24 1.3 Physical Explanation: Pattern Selection......................... 25 1.3.1 Laterally Unbounded System................................................ 25 1.3.2 Rectangular System, Laterally Bounded, Periodic Lateral Boundary Conditions ........... .. .. ........... .. ..... ...... .............. 28 1.3.3 Rectangular System, Laterally Bounded, NonPeriodic Lateral Boundary Conditions .............................................. .......... 31 1.3.4 Cylindrical System ................................................. ......... 33 1.3.5 A nnular System .................................................... .......... 36 1.4 Annulus vs. Cylinder ........................................ ............... ........... 37 1.5 Application ................................................................... .......... 40 2 LITERATURE REVIEW ........................................................... ........... 50 2.1 Single Fluid Layers ......................................................... .......... 50 2.2 Multiple Fluid Layers ........................................................ .......... 56 3 MODELING EQUATIONS................................................................. 61 3.1 N onlinear E quations.................. .......................................................... 61 3 .2 S calling ........................................................................ .. . .... 66 4 LINEARIZED EQUATIONS............................................................ 70 4.1 Linearization ..................................... ..................... .................... 70 4.2 Expansion into Normal Modes................................... ............ ...... 75 5 SPECTRAL SOLUTION METHOD.................................................. 80 5.1 Explanation of the M ethod ................................................. ......... 82 5.2 Application of the M ethod ................................................. ......... 89 6 EXPERIMENTAL DESIGN.............................................................. 98 6.1 G oals in Experim ental D esign............................................. ................ 98 6.2 Experimental Apparatus ...................................... ......................... 98 6.2.1 Test Section ..................................... ............................... 99 6.2.2 Bottom Temperature Bath............................. .................. 101 6.2.3 Top Temperature Bath.................................. 102 6.2 .4 Insulation ........................................................... . .... 103 6.2.5 Flow Visualization .................................. ............................ 103 6.3 Process C control System ................................................... .......... 104 6.4 Typical Experimental Procedure.................................................. 106 7 RESULTS AND DISCUSSION........................................................... 111 7.1 Cylindrical System s ............................ .. ...... .... ............................ 111 7.1.1 Constant Viscosity Computations: Cylinder.............................. 111 7.1.2 NonConstant Viscosity Computations: Cylinder....................... 113 7.1.3 Experiments: Cylinder ............................. .............. .......... 114 7.1.4 Comparison of Results: Cylinder....................................... 115 7.2 A nnular System s .................................................... .... ................... 116 7.2.1 Constant Viscosity Computations: Annulus.............................. 117 7.2.2 NonConstant Viscosity Computations: Annulus....................... 121 7.2.3 Experiments: Annulus......... ........... .............. 121 7.2.4 Comparison of Results: Annulus.......................................... 122 7.3 Error A analysis .............................................................. .......... 125 8 CONCLUSIONS AND POSSIBLE FUTURE STUDIES............................... 159 8 .1 S u m m ary ................................................................................... 1 5 9 8.2 Future Studies .............................................................. .......... 161 APPENDIX A THERMOPHYSICAL PROPERTIES......................................................... 163 B ADDITIONAL EXPLANATION OF EQUATIONS......................... 165 B. 1 Boussinesq Approximation ........................... ........... 165 B.2 Nonlinearities in the Governing Equations......................................... 166 B.3 Characteristic Velocity and Characteristic Time........... .............. 167 C DEVELOPMENT OF MATHEMATICAL MODEL FOR THE CONSTANT V ISC O SIT Y C A SE ............................................................................... 168 C 1 N onlinear E equations .......................................................... .............. 168 C.2 Scaling ................................................................................ .. 170 C .3 L inearization ......................................................... ........ ........... 172 C.4 Expansion into Normal M odes................................................... 174 D MATLAB PROGRAM GENERAL FLOWDIAGRAM ..................... 177 E EXAM PLE M ATLAB PROGRAM S......................................................... 178 E.1 Physical Properties ........................................................ ........... 178 E.2 D epths for Annular System ................................ .. .......... ......... 179 E.3 Main Program: Annular System with TemperatureDependent Viscosity....... 179 F LABVIEWTM PROGRAM GENERAL FLOWDIAGRAM.............................. 191 R E F E R E N C E S ............................................................................................. 19 2 BIOGRAPHICAL SKETCH ............................................................... ........... 194 7 LIST OF TABLES Table Page 51 Example of convergence of computed result with N, and Nz................ ................ 94 52 Comparison of rectangular, 2D, nostress results with rectangular 1D results......... 94 53 Comparison of calculated cylindrical results with results of Hardin et al................ 94 71 Constant viscosity computational results: critical temperature difference and onset flow pattern for sets of cylindrical dimensions considered in experiments.............. 129 72 Nonconstant viscosity computational results: critical temperature difference and onset flow pattern for sets of cylindrical dimensions considered in experiments....... 129 73 Experimental results: critical temperature difference and onset flow pattern for sets of cylindrical dimensions ......................................................... .......... 129 74 Summary of computational and experimental results for cylindrical systems.......... 130 75 Constant viscosity computational results: critical temperature difference and onset flow pattern for sets of annular dimensions considered in experiments................. 130 76 Nonconstant viscosity computational results: critical temperature difference and onset flow pattern for sets of annular dimensions considered in experiments.......... 130 77 Experimental results: critical temperature difference and onset flow pattern for sets of annular dimensions ............................................... .............. .......... 131 78 Summary of computational and experimental results for annular systems.............. 131 Ai Thermophysical properties ....................................................... ........... 164 LIST OF FIGURES Figure Page 11 Simplified system diagram: heated from below........................................... 43 12 Rayleigh convection ................................................................ .......... 43 13 M arangoni convection .............................................................. .......... 44 14 General diagram: critical Rayleigh number vs. disturbance wavelength................ 44 15 Crosssectional velocity profiles of flow patterns with different numbers of convective rolls ...................................................................... ......... 45 16 General stability diagram for buoyancydriven convection in a rectangular, laterally bounded system, with periodic boundary conditions at lateral walls, for varying aspect r a tio ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...... ... ... ... ... ... ... 4 5 17 Crosssectional velocity profiles of fourroll flow patterns with nonuniform roll size... 46 18 General stability diagram for buoyancydriven convection in a rectangular, laterally bounded system, with nonperiodic boundary conditions at lateral walls, for varying aspect ratio ............................................................. .............. ......... 46 19 Diagrams of m = 0, 1, 2, 3 flow patterns: "U": upward flow, "D": downward flow... 47 110 General stability diagram for buoyancydriven convection in a cylindrical container, for varying aspect ratio ............................................................. .......... 48 111 Convective flow pattern in an annulus (Stork & Muller 1974)........................... 48 112 Example: computed flow profile for annular system, crosssectional view............. 49 113 Diagram of crystal growth system ................... ... ................ ....... ..... ..... 49 31 System diagram: cylindrical and annular systems........................................ 69 51 Discretization nodes in a cylinder ................................................. .......... 95 52 Exam ple: grid point spacing .......................................................... .......... 95 53 Diagram of matrix/vector arrangement of discretized problem........................... 96 54 Examples of flow patterns and parenthetical notation for cylindrical systems........... 97 61 Simple experiment diagram ....................................................... ........... 108 62 Photograph of experimental apparatus....................................... ...... .... 108 63 Photograph of flowthrough top water bath of experimental apparatus................. 109 64 Crosssectional diagram of test section................................................ 109 65 Process control system layout ..................................................... .......... 110 71 Constant viscosity computational results: Racrt vs. A for cylindrical systems......... 132 72 Constant viscosity computational results: Racrt vs. A for cylindrical systems, closeup view ........................................................................ .......... 133 73 Computed velocity profile, crosssectional view, Lz = 6.85 mm, Lr = 8.75 mm, m = 0:(0,2), [ A T]crit = 5.80 C .................................................... .......... 134 74 Computed velocity profile, crosssectional View, Lz = 7.18 mm, Lr = 11.51 mm, m = 1:(1,3), [ AT]rit = 5.04 C .................................... ................ .......... 135 75 Computed velocity profile, crosssectional view, Lz = 6.53 mm, Lr = 11.28 mm, m = 1:(1,3), [ A T]crit = 6.43 C .................................................... .......... 136 76 Photo of onset flow pattern, Lz = 6.85 mm, Lr = 8.75 mm, m = 0:(0,2).................. 136 77 Photo of onset flow pattern, Lz = 7.18 mm, Lr= 11.51 mm, m = 1:(1,3)................ 137 78 Photo of onset flow pattern, Lz = 6.53 mm, Lr = 11.28 mm, m = 1:(1,3)................ 137 79 Constant viscosity computational results: Racrt vs. S for annular systems with A = .7 5 .............................................................................................. 1 3 8 710 Constant viscosity computational results: Racrt vs. S for annular systems with A = 1 .2 8 ............................................................................................. 1 3 9 711 Constant viscosity computational results: Racrt vs. S for annular systems with A = 1.28, closeup view ............................................................ .......... 140 712 Constant viscosity computational results: Racrt vs. S for annular systems with A = 1 .6 0 ............................................................................................. 1 4 1 713 Constant viscosity computational results: Racrt vs. S for annular systems with A = 1.60, closeup view ............................................ ................ .......... 142 714 Constant viscosity computational results: Racrct vs. S for annular systems with A = 1 .7 3 ............................................................................................. 1 4 3 715 Constant viscosity computational results: Racrt vs. S for annular systems with A = 1.73, closeup view ........................................................... ........... 144 716 Constant viscosity computational results: Racrct vs. S for annular systems with A = 2 .9 0 ............................................................................................. 1 4 5 717 Constant viscosity computational results: Racrct vs. S for annular systems with A = 2.90, closeup view ........................................................... .......... 146 718 Constant viscosity computational results: Racrct vs. S for annular systems with A = 3 .4 0 ............................................................................................. 1 4 7 719 Constant viscosity computational results: Racrct vs. S for annular systems with A = 3.40, closeup view 1 ......................................................... .......... 148 720 Constant viscosity computational results: Racrct vs. S for annular systems with A = 3.40, closeup view 2 ......................................................... .......... 149 721 Computed velocity profile, crosssectional view, Lz = 6.85 mm, Ro = 8.75 mm, S= .16, m = 0:(0,1), [AT]crit= 7.08 C ......................................................... 150 722 Computed velocity profile, crosssectional view, Lz = 6.85 mm, Ro = 8.75 mm, S= .30, m = 3:(3,0), [AT]crit= 7.67 C ......................................................... 151 723 Computed velocity profile, crosssectional view, Lz = 7.18 mm, Ro = 11.51 mm, S= .12, m = 0:(0,1), [AT]crit= 5.57 C ......................... ................................ 152 724 Computed velocity profile, crosssectional view, Lz = 7.18 mm, Ro = 11.51 mm, S= .50, m = 4:(4,0), [AT]crit = 7.09 C..................................... 153 725 Computed velocity profile, crosssectional view, Lz = 6.53 mm, Ro = 11.28 mm, S= .10, m = 2:(2,1), [AT]crit= 7.21 C..................................... 154 726 Computed velocity profile, crosssectional view, Lz = 6.53 mm, Ro = 11.28 mm, S = .40, m = 4:(4,0), [ A T]crit = 8.07 C..................................... 155 727 Photo of onset flow pattern, Lz= 6.85 mm, Ro= 8.75 mm, S= .16, m = 0:(0,1) ....... 155 728 Photo of onset flow pattern, Lz = 6.85 mm, Ro= 8.75 mm, S= .30, m = 3:(3,0) ....... 156 729 Photo of onset flow pattern, Lz = 7.18 mm, Ro = 11.51 mm, S= .12, m = 2:(2,1)........ 156 730 Photo of onset flow pattern, Lz = 7.18 mm, Ro = 11.51 mm, S= .50, m = 4:(4,0)........ 157 731 Photo of onset flow pattern, Lz= 6.53 mm, Ro = 11.28 mm, S= .10, m = 2:(2,1)........ 157 732 Photo of onset flow pattern, Lz= 6.53 mm, Ro = 11.28 mm, S= .40, m = 4:(4,0)........ 158 D1 MATLAB program general flowdiagram..................................... .......... 177 Fi LabVIEW TM program general flowdiagram.................................... .......... 191 LIST OF ABBREVIATIONS AND SYMBOLS a volumetric thermal expansion coefficient (C 1) 93 unit vector in zdirection E amplitude of small perturbation ic thermal diffusivity (m2/s) A eigenvalue of the eigenvalue form of the problem shown in Chapter 5 p/ dynamic viscosity (kg/(m*s)) P/R reference dynamic viscosity (kg/(m*s)) v kinematic viscosity (m2/s) p density (kg/m3) PR reference density (kg/m3) U dimensionless inverse time constant r transpose symbol for a matrix form of the xdirection dependence in a laterally bounded, rectangular system with periodic lateral boundary conditions 0 dimensionless group which appears in the motion equation when considering the dependence of viscosity on temperature A aspect ratio for cylindrical and annular systems As aspect ratio for rectangular systems A lefthandside matrix of the eigenvalue form of the problem shown in Chapter 5 B a constant related to the dependence of viscosity on temperature (C 1) B righthandside matrix of the eigenvalue form of the problem shown in Chapter 5 CP constantpressure heat capacity (J/(kg*OC)) Cv constantvolume heat capacity (J/(kg*OC)) D differentiation matrix f function representing the temperature dependence of viscosity F function representing the temperature dependence of viscosity, in terms of dimensionless temperature g magnitude of gravity (m/s2) g gravity, in vector form (m/s2) k thermal conductivity (W/(m*OC)) L.I.M.S.I. Laboratoire d'Informatique pour la Mecanique et les Sciences de l'Ingenieur Lr radius in a cylindrical system, or annular gap width in an annular system (m) Lx horizontal depth of fluid layer (m) Lz vertical depth of fluid layer (m) m azimuthal wave number Npr Prandtl number Nr parameter setting the number of discretization nodes in the rdirection Nx parameter setting the number of discretization nodes in the xdirection Nz parameter setting the number of discretization nodes in the zdirection N.A.S.A. National Aeronautics and Space Administration p modified pressure (Pa) p characteristic pressure (Pa) P pressure (Pa) Ra Rayleigh number Racrit critical Rayleigh number Ro outer radius of annular system R, inner radius of annular system S radius ratio for annular system S stress tensor t time (s) t characteristic time (s) tel time elapsed (s) (Appendix F) tseg,T time segment that each Tb setpoint remains in effect for (s) (Appendix F) tseg, VCR time interval between repetitions of the video cassette recorder's recording cycle (s) (Appendix F) T temperature (C) Tb temperature at bottom wall of system (C) Tb,sp setpoint value for Tb (C) (Appendix F) Tb,spj forj = (1, 2, 3, ...), series of setpoint values input for Tb (C) (Appendix F) TR reference temperature (oC) Tt temperature at top wall of system (oC) A T vertical temperature difference in fluid layer (oC) [ A T]crit critical vertical temperature difference for convection (oC) [ A T]guess guessvalue for critical vertical temperature difference for convection (oC) vj component of velocity in the jdirection Velocity (m/s) V characteristic velocity (m/s) X eigenvector of the eigenvalue form of the problem shown in Chapter 5 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 CONVECTIVE INSTABILITY IN ANNULAR SYSTEMS By Darren McDuff December 2006 Chair: Ranga Narayanan Major Department: Chemical Engineering Since natural convection can occur in such a wide range of systems and circumstances, much remains to be learned about the convective behaviors that appear in some systems. For example, convection in fluid layers of shapes other than cylinders or rectangles can be quite interesting, and it is this which is investigated in this research. This study focuses on Rayleigh convection, which is natural convection caused by buoyancy forces. The convective behavior of a system (referring to a fluid layer as the "system") depends upon the shape and dimensions of the system, the vertical temperature difference across the system, the thermophysical properties of the fluid comprising the system, and the characteristics of the disturbance given to the system. In particular, this research is concerned with the effects of the size and shape of the system. One important application of this study is in semiconductor crystal growth processes, such as the vertical Bridgman growth method. This growth method is commonly used, for example, in growing the semiconductor crystal, leadtintelluride. Crystal growth is usually thought of as a process in a cylindrical container. It is of interest to learn how the convective behavior of the system would differ if annular fluid layers were employed. The onset conditions for buoyancydriven convection (critical vertical temperature difference and flow pattern) in bounded, vertical, cylindrical systems and vertical, annular systems were determined both experimentally and theoretically. In the experiments, onset conditions were determined by flow visualization with tracer particles. The theoretical analysis involves bifurcation theory and computations using spectral methods. Computations show that the physics of the annular system are similar to those of the cylindrical system. More rolls of convecting fluid are included as the container is widened, and rolls may exist with either an azimuthal or radial alignment. The relations between the annular container dimensions and the onset conditions for convection are determined and shown. The agreement between computations and experiments is generally good. Some discrepancy in this agreement arises in cases in which the annular gap width is small, and a likely reason for this is explained. CHAPTER 1 INTRODUCTION This chapter provides an introduction to the physics of the convective phenomena investigated in this research. It then proceeds to a brief explanation of the dimensionless Rayleigh number which is very meaningful in regard to buoyancydriven convection. After this, an explanation is given regarding the competing physical phenomena which determine the number of convective cells that appear at the onset of convection and their sizes. Next, an explanation of why it is interesting to study the behavior of an annular system as compared to a cylindrical one is given. Lastly, an example of an industrial application to which this study has relevance is discussed. 1.1 Physics of Convection Many types of convective flows are possible, and many of them have been well studied. Convection can be classified as either forced or natural. In forced convection, flow is caused by external means. Flow of a liquid being pumped across a flat plate, for example, is forced convection. The flow that will be considered in this study, however, is of the other variety. Natural convection is convection in which flows are caused by forces that interact with thermophysical property variations. For example, natural convection in a system can be caused by buoyancy, interfacial tension gradients, or a combination of the two. When buoyancy is the cause of flow, this is called Rayleigh convection; if interfacial tension gradients cause flow, this is called Marangoni convection. If a system (here, "system" refers to a fluid layer) is going to convect, then the temperature difference corresponding to the onset of convection is determined by the direction from which the system is heated, the shape and dimensions of the system, the thermophysical properties of the fluid, and the characteristics of the disturbance given to the system. Regarding these factors, it is noteworthy that the fluid's thermophysical properties themselves may change significantly across the fluid layer depending on the magnitude of the vertical temperature difference to which fluid layer is subjected. This research focuses on the effects of the shape and dimensions of the system on convective behavior. Unless otherwise noted, the heating of the system referred to is uniform on the top and bottom surfaces. Also, for the physical explanations that follow, until otherwise noted, the systems discussed may be assumed unbounded in the lateral direction. 1.1.1 Rayleigh Convection A system capable of convection, in general, could be heated from below or above and include any number of liquid and vapor phases. For a physical explanation of Rayleigh convection, though, all that need be considered is a single layer of fluid, unbounded in the horizontal direction, heated from below, with the top and bottom walls held at constant temperatures. Understanding the origin of convection in this system will be sufficient to predict what is likely to happen in more complex systems. More importantly, a single fluid layer, heated from below, is precisely the system examined in this research. Generally, the density of a fluid decreases as temperature is increased. This means that the density at the top of the fluid would be higher than that at the bottom of the phase. Consider a zaxis along the vertical direction in the fluid layer, with the top surface at z = Lz, and the bottom surface at z = 0. This simple system is shown as Figure 11. This topheavy fluid layer is potentially unstable to gravity (shown in many figures as g) if subjected to a perturbation. If a particular disturbance to the system causes an element of fluid near the top of the layer to be displaced downward, then it will have a tendency to continue moving downward toward its gravitationally appropriate resting place. As this element of fluid moves downward, it displaces the lighter fluid surrounding it upward due to mass conservation. The downwardmoving fluid elements are heated once they reach the bottom of the system, and the upwardmoving fluid elements are cooled once they reach the top of the system, and so the process repeats. This process results in rolls of circulating fluid, or convective "cells," that form cellular flow patterns. Since the lowest density in the phase is at the bottom, a fluid element displaced downward would continue to move downward to z = 0 if its motion were not resisted in any way. Flow is resisted, however, by the dynamic viscosity of the fluid. Additionally, if the kinematic viscosity and thermal diffusivity of the fluid are sufficiently high, then the mechanical and thermal perturbations may quickly die out. It is when the vertical temperature difference across the system is increased beyond a certain critical value that this convectionRayleigh convectionwill occur; at that point, gravitational instability overcomes the damping effects of kinematic viscosity and thermal diffusivity. Rayleigh convection is stronger in a phase with lower kinematic viscosity and lower thermal diffusivity; this is because mechanical and thermal disturbances are damped out less quickly in a phase having lower values for those properties. Now, it might appear that Rayleigh convection results simply from large enough vertical temperature gradients. In truth, however, even with a sufficiently large vertical temperature gradient present, Rayleigh convection can begin only if disturbances with transverse (horizontal) variation are imposed on the system. A fluid element displaced by a disturbance must "feel" a relative difference in the density of the fluid horizontally neighboring it, and this can occur only as a result of disturbances with transverse variation. Note that the disturbance which causes the onset of convection need not necessarily be a mechanical one. A thermal disturbance with transverse variation could just as easily cause convection to begin by making a certain fluid element warmer or cooler than the fluid horizontally adjacent to it. Figure 12 shows a buoyancydriven convective flow pattern occurring in a fluid phase. 1.1.2 Marangoni Convection Convection can occur in another way, as well. The flow can be driven by interfacial tension gradients rather than buoyancy forces. This type of flow is called Marangoni convection. Since Marangoni convection is driven by gradients in interfacial tension, it is a phenomenon that can occur only in systems possessing multiple fluid phases which contact each other at interfaces. This research is concerned with the simpler case of a single fluid layer confined in all directions by solid walls. Thus, Marangoni convection could not occur in the system being researched. Still, a brief explanation of Marangoni convection will now be given for the sake of giving a more complete explanation of the physics of convection. To this end, consider a system of two vertically stacked fluids, in which one vertical boundary is heated so as to subject the system to a vertical temperature gradient. Again, assume the system is unbounded in the horizontal direction. In this hypothetical system, both the top and bottom walls are kept at constant temperatures and Marangoni convection can occur in this system regardless of which vertical boundary of the system is warmer. For the sake of making an explanation by example, assume that the bottom wall of the system is warmer than the top. Like Rayleigh convection, Marangoni convection results from perturbations with transverse variation. If the interface between the fluids is perturbed, then some interfacial regions could be pushed upward and some downward. This gives rise to a transverse temperature gradient along the interface. Consider an interfacial region pushed upward. The fact that the system is heated from below means that this region now experiences a lower temperature than the regions of the interface adjacent to it. Typically, surface tension, like density, increases with decreasing temperature. Thus, the interfacial region that has been pushed upward has an increased interfacial tension and pulls fluid at the interface toward it. When this happens, fluid from the bulk phases must move in to replace fluid being pulled away from regions of the interface that lowered in interfacial tension (when the system is heated from below, these regions are the troughs of the wavy interface). The flow patterns in Marangoni convection, like those in Rayleigh convection, are cellular. Since the interfacial tension gradients are generated by temperature differences along the interface, a higher vertical temperature difference across the system favors Marangoni convection. Figure 13 is a simple illustration of Marangoni convection. The replacement fluid arriving from the bulk phases is warm in the case of the lower fluid, and cool in the case of the upper fluid; this supply of replacement fluids can either strengthen or weaken the temperature gradient which fuels Marangoni convection. The effect that the fluid moving to the hot regions of the interface from the bulk phases has on this temperature gradient is dependent upon the amount of fluid from each phase that arrives at those locations and on the thermophysical properties of each phase. As with Rayleigh convection, dynamic viscosity resists flow and can damp the process. Marangoni convection, like Rayleigh convection, is stronger in a phase with lower kinematic viscosity and lower thermal diffusivity because a phase having lower values for those properties damps out disturbances less quickly. On the other hand, if these two properties have high enough values, a disturbance at the interface will die out rather than amplify to cause convection. Of course, the phase with a lower kinematic viscosity is not necessarily the phase that has a lower thermal diffusivity. Thus, the combination of these two properties determines which phase will convect more strongly. The phase which has a kinematic viscositythermal diffusivity combination more favorable to convection is the phase from which it will be easiest for fluid to move in and replace that which is lost at hot interfacial regions during Marangoni convection; it is primarily from this phase that fluid is supplied to the hot regions of the interface. If this more strongly convecting phase is the warmer phase, then the fluid it supplies to the hot interfacial regions is even warmer than what was originally present, and Marangoni convection is reinforced. In the opposite scenario, cold fluid moves in to the hot regions of the interface; this cools those regions down and counteracts Marangoni convection. In the case where system of two stacked fluids is comprised of a vaporliquid bilayer, rather than a liquidliquid bilayer, it should be noted that a liquid can generally transport more heat by convection than can a vapor. This is quantitatively shown when a comparison is made between the product of density and heat capacity (p Cp) for liquids and vapors. Though the thermal conductivities of vapors and liquids may differ, this difference may be considered less significant. The reason for this is that the fluids are already conducting heat in the motionless initial state of the system (base state). While the effect of the difference in density multiplied by heat capacity becomes important once convection begins, the effect of the difference in thermal conductivity is relatively unchanged. By this line of reasoning, replacement fluid arriving to troughs on the interface from the liquid phase of a vaporliquid bilayer system is likely to have more of an effect on the temperature gradient along the interface than replacement fluid arriving from the vapor phase of the bilayer. The same reasoning (involving the product of density and heat capacity) also can be used in deciding which fluid will more strongly affect the interfacial temperature gradient in a liquid liquid bilayer system. It is also very interesting to note that Marangoni convection can occur in the system even if the interface remains flat. If a nonmechanical, thermal perturbation is given to the interface, then some regions become hot and some cold. Fluid from the hot regions is pulled toward the cold regions by interfacial tension, and fluid at the cold regions is displaced. As this happens, fluid from the bulk must move in to replace that lost at the hot interfacial regions this marks the start of convection. 1.2 Rayleigh Number The focus of this explanation will get back on its main course and address Rayleigh convection once again. Scaling of the equations that model a fluid layer heated from below leads to the following definition of the dimensionless Rayleigh number, which is shown as Equation 1.1. The Rayleigh (Ra) number is a quantity that relates the factors determining whether or not buoyancydriven convection occurs. a*g*L3 *AT Ra = (11) v*K The length scale associated with Rayleigh convection is the vertical phase depth, shown as Lz. In this definition, a is a volumetric thermal expansion coefficient, g is the magnitude of gravity, AT is the vertical temperature difference across the layer, v is the kinematic viscosity of the fluid, and K is the thermal diffusivity of the fluid (which is equal to the thermal conductivity divided by the product of the density and heat capacity). A higher Rayleigh number corresponds to a system that is more susceptible to buoyancy driven convection. Notice that the factors which favor buoyancydriven convection are all grouped in the numerator, while those which oppose the flow form the denominator. The thermal expansion coefficient represents the degree to which the density of the fluid changes when it is heated, so of course a higher value for this works in the favor of convection. Gravity and the vertical temperature difference are the most obvious of the two components needed to drive Rayleigh convection. A larger vertical temperature difference makes the system more unstable to Rayleigh convection. Having a larger vertical phase depth means that the fluid feels a weaker effect from the noslip conditions at the top and bottom walls, and thus is less mechanically constricted, and can more easily flow in a convective pattern. This effect is quite strong, and accordingly, the Rayleigh number is proportional to the vertical phase depth raised to the third power. Kinematic viscosity and thermal diffusivity act to dissipate any disturbances given to the system and prevent convection. Thus, these two thermophysical properties are located together in the denominator of the Rayleigh number. 1.3 Physical Explanation: Pattern Selection Some of the most interesting facets of the convective behavior studied in this research relate to the competition of certain physical phenomena to determine precisely which convective flow patterns arise at the onset of convection. 1.3.1 Laterally Unbounded System Imagine, once again, the system shown in Figure 11. The system is constantly subjected to small mechanical and thermal disturbances. As mentioned earlier, only disturbances with transverse variation can cause the onset of convection. The disturbances that the system receives are of an extremely large range of wavelengths, especially considering the case in which the system is laterally unbounded. Still, only one of these many disturbances actually causes the onset of convection. It is this disturbance to which the system is most unstable. It is the wavelength of this disturbance which determines the wavelength of the onset convective flow pattern. It is interesting to ask what makes disturbance of one wavelength more stabilizing or destabilizing to the system than another. For the moment, consider two possibilities in terms of the wavelength of a periodic, waveshaped disturbance. One possibility is that the disturbance is of very short wavelength, meaning that it is quite jagged and choppy in appearance. The second possibility is that the disturbance is of very long wavelength, and so it is very gently, gradually sloping, with its minimumpoints and maximumpoints spread far apart from one another. The disturbance of very short wavelength possesses a great deal of transverse variation, since it repeats its lateral oscillations so frequently; this high degree of transverse variation favors convection. However, since the adjacent peaks and troughs of the waveshaped disturbance are so close to one another, they easily diffuse their momentum and thermal energy into the horizontally neighboring fluid, which dissipates the disturbance and stabilizes the system to convection. For example, a fluid element which initially was heated and lowered in density relative to the fluid laterally adjacent to it, as a result of a disturbance of very short wavelength, would quickly match the density of the laterally adjacent fluid due to rapid diffusion of heat from the fluid element. Thus, insofar as diffusion of heat and momentum is concerned, the system is quite stable to disturbances which are of very short wavelength. The system reacts to a disturbance of very long wavelength a bit differently. A disturbance of very long wavelength is so laterally spread out that the sort of rapid diffusive stabilization just described does not play a role. Consequently, the lateral density differences generated by a disturbance of very long wavelength remain present rather than dissipating into uniformity, which is a condition that favors buoyancydriven convection. However, when the disturbance is of very long wavelength, the adjacent peaks and troughs of the waveshaped disturbance are so far apart from one another that the local density differences between adjacent fluid elements (the transverse variations which are needed to drive buoyancy driven convection) are too weak to cause and sustain buoyancydriven convection. If a disturbance of very long wavelength were able to form a corresponding cellular, convective flow pattern of a very long wavelength, the flow pattern would not be able to be sustained; this is because it would simply require too much energy to drive the fluid in the convective cells across the very long, nearly horizontal paths which would need to be traveled within each convective flow cell. So, in summary, shortwavelength disturbances are not very destabilizing because they rapidly degrade by diffusion, and longwavelength disturbances are not very destabilizing because they are spread out so widely in the lateral direction that the variations they induce in the thermophysical properties of the fluid are not strong enough to drive flow. Thus, of the wide range of disturbances to which the system is subjected, the one disturbance which is most destabilizing to the system and causes the onset of convection will possess some intermediate wavelength, which is determined by the competing physical behaviors associated with short wavelength and longwavelength disturbances. The vertical phase depth, of course, plays a large role in governing which disturbances are most unstable, as well. At this point, one could imagine the appearance of a simple graph relating the stability of the system, in terms of the critical Rayleigh number (Rayleigh number corresponding to the critical vertical temperature difference at which convection begins), to disturbances of different wavelengths. The graph (Figure 14) can be drawn using a general, dimensionless wavelength and the dimensionless critical Rayleigh number, which is Raccrt; numerical values for either quantity are not necessary at this point. Once the system is at or above the critical vertical temperature difference and is subjected to the right disturbance, buoyancydriven convection begins. Still considering a onedimensional system, unbounded in the lateral direction, the flow pattern at the onset of convection possesses the wavelength of the disturbance which caused convection. 1.3.2 Rectangular System, Laterally Bounded, Periodic Lateral Boundary Conditions Now, consider a twodimensional system. The twodimensional system to be considered is the system that would be obtained if two lateral walls were simply added to the left and right sides of the system of Figure 11, with a nostress condition on velocity and an insulating condition on temperature imposed on each lateral wall. The horizontal width of this new system is called Lx and the vertical depth is called L This system can be examined to understand the interesting physics governing convective pattern selection. Due to the stressfree conditions on velocity and the insulating conditions on temperature at the lateral walls, the xdirection dependence of the velocities and temperatures in the system, which shall be called Vi, could be expressed in terms of cosines, in the form n Kx SCOS n = 0, 1,2, 3, .... (1.2) The number of convective cells which may form in the system, then, is related to Lx through this form of xdirection dependence, as is the size of the cells which form. Notice that n is an integer, so not just any pattern may arise in a laterally bounded system. Only those patterns which can be physically accommodated by the lateral width of the system may form, and so only disturbances with those shapes are physically admissible to the system. Of the set of disturbances that are physically admissible to the system, only one will be the most destabilizing and, when the critical vertical temperature difference has been exceeded, cause the onset of convection. In a laterally bounded system, for the same reasons as in a system without side walls, disturbances that are laterally very short or very long are not very destabilizing. Thus, again the shape of the most destabilizing disturbance and onset flow pattern is determined by competition between physical behaviors like those that were explained for shortwavelength and long wavelength disturbances. The size and shape of the onset flow pattern, of course, are heavily dependent on the dimensions of the system. In a system that is laterally very wide, the onset flow pattern would likely not be, for example, one very wide flow cell. A flow pattern such as that would not be sustainable for the physical reasons like those explained for disturbances of very long wavelength, which were already explained. The size of the system is typically described using the dimensionless aspect ratio. For a rectangular system, the aspect ratio is the ratio of the width of the system to the vertical height of the system, and it may be called Az. For a cylindrical system, the aspect ratio is the ratio of the radius of the system to the vertical height of the system, and it may simply be called A. For an annular system, the aspect ratio may also be called A, and it is the ratio of the outer annular radius (Ro) to the vertical height of the system. As the aspect ratio is varied, obviously the stability of the system to any given disturbance is varied with it. For example, a waveshaped disturbance which induces a convective flow pattern consisting of two rolls of convecting fluid would be most destabilizing at some particular value of A but if the aspect ratio were slightly increased (making the container laterally wider or vertically shorter) or decreased, then that same disturbance would not be as well accommodated by the lateral width of the system, and consequently would not be as well suited to destabilize the system. The reason for this is related to physical behaviors like those described for disturbances of very short and very long wavelengths, as described for the laterally unbounded system. What is meant by a convective pattern consisting of two rolls of convecting fluid is clarified in the following crosssectional velocity profiles of a convecting, laterally bounded system (Figure 15). Continuing with the example, continuously increasing the aspect ratio would make the twoconvectiverollinducing disturbance less and less destabilizing to the system, and thus make the corresponding twoconvectiveroll onset flow pattern more and more difficult for the system to form and maintain; as an onset flow pattern becomes more and more difficult to form and maintain, the critical vertical temperature difference corresponding to it increases, as does the corresponding critical Rayleigh number. At some sufficiently large aspect ratio, the system would find a disturbance that induces more than two convective rolls to be more destabilizing, and the corresponding onset flow pattern to be more energetically favorable. Supposing that increasing the Axz value caused the example system to become more unstable to a disturbance with that induced four convective rolls rather than two, the corresponding onset flow pattern would consist of four rolls of convecting fluid. The stability of the system to the fourconvectiverollinducing disturbance, and the degree of difficulty that the system would have in maintaining its corresponding onset flow pattern, are again dependent on physical behaviors like the shortwavelength and long wavelength behaviors described for laterally unbounded systems. As the Axz value is further and further increased, this type of patternswitching continues. If the A, value of the example system were further increased, the system would eventually reach an Az value at which it would be more favorable for the system to form and maintain an onset flow pattern with more than four rolls of convecting fluid. Supposing that the next flow pattern in the example consists of six rolls of convecting fluid rather than four, then difficulty that the system has in forming and maintaining this sixrollpattern would once again be determined by physical behaviors like the shortwavelength and longwavelength behaviors that have been explained earlier. Resulting from these physical behaviors, which correspond to the size and shape of an imposed disturbance, the exact critical temperature difference, and critical Rayleigh number, corresponding to each different flow pattern is dependent on the Ax value and has a minimum at some certain Axz value that is optimal in terms of how well a system with that aspect ratio physically accommodates the flow pattern. When a system is laterally unbounded, or a system has lateral walls that are stressfree and insulated (which means the system behaves essentially like a laterally unbounded system), then the minimum possible critical Rayleigh number for any possible flow pattern is the same. Based on all of this, a graph showing the relation between the critical Rayleigh number and the value of Az is given below as Figure 16. Note that in a laterally unbounded system, or in a laterally bounded system with periodic lateral boundary conditions, the convective patterns that can form are comprised of convective rolls which are all of equal lateral width. As an example of what is meant by this, see the flow pattern diagrams in Figure 15. In this figure, the roman numerals represent flow patterns with different numbers of convective rolls. Flow pattern "III" includes more convective rolls than "II" and flow pattern "II" includes more convective rolls than "I" does. Basically, aspect ratios at which the slope of the curve changes from positive to negative are aspect ratios at which there is a transition to a higher number of convective rolls being present in the onset flow pattern. 1.3.3 Rectangular System, Laterally Bounded, NonPeriodic Lateral Boundary Conditions When running experiments or considering practical applications, one does not encounter systems with periodic boundary conditions at their lateral walls. The more realistic case of a laterally bounded, rectangular system with nonperiodic lateral boundary conditions will now be detailed. The behavior in this system is essentially the same as what was just explained for the rectangular system with periodic lateral boundary conditions, but there are a couple of key differences: the convective rolls which form at the onset of convection are not necessarily of identical size to one another, and the minimum critical Rayleigh number for each possible flow pattern (considering different values of Axz) decreases as Axz increases rather than maintaining a constant value. Since the convective rolls which form at onset are not necessarily of identical size to one another, the flow pattern with four convective rolls in Figure 15 could instead appear as shown in Figure 17. The second key difference of this case from the case in which the lateral boundary conditions are periodic is that for higher and higher A, values, the minimum critical temperature difference values become lower and lower. The reason for this is that the noslip effects of the side walls do less to stabilize the system as the aspect ratio is increased. As A, approaches infinity, the critical temperature difference values for the flow patterns approach values corresponding to the well known critical Rayleigh number of 1708 that is associated with the onset of buoyancydriven convection in laterally unbounded systems (Davis 1967). For the two dimensional, laterally bounded system, noslip conditions at the lateral walls, here is a general plot of the stability of the system, in terms of the critical Rayleigh number, to disturbances corresponding to different onset flow patterns; the plot is simply based on the reasoning explained above. Again, in this figure, the roman numerals represent flow patterns with different numbers of convective rolls, where pattern "III" includes more convective rolls than "II" and pattern "II" includes more convective rolls than "I." Recall that aspect ratios at which the slope of the curve changes from positive to negative indicate transitions to onset flow patterns with higher numbers of convective rolls being present. 1.3.4 Cylindrical System So far, the example systems discussed have been assumed to be laterally bounded, rectangular systems or laterally unbounded systems. Radial walls in cylindrical and annular systems significantly affect the critical conditions for the onset of convection. This research addresses primarily cylindrical and annular systems. To prevent any possible confusion, note that throughout this dissertation, "cylindrical system" (or sometimes the term "open cylindrical system" may be used) refers to a fluid system bounded by a cylindrical container, in which the entire container is filled only with fluid, and does not include a centerpiece as an annular container does. From this point onward, unless otherwise noted, the systems discussed will be within vertical, cylindrical containers (meaning that the direction of gravity is parallel to the axis cutting through the radial center of the cylinder) with finite outer radii or within vertical, annular containers. To be completely clear, a vertical, annular container is one which has an annular crosssection when viewed from above. Increasing the aspect ratio, A, in a cylindrical system has the effects on pattern selection which were just explained for the laterally bounded, rectangular system. As A is increased, the onset flow patterns include more and more convective rolls, so as to occupy the system in a more energetically favorable fashion. Also, as A approaches infinity, the value of the critical Rayleigh number approaches the value of 1708, which is wellknown to be associated with laterally unbounded systems. There is, however, another aspect to pattern selection in cylindrical systems, and it will be explained now. In a rectangular system, the number of convective rolls which form can be related to the imposed disturbance's size and shape in the xdirection. For rectangular, laterally unbounded systems, the corresponding onset flow pattern can then be described by a wave number in the x direction, inversely proportional to the wavelength of the flow pattern in the xdirection. In all systems considered in this research, whether rectangular, cylindrical, or annular, it is most energetically favorable for the system to form onset flow patterns with only one convective roll in the vertical direction; thus, there is no need to consider a zdirection wave number, zdirection wavelength, or the size and shape of a disturbance in the zdirection for this research. The flow pattern in a rectangular system, then, can be described simply by describing its size and shape in the xdirection. In a cylindrical system, flow patterns may have different numbers of convective rolls in the radial direction, like the patterns in rectangular systems could have different numbers of convective rolls in the xdirection. However, flow patterns in cylindrical systems may also have differing numbers of convective rolls in the azimuthal direction. The flow patterns in a cylindrical system must, of course, be periodic. The periodicity of the flow patterns in the azimuthal direction, and thus how many convective rolls exist along the azimuthal direction, is characterized by an azimuthal wave number, called m. For example, the flow pattern described by the wave number m = 2 is one which repeats twice as one progresses once through the fluid layer in the azimuthal direction. For an m = 2 flow pattern, progressing 900 around the system in the azimuthal direction, from an arbitrary starting point, leads to a point that has the opposite vertical velocity of the starting point. Likewise, for an m = 1 flow pattern, progressing 1800 around the system in the azimuthal direction, from an arbitrary starting point, leads to a point that has the opposite vertical velocity of the starting point. An m = 0 flow pattern, though, is axisymmetric in the azimuthal direction. Figure 19 gives examples of flow patterns possessing some of the lowest m values. In these diagrams, the letter "U" indicates upward velocities and the letter "D" indicates downward velocities. The intensity of the shading indicates the magnitude of the velocities, with darker shades representing higher velocities. At certain aspect ratios, the a cylindrical system favors certain azimuthal periodicities, and thus, certain m values. For each m value, the system may have a differing number of convective rolls in the rdirection, and how many there are depends once again upon physical behaviors like the shortwavelength and longwavelength behaviors that have already been explained. The relation, for each m value, between the number of convective rolls, the critical Rayleigh number needed to obtain a particular number of rolls, and the aspect ratio appears much like what is shown in Figure 18. Generally, in the cylindrical cases, only m = 0, 1, or 2 patterns were favored as onset patterns; this means that in most cases, the critical vertical temperature difference corresponding to the onset of convection was lower for m = 0, 1, or 2 patterns than for patterns with other azimuthal wave numbers. As the aspect ratio approaches infinity, the critical Rayleigh numbers for flow patterns of all m values approach the unboundedsystemvalue of 1708. A general graph showing the relation of flow patterns and their critical Rayleigh numbers to the aspect ratio in a cylindrical system is given as Figure 110. This graph was obtained by performing a series of computations, to determine the critical Rayleigh number corresponding to several values of the aspect ratio, using the computational program developed for obtaining computational results in this research. The computational solution method is detailed in Chapter 5, and the computational results are shown extensively in Chapter 7. On this diagram, starting from a given aspect ratio and moving upward along the graph is analogous to increasing the vertical temperature difference across the system. Once the first of the m = 0, 1, 2, or 3 curves has been reached, then this is the critical Rayleigh number corresponding to the critical vertical temperature difference. The m value of the bottommost curve at a given aspect ratio is the m value for the onset flow pattern at that aspect ratio; this azimuthal wave number, since the critical temperature difference and critical Rayleigh number corresponding to it are lowest, is the azimuthal wave number of disturbances to which the system is most unstable. Notice that the curve for each particular m value looks much like the one in Figure 18, which reflects the effects of physics related to the sizes and shapes of different disturbances on the stability of the system. Graphs like Figure 110 are extremely helpful in research on buoyancydriven convection. The reader is asked to notice that the graphs presenting computational results for buoyancydriven instability in cylindrical systems in Chapter 7 are just like the one shown in Figure 110. Figure 110, interestingly, is the correct graph for cylindrical systems of any dimensions, and it holds true regardless of the set of thermophysical properties possessed by the fluid in the system. This will be further explained at the end of Chapter 4. Note that the stability diagrams in this research would, however, differ very slightly depending on whether or not the variation of viscosity with temperature is being considered. 1.3.5 Annular System Nearly everything that was just said for cylindrical systems in Section 1.3.4 applies for annular systems, as well. An exception is that as the outer radius of the system approaches infinity (and thus, the A value approaches infinity), one would expect the critical Rayleigh number to approach a value other than 1708, which corresponds to a laterally unbounded system, since an inner radial wall is present. Patterns in annular systems are described primarily by their azimuthal wave number, m. The patterns in annular systems may also include more than one convective roll in the radial direction, but this tends to happen only when the ratio of the inner annular radius (R1) to the outer annular radius (Ro) is quite small. This ratio of radii is a new dimensionless parameter affecting pattern selection which was not present in the consideration of cylindrical systems. The radius ratio shall be called "S" and will be given by Equation 1.3. S = R, (1.3) R The effect of the radius ratio is the key difference between the annular system and the cylindrical system in terms of the formation and selection of flow patterns. Stability diagrams of the relations between Racrt, the aspect ratio, and the radius ratio (S) for different wave numbers, and valid for systems with any set of thermophysical properties, could be made for the annular system, as well; they would appear somewhat like the one for the cylindrical system, shown as Figure 110. Such diagrams are shown in Chapter 7. As annular systems are the main focus of this research, they will be explained much more thoroughly in the following pages. 1.4 Annulus vs. Cylinder Certainly, there has been a good deal of research already done on natural convection in annular systems (Stork & Muller 1972, and Littlefield & Desai 1990). However, much of that work is considerably different from the current research. For example, many of those studies involve horizontal annuli rather than the vertical annuli addressed in the current study. Secondly, of those past works which do address convection in vertical annuli, many consider the case in which heat is supplied to the system from the inner rod of the annulus rather than the case considered in this study, which is that heat is supplied to the system from its bottom wall. The case in which heat is supplied from the bottom wall is highly relevant to industrial processes, and this will soon be discussed further. The radial walls of all cylindrical and annular systems considered in this research are assumed to be insulating. This assumption shall be maintained throughout this paper, as the radial walls in the experiments conducted for this research were made to be insulating. Whether the radial wall is treated as conducting or insulating does have an effect on what the critical vertical temperature difference is for the onset of convection. Generally, the system becomes unstable at a lower vertical temperature difference if the radial walls are insulating. This makes sense because insulating radial walls do not allow a system to dissipate thermal disturbances as easily as conducting radial walls would, and making the dissipation of thermal disturbances more difficult makes the system more unstable to convection. As one might expect, the difference between the insulating and conducting cases is more pronounced in systems that span a smaller distance in the radial direction. In such systems, the radial walls are closer to the interior of the system and so they have a relatively stronger effect on the behavior of the system. Also, whether the radial walls are insulating or conducting, it is obvious that if the system is of smaller radial extent, then the system will be more stable to convection than a wider system of matching vertical depth; this is because of the relatively larger effect of the noslip condition on velocity imposed at the radial walls. The addition of an inner block of circular crosssection to the center of a cylindrical system transforms the cylindrical system into an annular system. The inner block of the annular system is insulating, just like the outer radial wall, and varying the radius of this circular inner block changes the gap width occupied by fluid within the annular system. Changes in this gap width cause significant changes in the critical vertical temperature difference and flow pattern at the onset of convection in the annular system. Primary goals of this research were to better understand the phenomenon of buoyancydriven convection in terms of how it differs in an annular system compared to a cylindrical system, and how the onset conditions for convection in the annular system differ as the annular gap width varies. As one may expect, the addition of a circular inner block to the cylindrical system does, because of the additional radial wall and corresponding noslip condition on velocity that it introduces, impart to the system a higher stability to convection. That is to say, a cylindrical system of a certain vertical depth, which requires a certain vertical temperature difference for the onset of Rayleigh convection, would require a higher vertical temperature difference to convect if some small, circular block were introduced to the center of the system, converting it to an annular system. The circular centerblock present in an annular system also affects the flow patterns which may appear at the onset of convection. Stork and Muller (1974) have shown experimentally that in a onefluid annular system, heated from below, convective cells can form such that they fill the annular gap in an azimuthally aligned arrangement that is rather reminiscent of the arrangement of spokes on a bicycle tire. An example of this type of pattern presented in the work of Stork and Miller is shown as Figure 111. In this photograph of one of their experiments, the fluid convecting is a silicone oil, and the flow is visible due to the presence of aluminum tracer powder in the oil. When the flow is horizontal, the aluminum particles align in such a way as to reflect more light, making those sections of the oil appear lighter in color; regions of the oil which are flowing vertically appear darker in color. As the radius of the centerblock in the annulus is increased, the system becomes more stable to buoyancydriven convection, meaning that a higher vertical temperature difference across the system is needed for convection to begin; this makes sense because as the radius of the centerblock is further and further increased, there is more and more wallsurface area (and more of its accompanying noslip effect) compared to the amount of fluid in the system. Also, as the radius of the centerblock in the annulus is changed, the number of azimuthally aligned cells changes. Generally, a greater number of azimuthally aligned cells form at the onset of convection as the radius of the centerblock approaches the outer wall of the annulus, and a smaller number of azimuthally aligned cells form at the onset of convection as the radius of the centerblock is decreased. This research has shown that once the radius of the centerblock decreases past a certain value, new rolls of convecting fluid at the onset form in radial alignment (meaning that they span the entirety of the azimuthal direction and that they are concentric to the centerblock of the annular container) in addition to the original azimuthal arrangement. An example of this, generated by the computations done for this research, is shown in Figure 112. In this flow profile, the velocity is scaled, negative velocities represent downward flow in the z direction (which extends into and out of the plane of the diagram), and positive velocities represent upward flow in the zdirection. This crosssectional velocity profile represents a cross sectional region approximately half of the distance along the vertical depth of the system that it is computed for. If the radius of the centerblock of the annular system is further and further decreased, the radial alignment of the convective cells becomes more prominent compared to the azimuthal alignment of the cells, and the flow patterns become progressively more like those seen at the onset of convection in open, cylindrical systems. This transition in onset flow patterns was of great interest in this research. 1.5 Application Aside from its applications to the regular appearances of Rayleigh convection in nature, this study finds application in the drying of paint films, in smallscale fluidics, and notably, in semiconductor crystal growth. The vertical Bridgman crystal growth method is commonly used for growing semiconductor crystals such as leadtintelluride. Vertical Bridgman growth involves a liquid semiconductor melt layer, which lies atop the solid crystal phase being formed by the solidification of the melt. The arrangement is heated from above so that the crystal at the bottom of the system grows upward (see Figure 113). Since solidification is occurring at the liquidsolid interface, the temperature at that location is constant. Of course, a system that is heated from above, such as this crystal growth system, will not be destabilized by the vertical temperature gradient in the system. Still, buoyancydriven convection is possible in this system for another reason. As solidification occurs at the solid liquid interface, some portion of the species that comprise the liquid melt is rejected into the liquid phase, near the interface. The species comprising a liquid melt, such as leadtintelluride, are not all equal in weight. If the species rejected the most into the liquid near the interface is relatively light in weight, and makes the liquid region near the interface less dense than the remainder of the liquid melt phase above the interfacial region, then a topheavy system results. Thus, due to solutal gradients, buoyancy forces can cause convection in the crystal growth system just like the systems previously described in this chapter. Convection in this system can cause impurities from the ampoule to be transmitted to the solidliquid interface and this can cause flaws on the face of the forming crystal. Also, convection in the system can cause the heat transfer at the solidification interface to be non uniform. For these reasons, it is desired to know the conditions corresponding to the onset of convection in this system so that it can be better understood, predicted, controlled, and prevented. It should be mentioned that, in some cases, convection in a crystal growth system may not be a problem, and may even be desirable. If a crystal is grown at relatively low temperatures (unlike, for example, the high value of 1250 C at which gallium arsenide is grown), then convection helps transfer heat near the solidification interface and relieves thermal stresses there, preventing the growing crystal from fracturing. Also, the transport of impurities from the walls of the crystalgrowth container to the solidification interface does not tend to become a problem in these lowertemperature systems (in hightemperature crystal growth systems, it can relate to constituents being released from the crystal growth container by the degradation of the container). Thus, in some lowtemperature crystal growth systems, convection could clearly be beneficial. In vertical Bridgman crystal growth, the depth of the liquid phase is continuously changing as the melt solidifies into the growing crystal, consuming a portion of the melt. As the phase depth changes during crystal growth, so does the possible convective flow behavior. Thus, being able to predict the convective behavior of the system at differing fluid depths is important. Some crystal growth processes, such as the liquidencapsulated growth of gallium arsenide, require very high temperatures. In this research, convective behavior will be studied at more moderate temperatures, as the physics driving the convection are the same regardless of the temperature of the system. Typically, crystal growth is conducted using cylindrical fluid layers. It is interesting to ask how the convective behavior of a crystal growth system would differ if annular fluid layers were employed rather than cylindrical layers. One might imagine that the presence of an extra solid wall at the center of the fluid layer would, by means of the noslip condition it imposes on velocities at that location, impart to the crystal growth system a greater stability to convection. The additional center wall in an annular system would allow the system to suppress more disturbances than a cylindrical system. Furthermore, if the crystal were grown in an annular container with a sufficiently large annular centerpiece (relative to the outer annular radius), then any convective patterns which did arise would be very predictable and uniform sets of only azimuthally aligned cells, which could be easier to work with than the patterns that arise in cylindrical containers. Tcold zLX X Higher Density Fluid Layer Lower Density I hot Figure 11. Simplified system diagram: heated from below. z Lz Z 0 Figure 12. Rayleigh convection. Upliase .......... ............... "H O" T 6* 11 IIIIIIIIIIIIIIIIIIIIIlllllllllliiI . I", Phase I I Figure 13. Marangoni convection. Wavelength Figure 14. General diagram: critical Rayleigh number vs. disturbance wavelength. z t_O Two Convective Rolls zt Four Convective Rolls Figure 15. Crosssectional velocity profiles of flow patterns with different numbers of convective rolls. Racrit t I ' I I ' ~I II III Axz Figure 16. General stability diagram for buoyancydriven convection in a rectangular, laterally bounded system, with periodic boundary conditions at lateral walls, for varying aspect ratio. t% 4 x x zt x Four Convective Rolls Figure 17. Crosssectional velocity profiles of fourroll flow patterns with nonuniform roll size. 1708 Axz Figure 18. General stability diagram for buoyancydriven convection in a rectangular, laterally bounded system, with nonperiodic boundary conditions at lateral walls, for varying aspect ratio. m=0 @ D m=2 1 m=3 a, Figure 19. Diagrams of m 0, 1, 2, 3 flow patterns: "U": upward flow, "D": downward flow. Critical Ra vs. Aspect Ratio, Cylinder ....... m = 2 ::::::::::::::::::::::::::::::: m = 3 Figure 110. General stability diagram for buoyancydriven convection in a cylindrical container, for varying aspect ratio. Figure 111. Convective flow pattern in an annulus (Stork & Miuller 1974). 2900 2700 2500 2300 2100 1900 1700 0.7 m = 0 1.2 1.7 2.2 2.7 3.2 Aspect Ratio  m= 1 Example Flow Pattern in Annular System 0.01 0.008 0.006 0.004 0.01 0.01 0.008 0.006 0.004 0.002 0 0.002 0.004 0.006 0.008 0.01 xdirection 0.12 0.1 0.08 0.06 0.04 0.02 0 0.02 0.04 0.06 0.08 Velocity (scaled) Il.i i Velocities are Downward in zdirection and Positive Velocities are Upward in zdirection Figure 112. Example: computed flow profile for annular system, crosssectional view. Hot Liquid Melt Solid Crystal Cold Figure 113. Diagram of crystal growth system. CHAPTER 2 LITERATURE REVIEW Rayleigh convection in a single fluid layer has been a wellstudied problem in fluid mechanics for many years, and it has branched out into many newer problems over the years. Systems with multiple fluid phases, in which Marangoni convection exists in addition to Rayleigh convection, have been studied, as well. Also, researchers have examined differences in convective behavior seen in systems of differing shapes and sizes. It would be quite difficult to give a complete history of tremendous volume of research that has been done in this general field. In this section, a brief summary of some of the past research that may be considered most relevant to this study will be given. In Benard's classic experiments (1900), cellular, hexagonal flow patterns in vertically thin fluid layers (about .5 mm to 1 mm deep), heated from below, and open to air at the top surface were found (Benard 1900). Inspired by these experiments, about fifteen years later, Rayleigh investigated the instability that arises in a fluid layer heated from below (Rayleigh 1916). Originally, these researchers attributed this flow instability to buoyancy forces the phenomenon which would be termed Rayleigh convection. However, it turned out that the cells originally observed by Benard were caused more by surfacetension gradients than buoyancy differences. Also, it was found that the velocity and heat transfer boundary conditions on the fluid layer, as well as the size and shape of the system, had quite a significant impact on the critical conditions for the onset of convection. The early works of Benard and Rayleigh inspired many more researchers to investigate the convective behaviors of fluid layers. 2.1 Single Fluid Layers Single fluid layers can be studied in the presence of the Marangoni effect or in its absence. The former is more relevant to the computations and experiments in this thesis. Some of the earlier works in this area were that of Schmidt and Milverton in 1935 (Schmidt & Milverton 1935) and that of Silveston (Silveston 1958) in 1958. Both of these works involved experiments in which a single layer of liquid was heated from below and the critical vertical temperature difference for the onset of convection was determined. Schmidt and Milverton determined the critical vertical temperature difference corresponding to the onset of convection by measuring the heat transfer through the liquid layer, while increasing the vertical temperature difference, and then noting vertical temperature difference at which the heat transfer increased to a value higher than that of the initial state of pure conduction; this increase indicated the added heat transfer which accompanied the onset of convection. Silveston determined the critical vertical temperature difference in experiments not only by heat transfer measurements, but also by optical observations. In 1961, in his book, Hydrodynamic andHydromagnetic Stability, Chandrasekhar made a thorough mathematical analysis of the instability that arises solely due to buoyancy in a fluid layer heated from below (Chandrasekhar 1961). In 1967, Davis published a wellknown computational study of the Rayleigh convection that can occur in a bounded, threedimensional, rectangular box (Davis 1967). He used a Galerkin numerical scheme to model the behavior of the system, and determine the effects of lateral walls on the convective behavior. He also observed the widely known effect of greatly increasing the horizontal dimensions of the container, which is to cause the critical Rayleigh number to rapidly decrease to 1708. Lastly, as Chandrasekhar had done in his previously mentioned 1961 book, Davis explained that moderatelysized convective cells are preferred at onset to tall, narrow cells or wide, flat cells; this is because tall, narrow cells dissipate too large an amount of energy, and wide, flat cells require fluid particles to travel too great a horizontal distance before they can fall and release their potential energy. This sort of reasoning was introduced in Chapter 1. In a 1970 paper (Charlson & Sani 1970) and an extension of it published in 1971 (Charlson & Sani 1971), Charlson and Sani did a mathematical study to determine the critical temperature difference for the onset of axisymmetric (1970) and non axisymmetric (1971) flow patterns in cylindrical fluid layers heated from below. Their study included quite a wide range of aspect ratios (ratio of container radius to container height) and addressed the cases of both insulating and conducting side walls. In 1972, Stork and Muller published an experimental study in which they studied a system quite similar to the one mathematically analyzed by Davis in 1967 (Stork & Muller 1972). They varied the lateral dimensions of their rectangular system and found their critical temperature difference results to be generally below those of Davis. They attributed this discrepancy between experiment and theory to the effect of imperfect experimental control of the heat transfer boundary conditions at the lateral walls. In Koschmieder's 1974 experimental study, the heat flux through an oil layer before and after the onset of convection was monitored; he was, however, unable to get an accurate measure of the heat transfer through the oil layer for cases with low critical, vertical temperature differences (Koschmieder 1974). Like this research, the 1974 experimental study of Stork and Muller addressed Rayleigh convection in annular systems. As in their 1972 work, they examined a fluid layer of 10 mm vertical depth, comprised of silicone oil mixed with aluminum tracer powder, with the temperature at the top and bottom of the layer carefully controlled. In some experiments, they considered simply a system of cylindrical crosssection, while in others, they added a center piece to the cylinder to create an annular system. In all of their experiments, they detected the onset of convection by visual observation of the fluid layer and then obtained the critical vertical temperature difference by reading appropriately placed thermocouples. For the experiments on cylinders, onset conditions were determined for different aspect ratios; for the experiments on annuli, the onset conditions were determined for several different annular gap widths. As mentioned earlier, their experiments on annuli resulted in flow patterns that looked rather like spokes of a wheel, with several rolls of convecting fluid (always an even number of them) aligned in the azimuthal direction. As the annular gap width was increased (meaning, for example, the centerpiece of the annulus was made to be of smaller diameter), some quite interesting flow behavior was observed; for some large enough gap widths, there were more flow cells along the outer wall of the annulus than the inner wall. One important difference between the research in this thesis and the work of Stork and Muller is that this research focuses on the transition in the types of onset flow patterns that are observed when the annular gap width becomes sufficiently large (see Section 1.4), whereas this transition seemed to be more of a secondary observation in the work of Stork and Muller. In this research, annular systems of certain dimensions were designed specifically for investigating this transition in onset flow patterns. Sets of annular dimensions appropriate for this investigation were able to be chosen based on computations for the critical conditions, run for annular systems. Another important difference between this research and that of Stork and Muller is that this research includes a computation to accompany the experiments, and their work did not. In 1990, Hardin, Sani, Henry, and Roux published a computational study in which they determined the conditions at the onset of Rayleigh convection for cylindrical systems of several aspect ratios (Hardin, Sani, Henry, & Roux 1990). Also in 1990, Littlefield and Desai published a theoretical study in which they found the critical conditions for Rayleigh convection for an annular system with conducting side walls and with top and bottom walls that were assumed to be flat and free (Littlefield & Desai 1990). This research differs from the work of Littlefield and Desai in multiple ways. In this research, the computations address the more realistic case, with regard to a singlefluidlayer system, in which there is a noslip condition on velocity at the walls; additionally, the radial walls in this research were considered as insulating rather than conducting, and the variation of viscosity with temperature was accounted for. Another important difference of this research from the work of Littlefield and Desai is that it includes experiments which may be compared with the computations. The calculations of Littlefield and Desai matched qualitatively well with the experimental study of annular systems done by Stork and Muller. Like the experimental study of Stork and Miller, their work showed that as the ratio of the inner radius to the outer radius becomes larger (meaning that the annular gap becomes narrower), more and more convective cells form in the azimuthal direction at onset (larger azimuthal wave number at onset). They reached the conclusion that as the ratio of the inner radius to the outer radius (the "radius ratio" which was earlier named S for this research) becomes larger, the critical Rayleigh number (and critical vertical temperature difference) approaches the value corresponding to a vertical channel; they also determined that as the radius ratio becomes smaller, the critical Rayleigh number approaches the value corresponding to a vertical cylinder. Littlefield and Desai noted that the convective patterns which develop in the annular system have a higher tendency to include new cells in the azimuthal direction than in the radial direction. This can be explained, they say, by the physical argument that new cells which form in the azimuthal direction are more uniform in size (compared to one another and to cells already present) than new cells which form in the radial direction, and thus they form more easily because of the lower velocities and viscous stresses that they involve. At this point, to be clear, a list of the ways in which this research differs from the previous works on Rayleigh convection in singlefluidlayer annular systems is given. > This research includes corresponding computational and experimental work. > In this research, dimensions of annular systems investigated computationally and experimentally were chosen particularly to focus on the transition in the types of onset flow patterns which occurs when the ratio of the inner annular radius to the outer annular radius is sufficiently small, and to determine what value for this ratio corresponds to the transition (when this ratio is small enough, the onset flow patterns can have one or more convective rolls in the radial direction, while they generally have no convective rolls in the radial direction when this ratio is larger). > Computations for this research include the variation of viscosity with temperature. > Computations for this research address the realistic case in which there is a no slip condition on velocity at the walls. > Computations for this research are for the case in which the radial walls are insulating. Another study of Rayleigh convection in an annular system was presented in 1990 by Ciliberto, Bagnoli, and Caponeri (Ciliberto, Bagnoli, & Caponeri 1990). Their study was experimental and it involved observing the convective behavior of an annular layer of silicone oil by shadowgraph. The flow patterns observed in their study consist of an azimuthally aligned series of cells, and are consistent with the patterns seen by other researchers. In 1995, Zhao, Moates, and Narayanan published an interesting study of Rayleigh convection in cylindrical systems (Zhao, Moates, & Narayanan 1995). The study included both theoretical and experimental parts, which agreed well with one another. The experimental apparatus used was similar to the one being used in this research. The bottom temperature of a layer of silicone oil was controlled by circulating hot water below an aluminum plate, while the top temperature was controlled by circulating cooler water above a sapphire window. Flow was visualized by means of aluminum powder which was mixed in with the silicone oil. It is noteworthy that this study addressed the changes in convective behavior which arise when one accounts for the variation of viscosity with temperature rather than simply assuming it to be constant. It was seen that, depending on the aspect ratio of the system, the variation of viscosity with temperature could have a significant impact on not only the critical vertical temperature difference for the onset of convection, but also on the flow pattern at onset. This can be seen more later in this paper. The fact that the computations in this research account for the variation of viscosity with temperature when determining the critical conditions for convection is another noteworthy feature of this research, which may set it apart from previous works on annular systems. 2.2 Multiple Fluid Layers As mentioned, the studies on single fluid layers, for the fact that they exclude the Marangoni effect, are rather more relevant to this research. Still, some works on multiple fluid layers can be helpful in understanding what goes on in problems involving buoyancydriven convection in single fluid layers, too. Buoyancydriven convection, after all, is generally occurring in problems with multiple fluid layers, as well, even though the systems may be more complicated. Thus, a brief summary of it will now be given, and so will be presented a more complete picture of what sorts of research have been done on convective phenomena over the years. Following experimental evidence obtained by Block in 1956 (Block 1956), Pearson, in 1958, produced a paper which proposed a mechanism to attribute convective flow of the type observed by Benard to surface tension gradients rather than buoyancy forces (Pearson 1958). Pearson's study focused on a fluid layer, heated from below, laterally unbounded, with a rigid wall as its bottom surface, and a free surface at the top. A passive gas was the upper phase in the study. In that study, the temperature of the bottom wall was held fixed, while the temperature at the upper surface was governed by a heat transfer boundary condition. At the top and bottom surfaces of the fluid layer, Pearson considered the effects of having both conducting and insulating behaviors with regard to heat transfer. In Pearson's paper, a dimensionless number called "B" showed up as a preliminary form of what would, in later years, be called the Marangoni number. The Marangoni number is surface tensiondriven convection's analog to the Rayleigh number of buoyancydriven convection. Pearson's analysis assumed that the upper free surface could not deflect. In 1959, a paper by Sternling and Scriven presented physical mechanisms and a corresponding simple mathematical model to explain how several types of interfacial flows may develop in bilayer systems (Sternling & Scriven 1959). A 1964 paper by the same two authors extended Pearson's work by considering the same system but allowing deflection of the free surface (Scriven & Sternling 1964). Like Pearson's work, the work of Sternling and Scriven involved solving for the behavior of the system in response to a small, imposed disturbance. None of the works, that were just described, by Pearson, or Sternling and Scriven included an experimental component. The 1967 work of Koschmieder, however, was an experimental study of the convective behavior of a cylindrical layer of silicone oil, heated from below by a solid boundary, and cooled from above by a layer of air (Koschmieder 1967). In that study, the flow was visualized using aluminum powder, and it was found that convective flow in the system initially took the form of concentric circular rolls, and subsequently transformed into a hexagonal pattern. Koschmieder was able to make an accurate determination of the wavelengths of the flows, as well as their dependence on the dimensions of the fluid layer. Many studies involving multiple fluid layers were performed on systems comprised of one vapor layer atop one fluid layer. It is noteworthy that in 1972, Zeren and Reynolds published a study in which they examined Rayleigh and Marangoni convection in a liquidliquid system (Zeren & Reynolds 1972). The liquids considered were benzene and water. Their study included both theoretical and experimental parts. They felt that the presence of contaminants in the interfacial region could have affected their experimental results significantly. The work of Ferm and Wollkind, published in 1982, improved and extended the theoretical analysis given by Zeren and Reynolds (Ferm & Wollkind 1982). While it seems that most of the work done with annular systems has involved only single fluid layers, Bensimon, in 1988, published a study in which he experimentally examined the convective behavior in an annular layer of liquid with free boundary conditions at the top and bottom surfaces (Bensimon 1988). He arranged this system by placing a layer of silicone oil on top of a layer of mercury and leaving an air layer above the silicone oil. Flow visualization was accomplished in his study by the shadowgraph technique. In the shadowgraph images he presents in the paper, the flow patterns look quite similar to those that are seen in singlefluid layer, annular systems. Returning to the works of Koschmieder, an entirely experimental work was published by Koschmieder and Prahl in 1990 (Koschmieder & Prahl 1990). It investigated the tendency of wide fluid layers, when heated from below and open to air at the top surface, to convect in a pattern of hexagonal cells. They examined silicone oilair bilayers in small containers of differing shapes to find out what effects the differing container shapes had on the convective patterns that would form. They observed that as the widthtoheight ratio for the containers increased, more and more convective cells would form at the onset in order to fill the larger width. While determining the critical conditions by observation of the vertical heat transfer through the system is a nice method in terms of objectivity, Koschmieder and Prahl point out that in their experiments with small fluid layers, the voltage created by the heat sensor was too small in comparison with outside electrical noise to be useful. Thus, the more desirable method (which is the method used in the current research) was to optically observe the fluid layers and note the first appearance of fluid motion as the critical point. A 1992 experimental work by Koschmieder and Switzer examined surfacetension driven convection using shadowgraphy (Koschmieder & Switzer 1992). In a 1995 paper by Zhao, Wagner, Narayanan, and Friedrich, a theoretical study is presented that addresses Rayleigh and Marangoni convection in fluid bilayer systems both liquidliquid, and liquidvapor (Zhao, Wagner, Narayanan, & Friedrich 1995). A wide range of cases are considered, including heating from below, heating from above, and the case in which solidification is occurring at the bottom surface of a liquidliquid bilayer; the last case mentioned has strong similarity to crystal growth because of the solidification which it accounts for. The 1996 theoretical study of Dauby and Lebon addressed the Rayleigh and Marangoni effects in a liquidvapor bilayer system and included both linear and nonlinear analyses (Dauby & Lebon 1996). The linear analysis allows one to determine the critical vertical temperature difference at which convection begins, as well as the flow pattern at the onset. The nonlinear analysis provides more specific information on the flow behavior and can be used to predict supercritical behaviors. Most theoretical studies, up to this point, had assumed laterally unbounded systems. That assumption allows separation of variables, and thus an easier solution to the system of differential equations that model the problem. Of course, in actual experiments and other applications of these phenomena, lateral walls are present, and these walls can have important effects on the critical conditions and onset flow patterns. Thus, another point of interest in the study by Dauby and Lebon is that they formed their mathematical model to include noslip walls at the lateral boundaries. This is the most realistic boundary condition that could be enforced at those locations. Dauby and Lebon found their theoretical results to be in good qualitative agreement with the experimental results of Koschmieder and Prahl. In 1997 by KatsDemyanets, Oron, and Nepomnyashchy published a study, which examines convective behavior in trilayer fluid systems, and thus is also quite applicable to the science of crystal growth (KatsDemyanets, Oron, & Nepomnyashchy 1997). Johnson and Narayanan published a tutorial in 1999 which explained five different mechanisms by which convection in two vertically stacked fluid layers could couple (Johnson & Narayanan 1999). Their tutorial also discussed the situation in which, at certain aspect ratios, due to the effects of side walls and the dimensions of the system, a convecting system may find two flow patterns equally favorable from an energy standpoint; as a result, the system would oscillate between the two patterns. A point such as this is called a codimensiontwo point. CHAPTER 3 MODELING EQUATIONS To predict the onset conditions of convection requires a model that respects the physics of the problem. Such a model would utilize momentum, mass, and energy equations. The inputs to the model would be the vertical phase depth and thermophysical properties of the fluid, and the outputs would be the critical temperature difference needed for convection, as well as the associated pattern of the flow. The details of such a model are given in this chapter, while later in this thesis, an explanation of the numerical scheme used to solve the system of modeling equations is given. In Appendix C, the modeling equations are presented and developed more fully for the simple case in which the fluid is assumed to have a viscosity constant with respect to temperature. The modeling equations presented and developed throughout this chapter and Chapter 4, however, are for a more physically accurate mathematical model, in which the variation of viscosity with temperature is included. To get started the nonlinear equations that govern flow in a convecting layer of fluid shall be introduced. These equations are nonlinear primarily on account of the dependence of velocity on temperature, coupled with the effects of the interactions between velocity and temperature fields on heat transport. 3.1 Nonlinear Equations The domain equations used to model the convective behavior of the system are the momentum equation, the energy equation, and the continuity equation, which, respectively, are p +pW =VP+pg + V.S, (3.1) at pC +pCi VT= kV2T, (3.2) at tp V at (3.3) In these equations, p is the density of the fluid, t is time, i is the velocity of the fluid, P is the pressure of the fluid, g is gravity, S is the stress tensor for the fluid, Cv is the constant volume heat capacity of the fluid, Tis the temperature of the fluid, and k is the thermal conductivity of the fluid. The last term of the momentum equation has not been further simplified, at this point, because it will be needed in this form in order to properly account for the dependence of the fluid's viscosity on temperature. The stress tensor, S, can be expanded as s + (Vw ^ Vv~f~v (3.4) The use of this expansion, however, will be postponed for now. In this expansion, p is the dynamic viscosity of the fluid and r represents a transposed matrix. Since the velocity gradients in this system will not be very large, the energy equation has been simplified by neglecting the viscous dissipation. The system of equations is analyzed in cylindrical coordinates, which are clearly the natural choice for this problem; thus "r" denotes the radial direction, "0 denotes the azimuthal direction, and "z" denotes the vertical direction. The system may be treated as periodic in the 0direction (this is accomplished mathematically by the imo expansion of the variables into modes), with the oscillations of the form e in which m is the azimuthal wave number. Another approximation is to be made, and it affects the momentum equation. It is called the Boussinesq approximation and it addresses the variation of density with temperature. This approximation says that variation of density with temperature is negligible insofar as the change in momentum or mass is concerned, but that it does affect the acceleration in the system insofar as the body forces are concerned. The reasoning for introducing the Boussinesq approximation, and some notes on its applicability are explained in Appendix B as Section B. 1. Once the Boussinesq approximation is applied to Equation 3.1, it becomes R + R9 VV = VP+ pR(1 a(T TR)) + V. S. (3.5) In Equation 3.5, TR and PR are a reference temperature and the fluid's density at the reference temperature, respectively, and a is the volumetric thermal expansion coefficient of the fluid phase. Note, also, that the Boussinesq approximation results in the disappearance of the timederivative term of the continuity equation. At this point, multiple nonlinearities appear in the system of equations; a note about these nonlinearities is given in the appendix as Section B.2. An additional modification to the momentum equation is to define a modified pressure, which is  Vp = VP + pRg. (3.6) Using this, the momentum equation may be written as PR +PR V Vp pR(T TR )g + V.S. (3.7) at Now, the expansion of S shall be introduced to the momentum equation. The result of doing so is PR + PRV Vv = VP PRa(T R) + /V2V + 0t (3.8) As mentioned, a goal of this calculation is to include the effects of the variation of viscosity with temperature. Thus, the viscosity at any location can be written as S=/[f (T TR )], (3.9) wherefis some function of the difference between the temperature at that location and the reference temperature. In this equation, /R is the dynamic viscosity at the reference temperature. Note that, since p has the same dimensions as /R mustt be dimensionless. Using this expansion, Equation 3.8 may be rewritten as PR +p PRV*V =VppRa(TTR)g+Puf V32 + at (3.10) PR(Vf)(VV + (VVY) . The assumed form for viscosity's temperature dependence will be exponential, like the form assumed in the work of Zhao, Moates, and Narayanan (Zhao, Moates, & Narayanan 1995). The form of this exponential temperature dependence is f(T TR) eB(TTR), (3.11) in which B is a constant that can be determined using measurements of the fluid's viscosity obtained at a range of different temperatures. For temperatures in degrees Celsius, the dimensions of B would be  oC In this system of three domain equations, it appears that five unknowns are present. These are the velocities in the r, 0, and z directions, the pressure, and the temperature. The continuity equation can eliminate pressure by representing it in terms of velocities, and this reduces the number of unknowns to four. To solve for the convective behavior of a cylindrical system, constraints would be needed on the velocities in the r, 0, and z directions, as well as on the temperature, at the top wall (z = Lz), the bottom wall (z = 0), and the outer radial wall (r = Ro). Considering an annular system, however, four additional mathematical constraints are required; these constraints on the three components of velocity and temperature are applied at the inner radial wall of the system (r = R,). Figure 31 shows the geometry being considered. The symbol Lr is used to denote the radius in a cylindrical system, and also to denote the difference between Ro and R, (Lr = Ro R,) for an annular system. This means that for an annular system, Lr is the width of the annular gap that is filled with fluid. Thus, for the annular system, a total of sixteen constraints are needed, while only twelve are needed for the cylindrical system. These constraints are the boundary conditions. In their unsealed form, they are Vr =Vo =Vz =0 atz =0, (3.12) T = Tb at z= 0, (3.13) Vr = Vo = Vz = 0 atz=z (3.14) T=T, atz =L, (3.15) Vr V = Vz = 0 atr =R, (3.16) OT = 0 atr =Ro, (3.17) Ar and, when considering an annular system, the conditions at the inner annular wall, Vr =v =vz =0 atr =R,, (3.18) OT =0 atr =R, (3.19) 9r are needed. In these equations, vr, V0, and vz are the components of velocity in the r, 0, and z directions, respectively. The boundary conditions on velocity enforce noslip behavior at the walls, and state that no flow may pass through the walls. The conditions on temperature at bottom and top walls keep the temperature fixed at those locations; in these equations, Tb and T, are the constant temperatures at the bottom wall and top wall, respectively. The conditions on temperature at the radial walls represent the insulating behavior at those boundaries. 3.2 Scaling Now that the complete set of domain equations and boundary conditions has been listed, these equations will be made dimensionless. In scaling these equations, one has the option to choose several combinations of parameters to use in defining a characteristic velocity (V) and a characteristic time (T). A brief discussion of how to choose these parameters is in Section B.3 of Appendix B. Lengths are scaled using the vertical depth of the fluid. With that said, a list of the scaling relations used will now be presented. In these equations, the "hat" symbol indicates a dimensionless variable, and the "bar" indicates a characteristic value for a variable. The scaling relations are v=, (3.20) v T=T T, (3.21) AT AT =Tb Tt, (3.22) P, (3.23) p =R, (3.24) L t (3.25) t V = LV. (3.26) In the above equations, Lz is the vertical depth of the fluid phase. As mentioned earlier, TR is the reference temperature. In simple cases where the system of equations was solved without considering the dependence of viscosity on temperature, TR could be chosen to be the temperature at the top boundary of the system (Tt), and the reference dynamic viscosity, / R, would simply be the dynamic viscosity at that temperature. When the variation of viscosity with temperature is considered, however, TR is selected later, during the numerical solution of the system of equations, to correspond to the mean dynamic viscosity value of the fluid in its motionless base state, when the fluid layer is at the critical vertical temperature difference, just before the onset of convection (the viscosity of the fluid varies from top to bottom along the motionless liquid layer corresponding to the vertical temperature gradient). When the variation of viscosity with temperature is being considered, this mean dynamic viscosity value is the reference dynamic viscosity, /IR AT is the overall temperature difference across the system. In addition to these scalings, f(Equation 3.11) should be reexpressed in terms of dimensionless temperature. This new version off which will be called F, is F(T) = eB(AT)T (3.27) Applying these scalings to the domain equations, the following dimensionless (though the "hat" symbol will now be discarded) equations are obtained (Equations 3.283.30): pRL2 PRLv aL2g(AT)T 2 PR LZ+ RLz v Vv = Vp z (AT) T + FV2v + R//R t /R U V (3.28) (VF). (Vf + (Vf)T) LOT+ Ly f VT =V2T, K/CTat /C (3.29) 0=Vv, (3.30) in which K, the thermal diffusivity of the fluid, is equal to the fluid's thermal conductivity divided by the product of its density and heat capacity. Unless otherwise indicated, all thermophysical properties in these equations are considered at the reference temperature. The subscript "R" has not been included on the symbols for all of these properties; however, it has been left on the symbols PR and /R because, in those cases, it arose from the expansions by which the density and dynamic viscosity of the fluid were defined. This is the reason, for example, that the subscript "R" is not included on the kinematic viscosity in Equation 3.28, even though that property is simply a ratio of the dynamic viscosity and the density taken at the reference temperature. Observe that if a substitution of the ratio of K to Lz were made for V, then the coefficient of the temperature term in Equation 3.28 would be the Rayleigh number. Since many possible values could be assigned to 7, depending on the relative values of system parameters, iV has been left in the set of modeling equations. Again, a discussion regarding the definition of T is in Section B.3 of Appendix B. As for the scaling of the boundary conditions, the conditions on velocity remain unchanged in appearance, as do the conditions on temperature at the radial walls. However, the conditions on temperature at the top and bottom walls appear slightly differently depending on whether or not the problem is being solved considering the variation of viscosity with temperature. The reader is asked to refer to Appendix C for more information on the case in which this variation is not being considered. In cases in which the variation of viscosity with temperature is being considered, the dimensionless temperatures at the vertical boundaries are dependent on TR. TR, though, is dependent on the vertical temperature gradient in the fluid, which means that it needs to be determined along with the critical vertical temperature difference for convection by an iterative approach (described in Chapter 5). Thus, for the nonconstant viscosity case, the temperature boundary conditions at the top and bottom walls must be left in a more general form for now. In this general form, they are T =b R at z = 0, (3.31) AT T T S=at z = (3.32) AT This completes the presentation of the scaled domain equations and boundary conditions. Next, the equations shall be simplified by linearization and the removal of time dependence. Cylindrical System Annular System r = Ro r = 0 Fluidr = Ri r = 0 r = Ro Figure 31. System diagram: cylindrical and annular systems. Figure 31. System diagram: cylindrical and annular systems. CHAPTER 4 LINEARIZED EQUATIONS In this chapter, the process of simplifying the set of modeling equations by linearizing them around a motionless base state, and then eliminating their time dependence shall be shown. Since the goals of the calculations in this research are to determine only the onset conditions for convection (critical vertical temperature difference and flow pattern at onset) as opposed to the flow behavior at supercritical conditions, a linearized model is sufficient for the analysis. A nonlinear calculation would be necessary if the supercritical behavior of the system were to be determined. Time dependence is being removed from the equations because a system is independent of time at critical conditions (also known as marginal stability). Linearization and the removal of time dependence are accomplished by introducing a series perturbation expansion to the variables, then applying a further expansion to the perturbed variables (expansion into modes), and then setting the inverse time constant in the expansion equal to zero. As mentioned earlier, the same process is presented for the less complex case in which the fluid is assumed to have a constant viscosity with respect to temperature is included in Appendix C. 4.1 Linearization The system of equations shall be linearized around the motionless base state of the system. In this motionless initial state, heat is transferred within the system only in the vertical direction and only by conduction. The form of the expansions which shall be used in this linearization, considering, for example, the expansion of velocity, is 0 + + 2 + ac =s=0 2! a 2 =0 1 3 3V (4.1) 3 a.3 6 = 0 Replacing the expressions in brackets, this may be rewritten as v = v0 + +2v2 + 3v3 +.... (4.2) 2! 3! The subscript "0" denotes values pertaining to the motionless base state, in which heat is transferred only by conduction, and the temperature and pressure gradients are only in the vertical direction. Clearly, then, f0 is equal to zero. This fact allows the cancellation of some terms in the linearized forms of the modeling equations. If the system begins to flow due to some perturbation, then the remaining terms in this series (v3, v2, etc.), which are the perturbed variables, represent the flow behavior of the system. The amplitude of the perturbation that transforms the base state into the flowing, "perturbed" state is called E; this perturbation is taken to be extremely small. Considering a small perturbation allows a conclusion on the stability of the system to be reached using linearized equations. If a system is unstable to a small perturbation, it will, of course, be unstable to larger perturbations. If a system is stable to a small perturbation, then no such conclusion may be drawn; in this case, the system could simply require a larger disturbance to become unstable. The linearized equations (comprised of terms that are first order in E) are an approximation of the behavior of the system and can be used to determine the conditions corresponding to the onset of convection and the pattern of flow at the onset. If it is desired to determine the magnitudes of the velocities in the convecting fluid, as opposed to simply determining their values relative to one another, then a nonlinear analysis is necessary. A disturbance of magnitude E is mathematically applied to the system by expanding each variable in the modeling equations with the form shown in Equation 4.2. The set of linearized equations is then obtained by collecting only the terms which are first order in E. The resulting linearized domain equations are PRL2 Vp aLg(AT) FV L 2 Ot L L aT V2T (4.3) L\0 __ V 2T (4 .4) ht7 t K Oz 0 = V.. (4.5) In these equations, Fo is the same as F except that the exponential temperature value is now the base state temperature, To. Fo(To), then, is Fo (T) = e B (AT) (4.6) The linearized forms of the boundary conditions are S= 0 = Vlz = 0 atz = 0, (4.7) T, =0 atz=0, (4.8) v1, = V =lz = 0 atz =L, (4.9) T, = 0 atz =Lz, (4.10) V1, = o 0 VIz = 0 at r = Ro, (4.11) S= 0 at r = Ro, (4.12) ar and, when considering an annular system, V1r =V0 = Vz = 0 atr =R, (4.13) aT 0 atr =R, (4.14) 9r hold. Linearization simplifies the boundary conditions on temperature at the top and bottom walls because it is the base state values of temperature, To at the top wall and bottom wall, which are equal to the fixed values of temperature at those boundaries. Substituting for Fo in the momentum equation, and making use of VFo = VT = g (4.15) aTo aTo az in which 3 is the unit vector in the zdirection, brings the momentum equation into the form PRL Z _v a L + zg(AT) B(AT) + PtRt 8t U V T (4.16) (B(AT)eB(AT)T O z ), *(Vf1 + (Vf) 8z Note that all factors on the righthandside of Equation 4.15 are dimensionless since VFo, from which they arose, was already made dimensionless in Chapter 3. For a more detailed and well given explanation of this perturbation method, refer to the book titled Interfacial Instability by Johns and Narayanan (Johns & Narayanan 2002). The momentum equation (Equation 4.16) is now rewritten as the three scalar equations which are its components in the r, 0, and z directions. The V operators are expanded in cylindrical coordinates from this point onward. The scalar component equations comprising the vector form of the momentum equation are the rcomponent of the momentum equation, the 06component of the momentum equation, and the zcomponent of the momentum equation, which, respectively, are vt t Or 1 02 a2 r2 92 z2 (4.17) (B(AT)eB(AT)To LU Vt eB(AT)To ) + , )[r + r  r D0 1 a2 a2 r2 92 z2 T& & az az 17 I 2 Yio 1 &V r 00 (4.18) eB(AT)To I I, p ,  + vU t at 0z B(AT)TO 02 1 a 1 02 02 e (++ +_ Vl + Or 2 r Or r2 O 2 16 1 S ] 1 (4.19) (B(AT)e )(A) T) 2 + a g T, . 4.2 Expansion into Normal Modes The set of domain equations now includes Equations 4.4, 4.5, 4.17, 4.18, and 4.19. A second expansion (expansion into modes) will now be applied to the variables, as well. This new expansion separates the rdirection and zdirection dependencies of each variable from the 0  direction and time dependencies. In the example expansion shown as Equation 4.20, the new variable representing only the rdirection and zdirection dependencies is marked with a "prime" symbol. This expansion assumes a periodic form for the 06 direction dependencies, which is sensible, since the cylindrical and annular systems being considered are indeed periodic in the 06 direction. The nature of the periodicity in the 06 direction is described by an azimuthal wave number, called m. Note that this wave number is dimensionless since it arises from dimensionless variables. The time dependence of each variable is assumed to have an exponential behavior. This exponential behavior with respect to time is governed by an inverse time constant, called c c is dimensionless since it arises from dimensionless equations; if unsealed, its units would be (1/time). Using pressure as an example, the form of the expansion is = im + (420) /01 P 1 Z ) e (4.20) The form of 06 direction dependency assumed in the expansion can be used because the periodic spatial dependence of the system in the 0direction can be represented by a series of sines and cosines. When the system of modeling equations is subjected to a perturbation (E), the perturbation may either grow, resulting in the onset of flow, or it may die out if the system is not at critical or supercritical conditions. The value of a is an indicator of the system's response to a given perturbation. It should be noted that in certain cases, for which the system oscillates between convective flow patterns, the value of a is imaginary. It can be mathematically shown, however, that a is strictly real for the nonoscillatory cases considered in this research. If a disturbance applied to the system dies out, meaning that the system is stable, then a has a negative value. If the system is unstable, and will flow when subjected to a perturbation, a is positive. At marginal stability, however, when the system initially becomes unstable, the system is independent of time and the value of a is 0. The fact that a = 0 at marginal stability allows the elimination of all timederivative terms in the modeling equations once this expansion is applied. The final forms of the modeling equations (including all three components of the momentum equation, the energy equation, the continuity equation, and the boundary conditions), with this expansion applied, and with the "prime" symbols dropped from the newly defined variables, are (m2 + 1) (B(AT)eB(AT)To &0 & im 0 = p+ r eB(A)T)o a2 or2 2im r2 I 1 a +r 9 (B(AT)e(AT)To Dz B(AT)TO 2 e 1 r 2 r 1r r&r (B(AT)e B(AT) 8z a2 (2 (m2 +1) 2+ v + r2 le im r m2 2 2+ v + 2 2 1, 82 ag(AT) Dz V 0= + Or a2 or2 a2 + &2 r r 2im _ r2 1 1 (4.21) (4.22) (4.23) B(AT)To Sr O=L 9 To [ 92 1 9 92 m2 TTO V 0+ 2 2 T, (4.24) cK 9z 9r r 9r z2 r 8 ^ 1 im 8 0= + V1, + 1+ 1V, (4.25) Vlr, = Vl = Viz = 0 atz = 0, (4.26) T, = 0 at z = 0, (4.27) Vi. = V10 = Vz = 0 atz =Lz, (4.28) T, = 0 at z = L, (4.29) vl,. = V = ~ = 0 atr =Ro, (4.30)  = 0 at r = Ro, (4.31) Or and, when considering an annular system, v. = V0 = vz = 0 atr = R, (4.32) = 0 at r = R, (4.33) Or are needed. At this point, two important dimensionless groups which show up in the motion equation are Fo (see Equation 4.6) and Q = B(AT)eBATTO (4.34) Now, the entire collection of modeling equations (Equations 4.21 4.33) may be simultaneously, numerically solved on a computer as an eigenvalue problem, in which the eigenvalue is A T. Since c has been set equal to 0, the problem is being solved at the critical point for the onset of convection, and so the A Tvalue to be determined is the critical vertical temperature difference for convection, or [ A T]crit. How exactly the solution for the [ A T]crit is carried out is the topic of the next chapter. Notice that the Prandtl number (v/Ki) does not appear in any of these equations. This may be surprising since one might expect that ratio of thermophysical properties to play a role in determining the critical conditions for convection. Now, recall that in Chapter 1 it was shown that the graph of critical Rayleigh number versus aspect ratio, for different azimuthal wave numbers, was identical for any system, regardless of the system's thermophysical properties. This interesting fact can be clearly explained at this point. The kinematic viscosity and thermal diffusivity of the system affect how quickly disturbances may die out or grow. Thus, before the onset of convection, their values make c more or less negative depending on how stable the system is, and after the onset of convection, their values make U more or less positive depending on how rapidly a destabilizing disturbance grew and how strongly the system is convecting in response to it. Precisely at the onset of convection, however, which is the condition at which the modeling equations are being solved, the system is independent of time (o = 0), and so the thermophysical properties play no role in determining the critical conditions for convection. If, rather than simply changing the thermophysical properties of the fluid, the insulating boundary conditions on temperature at the radial walls were changed to conducting boundary conditions, the graph of the critical Rayleigh number versus aspect ratio would change (as explained in Chapter 1, the critical vertical temperature differences for convection are lower when insulating radial walls are considered rather than conducting radial walls). CHAPTER 5 SPECTRAL SOLUTION METHOD The goal of this chapter is to explain the spectral solution method used for computations in this research. The method is introduced and described in Section 5.1. Section 5.2 focuses on the application of the method. Rearranging the set of modeling equations to be solved numerically as an eigenvalue problem for A Tis not too difficult. To do so, one must isolate the term of the zdirection component of the momentum Equation 4.23 which includes T1. What is meant by isolating that term is arranging the equation so that that term is on the side of the equation opposite from all of the other terms. Supposing the side of the equation where the T, term is placed is chosen to be the righthandside of the equation, then the remaining terms of the zdirectioncomponent of the momentum equation, and all nonzero terms of all of the other domain equations and boundary conditions should be kept on the lefthandsides of those equations. In general, what is being done here is that the problem is being recast in the form of a generalized eigenvalue problem for A T. The form of this generalized eigenvalue problem is AX =ABX. (5.1) Here, A is the matrix containing the coefficients of the velocity components, temperature, and pressure from the lefthandsides of the domain equations and boundary conditions. B, as one would imagine, is the matrix of the coefficients of the velocity components, temperature, and pressure from the righthandsides of the set of modeling equations. X is the eigenvector; it is a columnvector including all three velocity components, temperature, and pressure. From this point onward, the subscript "l"'s used to indicate perturbed variables will be discarded, as nearly all variables referred to will be perturbed variables. The only exceptions to this would be the base state variables, which, for this reason, will still be denoted by a subscript "0" in all locations. The form of X is Vr Vo X= v. (5.2) T P A, of course, is the eigenvalue, which is A T. Arranging the system of modeling equations in this form is straightforward when the variation of viscosity with temperature is not being considered. For the nonconstant viscosity case, however, A Tis present in more locations than just the coefficient of T, in the zdirectioncomponent of the momentum equation; it is also present as an exponent in several terms. These terms, though, may be kept on the lefthandside of the modeling equations rather than being placed on the righthandside of Equation 5.1. If the equations are arranged in this way, A T can be determined by an iterative approach. In this iterative approach, a value must first be chosen for Tt, the temperature at the top of the fluid layer. In all experiments done for this research, this temperature was kept constant at 30.0 C, so that value was always substituted for Tt in the calculations. Next, a guessvalue, which will be called [A T]guess, is substituted for A Ton the lefthandside of Equation 5.1. Then, the eigenvalue (A T) on the righthandside of Equation 5.1 is solved for, and then its value is used to update the [A T]guess value on the lefthandside of Equation 5.1, and so on. The procedure is repeated until convergence, which occurs quickly. Since TR and '/R are dependent on the vertical temperature gradient prior to the onset of convection, their values are updated as [ A T]guess is updated. Note that updating the [ A T]guess value on the lefthandside of Equation 5.1 represents updating the value of viscosity, through its exponential temperature dependence. Also, note that in the case where the variation of viscosity with temperature is not accounted for, the viscosity value at 35.0 C, which is given in Appendix A, is used throughout the calculation. 5.1 Explanation of the Method Once the system of modeling equations has been written in the form of Equation 5.1, it can be used to numerically approximate the critical vertical temperature difference for convection, and the onset flow pattern. As is typical of numerical solution methods, the first step is to consider discretized versions of the modeling equations, which describe the behavior at a collection of individual points in the system as an approximation of the full behavior of the system. This means that each variable is written as a vector containing the values of that variable at a set of points within the system. For example, a cylindrical system could be discretized into sets of points (or "nodes") in the rdirection and zdirection as shown in Figure 51. Of course, the subscripted numbers here have nothing to do with the subscripted numbers that appeared during linearization in Chapter 4; they simply indicate the spatial arrangement of the nodes. In this case, the cylinder was discretized into rows of four nodes in the rdirection, and columns of four nodes in the zdirection. Clearly, these sets of nodes in the rdirection and z direction form a grid of nodes in which the location of any node can be specified with an index for the rdirection and an index for the zdirection. For this reason, the discretization nodes are also referred to as "grid points." Either direction could have been discretized into any number of nodes. Of course, considering a larger number of nodes leads to a more accurate discretized model. The numbers of nodes created in the rdirection and the zdirection are decided by the selection of parameters called Nr and N, respectively. For the example above, Nr = 3 and Nz = 3. Note that discretizing the rdirection and zdirection means that r and z are written as vectors now rather than scalars. In general, the sets of nodes in the rdirection and zdirection, respectively, can be written as ro F = r2 (5.3) rv z1 (5.4) ZNz Due to the inclusion of ro and zo, the numbers of nodes in the rdirection and zdirection are actually (Nr + 1) and (Nz + 1), respectively. It should be mentioned that, in this research, the treatment of the 0 direction was much simpler than the discretization applied to the rdirection and zdirection, due to the periodic form assumed for 06 direction dependencies; the treatment of the 0direction will be discussed later. For now, the reader is free to imagine that the discussion here pertains to simply systems axisymmetric in the 0 direction, even though, in actuality, the treatment of the 06 direction is such that everything presented here applies to systems of any azimuthal periodicity. As indicated by Figure 51, when considering a system in two or more spatial dimensions, the indices of the nodes in the rdirection and zdirection may be thought of as indices corresponding to different planes slicing through the system, normal to the rdirection and the zdirection, respectively. Thus, discretizing the rdirection into (Nr + 1) nodes the and z direction into (Nz + 1) nodes divides the system into a grid of nodes that includes a total of ((Nr + 1) (Nz + 1)) nodes. Again, note that the computational treatment of the 0 direction has not yet been addressed; this means that the discretized form of the system being discussed at the moment is only one plane in the 0 direction. Given this form of discretization, every variable in the system must be written as a vector of its scalar values at every node of the discretized system. Temperature, for example, would be written in discrete form as shown in Equation 5.5. To allow a shortened, more clear representation on paper, Equation 5.5 will show T in the transposed arrangement; note that it is actually a columnvector like F and f The transposed form of the discretized temperature is r =[T T IT I ... IT IT I . ro,zo r8,zo *' r,Zo ro,zz r,z >* *> T T I T 1 (5.5) rNr,z Z i ro,zz Trl,ZNz ,ZNrz (5.5) The raised ellipsis in Equation 5.5 indicates the continuation of this rowvector on the next line. The discretization nodes, although they may appear to be so in Figure 51, are not necessarily evenly distributed in the system. In fact, it can be quite beneficial to use nodes which are not evenly spaced. The sets of nodes used in this research are not evenly spaced; they are called the GaussChebyshevLobatto grid points and the GaussRadau grid points (Trefethen 2000). Both sets of grid points are more highly clustered near the boundaries of the system than near the center of the system. It is highly beneficial to use these sets of clustered grid points as opposed to evenly spaced grid points because the clustered grid points allow much quicker and more accurate numerical convergence. Using these sets of grid points is also beneficial because, in the systems considered in this research, much of the important convective flow behavior occurs near the edges of the system, and using sets of grid points more clustered near the edges of the system thus allows the convective behavior of the system to be more easily captured by the discretized modeling equations. GaussChebyshevLobatto grid points are a set of points spanning the range [1,1] (including the endpoints at both 1 and 1), and so the actual distances within the system must be rescaled in order to fit this range. The set of GaussRadau grid points, however, spans the range (1,1], meaning that it does include the endpoint at 1, but does not include the endpoint at 1. Simple diagrams showing the general appearance of the GaussChebyshevLobatto and Gauss Radau sets of grid points, considering the rdirection, and considering seven grid points, are shown below as Figure 52. To show how to determine the exact spacing of these clustered sets of grid points, the equations to generate the GaussChebyshevLobatto and GaussRadau points are shown below, considering the discretization of the rdirection. The equation to generate the GaussChebyshev Lobatto points is r =cos j = 0,1,...,N,, (5.6) TN) and the equation to generate the GaussRadau points is r1 = cos 2 j = 0,1,... Nr (5.7) 2Nr +1 I A negative sign has been added to the GaussChebyshevLobatto equation so that that set of points corresponds with Figure 51. For the GaussRadau points, r, has been marked with a superscripted "*" to indicate that the order of the set of points generated by Equation 5.7 actually needs to be reversed so that the points correspond with Figure 51. In all calculations done for this research, the GaussChebyshevLobatto grid points were used in discretizing the zdirection. In calculations for annular systems, the GaussChebyshev Lobatto grid points were used for the discretization of the rdirection, as well. In calculations for cylindrical systems, however, GaussRadau grid points were used in discretizing the rdirection. The reason for this will be explained now. The discretization in the rdirection is done along the radius of the cylinder, and not the diameter. In the set of modeling equations for the cylindrical system, though, there are no boundary conditions at r = 0 (the interior endpoint of the discretization in the rdirection). Thus, including that location in the discretization of the r direction would cause computational problems. This difficulty can be dealt with by excluding the location r = 0 from the discretization. The use of GaussRadau grid points for the discretization in the rdirection accomplishes this because this set of grid points excludes one endpoint. Given that all variables must be written as vectors of their scalar values at every node of the discretized system, in the manner shown in Equation 5.5, it is clear that the eigenvector, X, will actually be a concatenated columnvector comprised of the columnvectors for each system variable. This means that it is a vector in which the smaller columnvectors for each system variable are vertically arranged above and below one another, headtotail. In the discretized form of the problem, the matrices A and B of Equation 5.1 are actually large matrices comprised of many smaller matrices which operate on the discretized variables in X in accordance with the modeling equations. Each row of submatrices within A and B represents a different domain equation or boundary condition. Each submatrix in a row of submatrices representing a domain equation must operate on all grid points in the system, whereas each submatrix in a row of submatrices representing a boundary condition operates on only the grid points on the boundary where the condition is enforced. A diagram illustrating the general layout of the matrices and vectors of the discretized, generalized eigenvalue problem is shown below as Figure 53. All differential operators in the modeling equations must be expressed in matrixform so that they may operate on the discretized variables. The differentiation matrices can be derived from polynomial interpolation equations. When taking derivatives in a direction that was discretized with GaussChebyshevLobatto grid points, the differentiation matrix is different than the one that would be used for taking derivatives in a direction that was discretized with Gauss Radau grid points. An example differentiation matrix, for rdirection differentiation on Gauss ChebyshevLobatto grid points (as given by Equation 5.6) with Nr = 3 is Dr,N,=3 = 3.1667 4.0000 1.3333 0.5000 1.0000 0.3333 1.0000 0.3333 (5.8) 0.3333 1.0000 0.3333 1.0000 0.5000 1.3333 4.0000 3.1667 where the subscript of Dr reflects the fact that this matrix is for N = 3. More details about how to obtain differentiation matrices and how to set up and apply the spectral solution method being used in this research can be found in a helpful book by Trefethen, called Spectral Methods in MATLAB (2000). The discretization and discrete problem formulation described so far are an accurate description of how to treat the system of modeling equations on a given plane normal to the azimuthal direction (an rz plane). Dealing with the 0 direction dependence of the physical system is mathematically quite simple because of the periodicity of the system in the 0 direction. Recall that the periodic 0 direction dependence of the system was represented in the form shown in Equation 4.20. This results in the presence of the azimuthal wave number, m, in the final set of modeling equations. To account for the 06 direction dependence of the physical system, one must simply assume a value for m. To choose a value for m, of course, is to assume the exact periodic form (number of periods in the 0direction) of the convective flow pattern at onset. If, for example, m is assumed to be equal to 2, then the onset flow pattern being sought is one which repeats twice as one progresses a single cycle through the fluid layer in the 0  direction. Ifm = 0 is assumed, then the onset flow pattern is one which does not vary in the 0  direction at all; the m = 0 case is axisymmetric in the 06direction. Once a value is set for m, and thus the particular periodic form of the onset flow pattern has been assumed, the discrete system of modeling equations in matrixform (Figure 53) can be solved for the eigenvalue, A T. This can be done numerically with an eigenvalue solver (for example, in MATLAB) that accommodates generalized eigenvalue problems. Note that, in this research, the lefthandside matrix (the A matrix in Figure 53) was not of full rank. Another way of saying this is that the set of eigenvalues of the A matrix, itself, included some eigenvalues which were equal to zero. MATLAB was selected as the software to be used in performing the computations for this research. The MATLAB generalized eigenvalue solver used in this research, like most other generalized eigenvalue solvers that the author has encountered, required that the A matrix be of full rank. Consequently, if the matrices of discretized modeling equations were directly fed into a generalized eigenvalue solver, then the output would include several spurious values. The reasons for the appearance of the spurious values related to the redundancy of boundary conditions at some points. To read more about this type of complication, the reader is asked to see the 1993 work of Labrosse (Labrosse 1993). In order to eliminate spurious values from the output of the eigenvalue calculation, the A matrix was put through a preconditioning process involving singularvalue decomposition. In this process, the A matrix was decomposed into three matrices by singularvalue decomposition, then the rows in the matrices corresponding to any eigenvalues equal to zero could be located, and finally those rows were filtered out of further computations. The critical temperature difference and flow pattern at the onset of convection can be determined by comparing the A Tvalues (eigenvalues) obtained for a range of different m values. When making this comparison, the lowest A Tvalue corresponds to the m value for which the fluid layer is most unstable to buoyancy forces. The lowest A Tvalue and the corresponding m value, thus, represent the critical temperature difference and flow pattern at the onset of convection, respectively. Note, also, that in order to account for the variation of viscosity with temperature, it is necessary to use an iterative procedure (as described earlier in this chapter) to obtain the A Tvalue corresponding to each m value. A flowdiagram displaying the general procedure carried out in MATLAB to solve the system of modeling equations for onset conditions is given in Appendix D. A set of MATLAB programs were created to be used in the calculations for this research, and some examples of these programs are given in Appendix E. 5.2 Application of the Method As mentioned, considering a larger number of nodes makes the discretized model more accurate. If the number of discretization nodes used to model the physical system were small enough, an inaccurate critical temperature difference would be obtained from the eigenvalue calculation. If the number of nodes were, then, increased further and further, the critical temperature difference results given by the eigenvalue calculation would progressively converge toward the physically correct value. Eventually, when considering a sufficiently high number of discretization nodes, the critical temperature difference obtained would be the physically correct value, and further increases to the number of nodes considered would have no impact on the critical temperature difference result obtained. A nice quality of the spectral method used here, though, is that not too large a number of nodes are actually needed for the computational solution to converge to an accurate solution. In all computations done for this research, 17 nodes were used in the rdirection and in the zdirection; this was always more than enough nodes to obtain a converged solution. Convergence of the computational results was, of course, verified by repeating computations with even higher numbers of discretization nodes. To exemplify this convergence behavior, a table is given below (Table 51) which shows the convergence of critical temperature difference values obtained using the author's MATLAB program. The example system considered for this convergence table is a fluid layer in a cylindrical container, with rigid, noslip walls, insulated at the lateral boundary, and heated from below so as to facilitate Rayleigh convection. In particular, a case is being considered in which the size of the cylindrical system is such that the convective flow pattern at onset is a onecell pattern, axisymmetric in the 06direction. The convergence of the critical temperature difference ([ A T]crit) is shown with respect to both Nr and Nz. To confirm the validity of the set of computational programs written for this research, results obtained using the programs were checked against results for standard RayleighBenard problems (problems involving buoyancydriven convection in systems heated from below) and some results obtained by other authors. As one check of the validity of the computational programs, a comparison was made between the critical temperature difference results obtained for a twodimensional, rectangular case (xdirection and zdirection only) and a onedimensional, rectangular case (zdirection only). In the twodimensional case, the rectangular fluid layer was heated from below, with stressfree, insulated side walls, and constanttemperature, noslip walls as the top and bottom boundaries of the layer. The distance between the side walls for the twodimensional case is called Lx. In the onedimensional case, the fluid layer was heated from below, was bounded at the top and bottom edges by constanttemperature, noslip walls, and was unbounded in the lateral direction. In neither case was the variation of viscosity with temperature considered. The twodimensional case with stressfree side walls was expected to produce exactly the same critical temperature difference results as a corresponding onedimensional case, in which the system had the same vertical depth as the twodimensional case. As Table 52 shows, it did. At first it may appear unusual that the onedimensional calculation could produce two different critical temperature difference results for the same vertical depth (Lz = 7.2 mm). The reason this is possible has to do with the way the onedimensional calculation is carried out. As explained in Chapter 1, a twodimensional system with lateral walls can only physically accommodate a certain set of onset flow patterns, which is dependent on the lateral width of the system. Thus, if the critical temperature difference result for a twodimensional system with lateral walls is to be compared with the result for a corresponding onedimensional calculation, it is only the disturbances of proper shape and size to induce these physically admissible onset flow patterns which may be considered in the onedimensional calculation. It is because of this that the critical temperature difference results obtained from the onedimensional calculation for the two cases in which Lz = 7.2 mm are different. Now that results are being shown for systems of specified depths and boundary conditions, it is a good time to mention that the fluid used in all experiments and calculations for this research is Dow Coming 200 1 Stoke silicone oil. Thus, its thermophysical properties are used in all calculations presented in this paper. A list of these thermophysical properties is given in Appendix A. The dependence of temperature on the viscosity of the oil was determined using viscosity measurements taken with a ColeParmer 98936 series viscometer. The exponential equation for the viscosity of the oil as a function of temperature is included in Appendix A, as well as a description of the manner in which viscosity measurements were obtained. Except for the viscosity, which was measured as described, the thermophysical properties of the silicone oil were based on information provided by Dow Corning. As an initial check of the ability of the computational programs to model cylindrical systems, critical temperature difference results were computed for buoyancydriven convection in a cylindrical fluid layer, and compared to the results given in a 1990 publication by Hardin et al. (Hardin, Sani, Henry, & Roux 1990). The particular system considered for this comparison was a layer of silicone oil in a cylindrical container, heated from below, with noslip walls at all boundaries, constanttemperature conditions at the top an bottom boundaries, and insulation at the radial wall. Again, the variation of viscosity with temperature was not considered. Hardin et al. give several critical temperature difference results for this system, as well as the corresponding flow patterns. A comparison of the results produced by the author's computations for this system with those of Hardin et al. is given in Table 53. These comparisons show the validity of the MATLAB programs written for the computations in this research. The results of Hardin et al. were given in terms of the dimensionless Rayleigh number, so in order to compare with them, the author's critical temperature difference results in this table have been reexpressed in terms of the Rayleigh number, as well (as shown in Equation 1.1). Since the variation of viscosity with temperature was not considered in this case, the viscosity value used in the computation was simply the constant viscosity value at 35 C given in Appendix A. The aspect ratio is simply the ratio of the radius of the cylindrical container to the vertical height of the container (L/ Lz). The way in which convective flow patterns are named should be briefly discussed here, too. In general, they are named by their azimuthal wave number, m. For example, a ringshaped pattern axisymmetric in the 0direction is called an m = 0 pattern. Optionally, when considering cylindrical systems, the flow pattern may be referred to with a parenthetical notation, including two indices. A parenthetical system for naming the onset flow patterns in annular systems is given in Chapter 7. In the parenthetical notation for cylindrical systems, the first index is the azimuthal wave number and the second index represents the maximum number of convective rolls that can be counted across the diameter of the cylindrical test section. Diagrams to exemplify the use of this parenthetical notation for flow patterns in cylindrical systems are shown as Figure 54. The diagrams are taken almost directly from the 1995 paper by Zhao, Moates, and Narayanan (Zhao, Moates, & Narayanan 1995). In Figure 54, a top view of each flow pattern is given, and below the top view of each pattern is a side view of the same pattern. In the top views, the "X" indicates falling fluid and the "0" indicates rising fluid. Diagrams (a) and (c) in Figure 54 could simply be called m = 0 patterns if the parenthetical notation were not being used; likewise, diagrams (b) and (d) could be called m = 1 patterns. Now, the method of obtaining computational results has been explained, and the validity of the computational method has been demonstrated. The computational results were compared with the results from experiments, and an explanation of the experimental apparatus and the manner in which experimental results were obtained is given in Chapter 6. Table 51. Example of convergence of computed result with Nr and Nz. Nr Nz [AT]crit(C) 3 13 6.55 4 13 5.55 5 13 5.60 6 13 5.59 7 13 5.60 8 13 5.60 9 13 5.60 10 13 5.60 11 13 5.60 12 13 5.60 13 13 5.60 Nr Nz [A Ticrit (oC) 13 3 Inf 13 4 5.42 13 5 5.55 13 6 5.60 13 7 5.60 13 8 5.60 13 9 5.60 13 10 5.60 13 11 5.60 13 12 5.60 13 13 5.60 Table 52. Comparison of rectangular, 2D, nostress results with rectangular 1D results. [ A T]crit (OC), L. (mm) L, (mm) 2D, NoStress [ A T]crit (OC), 1D 5 5 12.87 12.87 7.2 7.2 4.31 4.31 23 7.2 4.33 4.33 9 9 2.21 2.21 Table 53. Comparison of calculated cylindrical results with results of Hardin et al. Predicted Flow Pattern, Author's Computation Predicted Flow Pattern, Hardin et al. Critical Rayleigh Number, Author's Computation Critical Rayleigh Number, Hardin et al. Aspect Ratio: .75, Lr = 4.5 mm, Lz = 6 mm (1,1) (1,1) 2590 2592 Aspect Ratio: 1, Lr = 6 mm, Lz = 6 mm (0,2) (0,2) 2260 2260 Aspect Ratio: 1.5, Lr = 9 mm, Lz = 6 mm (0,2) (0,2) 1895 1895 Aspect Ratio: 2.5, Lr= 15 mm, Lz = 6 mm (0,4) (0,4) 1780 1781 ro, Z3 t ri, r2,Z3 2 3,Z3 r z >  SroZ 1,Z2 r2,Z2 3, Z2 ro,zi r,z r2,ZI r3,Z ^^__ ___  i    o ro,Zo rZo0 r2,Zo rZO Figure 51. Discretization nodes in a cylinder. GaussChebyshev =6 Lobatto  r .866 r5 .866 ro 1 r2 .5 r3 0 r4 .5 r6 1 GaussRadau Nr=6 r, .7485 r5 .8855 ro .9709 r2 .3546 r3 .1205 r4 .5681 r6 1 Figure 52. Example: grid point spacing. LeftHandSide Matrix (LHS), A LHS Motion Equation, rcomponent, and velocity boundary conditions LHS Motion Equation, 0  component, and velocity boundary LHS Motion Equation, zcomponent, and velocity boundary conditions LHS Energy Equation, and temperature boundary conditions LHS Continuity Equation rV T P RightHandSide Matrix (RHS), B Figure 53. Diagram of matrix/vector arrangement of discretized problem. RHS Motion Equation, rcomponent, and velocity boundary conditions RHS Motion Equation, 0  component, and velocity boundary RHS Motion Equation, zcomponent, and velocity boundary conditions RHS Energy Equation, and temperature boundary conditions RHS Continuity Equation Vr T P [OO (a) L, / L, = 1, Pattern: (0,2) o0D (b) L,/L = 1.8, Pattern: (1,3) (c) L,r / L = 2.5, Pattern: (0,4) (d) L, / L, = .75, Pattern: (1,1) Figure 54. Examples of flow patterns and parenthetical notation for cylindrical systems. CHAPTER 6 EXPERIMENTAL DESIGN After computational results were obtained in this research, they were compared with experimental results for systems identical to those considered in the computations. The apparatus which was designed to run experiments for this research is the subject of this chapter. Any possible errors associated with the apparatus and its measurements are discussed in Chapter 7. 6.1 Goals in Experimental Design There were several requirements made of the experimental system. One requirement was that the top and bottom temperatures of the fluid layer (T, and Tb, respectively) were kept constant and uniform at desired values. By monitoring and controlling the top and bottom temperatures of the fluid layer, the vertical temperature difference across the fluid layer could be regulated. The radial walls in the experiments needed to be insulating so that heat did not escape the experimental test fluid into the air surrounding the experiment. Lastly, some means of flow visualization needed to be employed in order to determine when the fluid layer was or was not convecting, so that the critical conditions for the onset of convection could be determined. Below, the means by which these goals were accomplished are described. 6.2 Experimental Apparatus As mentioned, the test fluid in the experiments was Dow Coming 200 1 Stoke silicone oil (thermophysical properties given in Appendix A). This fluid was chosen as the test fluid for this research because, for systems in the size range of those examined in this research, its thermophysical properties are such that the critical vertical temperature differences for convection were easily obtainable in experiments. The only thermophysical properties of the silicone oil which vary significantly with temperature are the density (enough to facilitate Rayleigh convection) and viscosity. In the experimental apparatus, a lucite ring acted as the outer radial boundary of the silicone oil test section. A lucite middlepiece was added when the test section was annular. The dimensions of these lucite pieces, thus, set the sizes of the cylindrical and annular test sections. A copper plate at the top of a continuously stirred water bath was the bottom boundary of the test section, and a sapphire window at the bottom of a flow through water bath was the top plate of the test section. To create heatedfrombelow conditions, the temperature of the bottom water bath was adjusted to certain setpoint values while the top water bath was always kept at a constant, cooler temperature. A process control computer, running a Lab VIEWTM process control program written by the author, sent signals to turn a heater in the experiment on and off as needed in order to keep the bottom water bath at desired setpoint temperatures. The process control system built for the experiments in this research will be further discussed in Section 6.3. The test section was insulated to prevent heat loss in the radial direction. To allow flow visualization, a small amount of aluminum tracer powder was mixed into the silicone oil. So that the flow behavior of the experimental system could be recorded and reviewed, a digital camcorder was mounted above the test section. More detailed descriptions of the components and features of the experimental apparatus, as well as some notes on the experimental startup procedure, are given in the following subsections. A simple diagram of the apparatus is shown below as Figure 61, and a couple of photographs of the apparatus follow as Figures 62 and 63. 6.2.1 Test Section As mentioned, a lucite ring acted as the outer radial boundary of the test section which contained the silicone oil. The bottom boundary of the test section was a copper plate, which was the top surface of a stirred water bath that will be described shortly. The top boundary of the test section was a sapphire window, which was the bottom surface of a flowthrough water bath, which also will be described. The lucite rings typically ranged from about 15 mm to 30 mm in diameter and about 6 mm to 8 mm in height. The sets of experiments for cylindrical systems were run before the sets for annular systems; this procedure was used so that when the annular experiments were run, a nesting hole in the bottom copper plate could be used to anchor the bottoms of the lucite centerpieces, which were added to transform the cylindrical systems to annular systems. A reason for selecting lucite as the material from which to construct the outer boundary rings and the centerpieces for the annular systems is its thermal conductivity, which is very close to that of the silicone oil (see Appendix A). Since the system of silicone oil being examined in this research was subjected to vertical temperature gradients, so, too, were the radial walls of the system. Thus, the fact that lucite has a thermal conductivity close to that of silicone oil made it an appealing choice as the material from which to construct the radial walls because it ensured that the vertical temperature gradients at the radial boundaries of the test section would not differ from those in the interior. This, and the fact that the lucite outer radial wall was surrounded by an insulating coating (which shall be further explained later), ensured that heat would not flow into or out of the silicone oil in the test section in the radial direction. After placing the lucite ring onto the copper plate and before filling the test section with oil, a lucite clamppiece was screwed down on top of the lucite test section ring to tighten the ring down against the copper plate (the clamppiece screws into the bottom water bath). This was an important measure taken to prevent any leakage of air/oil at the bottom of the test section. Once the test section ring was tightened down to the copper plate, the test section was intentionally overfilled with silicone oil. This was done so that the sapphire plate of the flow through top water bath could more easily be pressed down onto the top of the test section without 