THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING STABILITY OF THE UNSTRATIFIED AND STRATIFIED FLOW BETWEEN TWO HORIZONTAL, INFINITE DISKS ROTATING AT SLIGHTLY DIFFERENT ROTATION RATES At len R.: Brunsvo l'd A dissertation Submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan Department of Mechanical Engineering 1972 July, 1972 IP-845

t" ) i

To My Wife ii

ACKNOWLEDGMENTS I wish to express my gratitude to Professor Charles M. Vest, Chairman of the doctoral committee, for his advice, encouragement, and generous devotion of time throughout the course of this investigation. I wish also to express my appreciation to Professor Arpaci for serving on the thesis committee, for his personal interest and, above all, for his excellence as a teacher and scholar. Professors Debler and Graebel, who served on the thesis committee, were most helpful in discussing the goals of the research and in realistically appraising the effort. Generous financial support from the Department of Mechanical Engineering as a teaching fellow and as a research assistant is gratefully acknowledged. Additional support in the form of fellowships from the Ole Evinrude Foundation, Cummins Engine, the Ford Motor Company, and the Department of H.E.W through NDEA, as well as a traineeship from the National Science Foundation is also acknowledged. The computer support of The University of Michigan Computer Center through the Mechanical Engineering Department and the assistance of the Industry Program of the College of Engineering is gratefully acknowledged. Mrs. Gerry Nowak typed the final manuscript and her efficiency and advice is acknowledged. I wish to thank my wife and children for their sacrifices, encouragement and patience. Finally, I wish to thank my parents for their encouragement through the years and for their dedication to their family. iii

TABLE OF CONTENTS ACKNOWLEDGMENTS......................................... iii LIST OF FIGURES.................................-* —.... v LIST OF TABLES......................................... vii NOMENCLATURE..................v............ viii INTRODUCTION............................................ 1 CHAPTER 1. LITERATURE SURVEY........................... 4 1.1 Base Flow 4 1.2 Stability 8 1.2.1 Unstratified Flows 9 1. 2.2 Stratified Flows 10 CHAPTER 2. THEORETICAL ANALYSIS........................ 14 2.1 Basic Equations 14 2.2 Governing Equations of the Unstratified Flow 17 2.2.1 Unstratified Base Flow 17 2.2.2 Validity of the Unstratified Base Flow Solution 21 2.2.3 Stability Equations of the Unstratified Flow 25 2.2.4 Limiting Cases of the Unstratified Flow 31 2.3 Governing Equations of the Stratified Flow 33 2.3.1 Stratified Base Flow 33 2.3.2 Stability Equations of the Stratified Flow 38 2.3.3 Limiting Cases of the Stratified Flow 42 CHAPTER 3. NUMERICAL ANALYSIS —SOLUTION TO THE EIGENVALUE PROBLEM......................... 45 3.1 Finite Difference Representation 46 3.2 Algebraic Eigenvalue and Determinant Solutions 47 3.3 Convergence and Numerical Procedure 49 CHAPTER 4. RESULTS AND CONCLUSIONS..................... 63 4.1 Unstratified Flow 63 4.2 Stratified Flow 74 4.2.1 Critical Curves 78 4.2.2 Neutral Curves and Eigenvalue Spectra 95 APPENDIX I. STABILITY EQUATIONS OF THE STRATIFIED FLOW IN THE LOCAL COORDINATE SYSTEM............ 103 APPENDIX II. FINITE DIFFERENCE SOLUTION TO THE EIGENVALUE PROBLEM OF THE STRATIFIED FLOW...... 110 BIBLIOGRAPHY.....................1....................... 118 iv

LIST OF FIGURES Figure Page 1 Horizontal, Infinite Disks - Rotating at Slightly Different Rotation Rates........... 15 2 Comparison of Linearized Base Flow Solution, Eq. 12, to Nonlinear Solution of Lance and Rogers (1962), Q2/Q1 = 2.................... 22 3 Base Flow Velocity Profiles, Eq. 12, at Small T = QL2/u................................... 23 4 Schematic Showing Local Coordinate System...... 27 5 Convergence, Unstratified Flow - Growth Rate Ar, Near Critical Point as Function of Orientation, Y, and No. of Intervals, N, (T = 100) 62 6 Convergence, Unstratified Flow - Growth Rate, Xr, Near Critical Point as a Function of Wave No., a, and No. of Intervals, N, (T = 100)................................... 62 7 Critical Curve, Unstratified Flow - R it as Function of T..oo...o... oo............ 65 8 Wave No., a, and Orientation, y, at R Crit as Function of T......................... 66 9 Vestige of the "A" Wave....................... 68 10 Growth Rates as Functions of Reynold's No. in Vicinity of "A" and "B" Wave Intersection.......................... 71 11 Critical Grashof No., Gr as Function of V112 = crit R = RT-1/2, (Pr = 6.0).......... 79 12 Wave No., a, (at Gr ) as Function of = RT,=crit R RT-1/2 (Pr = 6.0)..................... 80 ~~~~~~

LIST OF FIGURES (cont.) Figure Page 13 Orientation, yO, (at Gr ) as Function -l'2 crit of R = RT-1, (Pr = 6.0)...... 81 14 Critical Grashof No., Grcrit, as Function of R = RT-, (Pr = 1.0) 82 15 Wave No., a, (at Grcrit) as Function of lccrit R = RT, (Pr = 1.0)........... 83 16 Orientation, y, (at Grcrit) as Function of -1/2 crit R = RT-/, (Pr = 1.0).84 17 Critical Grashof No., Gr as Function of crit' R = RT /2, (Pr = 0.025).o***9ooo* e 90 18 Wave No., a, (at Gr it) as Function of -1/2 crit R = RT, (Pr = 0.025).oo*..o.....o 91 19 Orientation, y0, (at Gr crit) as Function of -1Y crit R = RT 1/2 (Pr = 0.025).................. 92 20 Wave Speed, c r/R, and Frequency,?, (at Gr - 1/2crit as Functions of R = RT, (Pr = 0.025, T = 50).. 93 21 The Stabilizing Effect of Prandtl No., Pr, (T = 50)...... 96 22 Neutral Curve, Gr vs yo, (Pr 1.0, T = 100, R = 10, a = 5.42)....................... 100 23 Neutral Curve, Gr vs a, (Pr = 1.0, T = 100, R = 10, y = 75.00)......................... 100 24 Neutral Curve, Gr vs y~, (Pr = 1.0, T = 100, R = 40, a = 5.56).......................... 101 25 Neutral Curve, Gr vs a, (Pr = 1.0, T = 100, R = 40, Y = 37.0~)........................ 101 26 Neutral Curve, Gr vs a, (Pr = 1.0, T = 20, R = 67.1, Y = -50.0~)...................... 102 27 Growth Rate as Function of Grashof No., (Pr =1.0, T = 100, R = 40, y = 80~, a = 5.56)....... 102 vi

LIST OF TABLES Table Page I Convergence, Benard Problem - Rayleigh No., IR, as a Function of Wave No., a, and No. of Intervals, N.............. 56 II Convergence, Rotating Benard Problem - Grashof No., Gr, as a Function of Wave No., a, and No. of Intervals, N [T = 11.18 (Ta = 500) and Pr = 5.0]........................... 56 III Convergence, Stratified Flow - Grashof No., Gr, as a Function of Wave No., a, Orientation, y, and No. of Intervals, N (Pr = 1.0, T = 100, R = 100)................................. 59 IV Convergence, Unstratified Flow - Critical RE as a Function of Intervals, N................... 59 V Critical Results, Unstratified Flow............... 64 VI Critical Results, Stratified Flow (Pr = 6.0)....... 75 VII Critical Results, Stratified Flow (Pr = 1.0)...... 76 VIII Critical Results, Stratified Flow (Pr = 0.025).... 77 Vii

NOMENCLATURE A base flow parameter defined in Eqs. (18), or coefficient matrix, Eq. (32) a dimensional wave number a thermal diffusivity a unit vector in local coordinate system, Appendix I B base flow parameter defined in Eqs. (18), or coefficient matrix, Eq. (32) b unit vector in local coordinate system, Appendix I C coefficient matrix, Eq. (32) c dimensionless wave speed referred to mean velocity at mid-plane c dimensionless wave speed referred to disk velocity, defined in Eqs. (21) D operator, d/dz or d/dS Fr Froude number, 2rO/g Gr Grashof number, aATgL /v g gravitational acceleration h scale factor in local curvalinear coordinate system defined in Eqs. (13) I unit matrix ii unit vector, x1 coordinate i imaginary number, /-1 j 1 unit vector, y1 coordinate unit vector, z coordinate L gap between disks N number of difference intervals P pressure viii

NOMENCLATURE (cont.) Pr Prandtl number, v/a R Reynold's number,' UOL/v R modified Reynold's number, R//T ffi Rayleigh number, GrPr r radius rO radius at local coordinate origin ds differential length in local coordinate system, Eq. (13) T Taylor number, QL2 /v Ta traditional Taylor number, 4T2 T1_2 temperature of bottom or top disk, respectively AT temperature difference, T2 - T! t time U0 reference velocity, sr0 u dimensionless velocity in x direction, Eqs. (15) u dimensionless base flow similarity velocity in x direction, Eqs. (17), or in zonal direction in cylindrical coordinates u1 dimensionless base flow similarity velocity in X1 direction, Eqs. (18) V velocity vector V radial velocity V zonal velocity Vz vertical velocity v dimensionless velocity in y direction, Eqs. (15) v dimensionless base flow similarity velocity in y direction, Eqs. (17), or in radial direction in cylindrical coordinates v1 dimensionless base flow similarity velocity in Yl direction, Eqs. (18) ix

NOMENCLATURE (cont.) w dimensionless velocity in z direction, Eqs. (15) X coordinate, Fig. 4 x local skewed coordinate, Fig. 4 x1 local zonal coordinate, Fig. 4 Y coordinate, Fig. 4 y local skewed coordinate, Fig. 4 Y1 inward directed local radial coordinate, Fig. 4 z vertical coordinate Z vertical coordinate Z complex. velocity in cylindrical coordinates, V /r + i Vr/r (Eq. 12) Z1 complex velocity in local x1, Y1 coordinate system, Eq. (18) Greek Letters a dimensionless wave number a modified dimensionless wave number, a//T ~a coefficient of thermal expansion Y local coordinate skew angle in radians but referred to in degrees, Fig. 4 Y' perturbation parameter A length of difference interval, Eqs. (31) 2 2 2 A operator, D a 6 centered difference operator, Eqs. (31) E rotation parameter, (Q2 - Q1]/2 e1,2 dimensionless rotation parameter at bottom or top disk, respectively, Eqs. (5) e dimensionless temperature,(T - T1iAT x

NOMENCLATURE (cont.) A complex eigenvalue, Eqs. (32) v kinematic viscosity stretched coordinate at disk, I(2 + z) p mass density, Eq. (1) po reference density at T1 a dimensional complex phase velocity perturbation stream function, Eqs. (I-7) Q rotation at mid-plane Q1,2 rotation of bottom or top disk, respectively xi

INTRODUCTION Stability criteria of laminar viscous flows have represented a large portion of the studies in theoretical fluid dynamics for many years. Because of the complexity of the transport equations which govern the physical systems, the major portion of these studies have been limited to the stability criteria of steady, two-dimensional laminar viscous flows (hereafter referred to as base flows). Finding a two-dimensional representation of a three-dimensional phenomena is often a significant contribution in a stability study. In some cases, this may be justified on a mathematical basis, Squire (1933). In others, the justification is based on experimentally observed phenomena. This will be the case in the present study. The approach taken in stability studies is to consider the effect of a small time-dependent perturbation on the base flow. In the resulting suitably nondimensionalized perturbation equations, if only the lowest order terms are considered, the study is referred to as the linear stability theory and, if higher order terms are considered, the finite amplitude or nonlinear theory. The present study is concerned with the linear theory. The disturbances (small perturbations) are usually represented in normal modes as 4(z) e i(a.xy)+at a and 8 are real and are the wave numbers in the x and y directions. a = a + ia is the complex phase velocity. This reduces the governing perturbation equations, which were a set of homogeneous partial differential 1

2 equations with homogeneous boundary conditions, to a set of homogeneous ordinary differential equations with homogeneous boundary conditions. This set is further reduced, by various methods, to a set of linear algebraic equations which then constitute an algebraic eigenvalue problem. a is treated as the eigenvalue. When a is positive, the disturbance grows and the flow is unstable. Either a new laminar flow regime then occurs or turbulence sets in. When o is negative, the flow is stable. To determine when a flow is unstable, one seeks the locus in the nondimensional parameter space which delineates the stable from the unstable modes. This is frequently referred to as neutral or marginal stability and occurs when ar - 0. To complete the stability study, one must find the minimum of the a E 0 loci in the space of the parameters which define the flow. When i = 0, the disturbance is stationary; when a # 0, nonstationary. Experimentally, when instability occurs, the following types of phenomena have been observed. Stationary: i) Two- and three-dimensional cellular patterns. Nonstationary: i) Two- and three-dimensional oscillating cellular patterns. This phenomena is commonly referred to as overstability. ii) Two- and three-dimensional (skewed) traveling waves. iii) A stationary or traveling vortex or cell which may or may not be skewed.

3 It is notable that the skewed traveling waves and the skewed traveling vortex, or cell, may, in some cases, be accurately represented as two-dimensional phenomena. In this dissertation, the stability criteria of the motion of both an unstratified and a stratified viscous layer of fluid contained between two horizontally infinite coaxial disks rotating at slightly different rotational rates is studied. The mean rate of rotation, Q, and the vertical length scale, L, considered are such that the Taylor number, T = L 2/v, falls in the low-to-moderate range (O < T < 100). The major causes of instability considered are: shear, Coriolis acceleration, and the buoyancy forces due to stratification. The three-dimensional base flow of the problem is modeled by a similarity solution which is well established for the unstratified flow and which is extended to include the parameter range considered in the stratified flow. Past experimental studies of problems governed by the combined effects of shear and stratification and shear and Coriolis accelerations have shown that, in both, instability results in a new flow field which is banded in structure. Concluding that the shear of the present problem will also result in a banded structure and with the aid of the base flow similarity solution, a simplified formulation of the present problem is established. A local viewpoint is taken (similar to the approach taken in boundary layers) which assumes banded disturbances. Both stationary and nonstationary instabilities are predicted. The results obtained are shown to be consistent with all the known limiting cases of the flow.

CHAPTER 1 LITERATURE SURVEY 1.1 BASE FLOW Batchelor (1951) first suggested that the steady laminar flow of an unstratified viscous fluid between two infinite, coaxial disks could be treated using the similarity approach of von-Karman (1921) for the flow resulting from a single or free disk. The similarity equations result by assuming the base flow is axisymetric and that the axial velocity, V, is a function only of the axial coordinate, z. The result is that the radial velocity, Vr, and the zonal velocity, V0, are such that Vr/r and V0/r are functions only of z. There results a set of nonlinear, ordinary differential equations governing the base flow. From the general nature of this set of equations (hereafter referred to as the similarity equations), Batchelor discussed the permissible flow regimes for various boundary conditions. Stewartson (1953) pursued some of the possible solutions to the similarity equations and furnished some preliminary experimental evidence of the validity of the approach. In his paper Stewartson also demonstrated that, for small shear (disks rotating in the same sense with slightly different rotation rates), the set of nonlinear differential equations is reduced to a linear set. Lance and Rogers (1962) gave numerical solutions to the similarity equations for various boundary conditions. Pearson (1965) numerically integrated the time-dependent similarity equations for 4

5 a number of boundary conditions imposed on an initially quiescent fluid. One or both of the disks' rotation rates was impulsively changed and the calculations were continued until a steady state was attained. Pearson's results are in general agreement with the results of Lance and Rogers. Mellor, Chapple, and Stokes (1968) pursued both numerical and asymptotic solutions for the case of one stationary disk and presented experimental results which verified the analytical results of the simplest flow regime. Rogers and Lance (1964), in a numerical study of the boundary layer on a disk of finite radius in a rotating fluid, found that for r/r0 < 1/2, where r0 is the outer radius, their results agreed with the numerical results of the similarity solution for the infinite disk given by Bodewadt (1940). The motion between finite, concentric disks has also been studied experimentally by SchultzGrunow (1935), Stewartson (1953), and Picha and Eckert (1958). The conclusion of these experimental studies is that for shrouded disks (when the fluid is reingested to the system), the resulting flow regimes are in agreement with the predictions of the similarity solutions of the infinite case. This is not always the case for unshrouded disks, especially when the disks rotate in opposite directions. The experimental results of Mellor, Chapple, and Stokes (1968) for one disk stationary show that, at large R, the unshrouded results are in good agreement with the similarity solution. The hot wire probes, however, indicate that even with the shrouding removed, the fluid was being reingested to the system.

6 The combined effect of rotation and stratification (without shear) is considered next. Zeipel (1924) showed that solid body rotation of a stratified body of fluid is not possible. From the viewpoint of a rotating reference, the fluid is driven horizontally both parallel to and perpendicular to any horizontal pressure gradient. When the temperature differences are small, the ratio of the inertial to the Coriolis accelerations (the thermal Rossby number) is small and the perpendicular component of the motion dominates. This phenomenon is commonly referred to as geostrophic flow and the zonal motion which occurs in the interior of the fluid is frequently referred to as the thermal wind. Thus, convection must occur in any rotating, stratified layer of fluid. However, as discussed by Greenspan (1969) in his monograph (pp. 13-14), there are many physical situations for which the time scales of the convective processes are much greater than the time scales of the rotation. The literature regarding the importance of this motion, considered pertinent to the present discussion, can conveniently be separated into two categories by considering the Taylor number, T = OL /v. The first is the low-to-moderate Taylor number range for which it is assumed the time scales of the convective processes are much greater than the rotational time scales. Chandrasekhar (1953), studying the instability of a rotating, stratified layer of fluid (hereafter referred to as the rotating Benard problem), assumed the base flow is one of solid body rotation. This implies that the "excess centrifugal force" caused by stratification is negligible. The rotating Benard problem is also discussed

7 in detail in Chandrasekhar's monograph (1961) where he gives reference to subsequent experimental studies. These references will be discussed in more detail in the next section. At this point, it is sufficient to note that there is ample evidence that the above assumption is valid in a large range of the parameters of physical interest. This range will be discussed explicitly in the next section. The second category, that of large Taylor number, is considerably more complicated. Although the present study is limited to the low-to-moderate Taylor number range, some of the large Taylor number results are relevant to the task of establishing a valid representation of the stratified base flow, as well as its limitations. It should be noted that, in this range of Taylor number, the parameter usually defined is the Ekman number, which is proportional to the inverse of the Taylor number, as defined herein. Millsaps and Pohlhausen (1952), considering the heat transfer from a single rotating disk in an incompressible fluid, demonstrated that the full Navier-Stokes equations are reduced to the similarity equations with the energy equation uncoupled from the momentum equations. Riley (1964) considered the heat transfer from a single rotating disk in a compressible fluid. After noting that, on a single disk the pressure is constant, assuming the viscosity is proportional to the temperature and making the von-Mises transformation, he demonstrated that the momentum equations are uncoupled from the energy equation. On the basis of this, Riley assumed the similarity equations are valid.

8 A number of authors have considered radially bounded convective motions for the large Taylor number range. These include Duncan (1966), Barcilon and Pedlosky (1967a, 1967b, 1967c) and Homsey and Hudson (1969). The important observation (relevant to the present study) is that even at large Taylor numbers, the similarity equations are often valid. However, this is not always true as the effect of stratification and the effect of the side walls can significantly change the nature of the flow. In general, the resulting flow fields are described as boundary layers (referred to as Ekman layers) near the disks, geostrophic interior flows and a rather complicated system of boundary layers near the side walls, which depend on the boundary conditions as well as the parameters of the governing equations. In the case of stable stratification, the Ekman layers can be suppressed if the stratification is sufficient. A very important conclusion from the above references is that, at large Taylor numbers, the effect of the side walls must be considered in order to represent real physical situations. This is a result of the boundary layer nature of the flow and is, therefore, not necessarily true at small-to-moderate Taylor numbers. No literature has been found which includes rotation, stratification in the same direction as the rotation vector and independently controlled shear as is the case in the present study. 1.2 STABILITY For convenience, the literature concerning unstratified and stratified phenomena are separated.

9 1.2.1 Unstratified Flows Since shear is considered in the present study, the nature of instability in a sheared layer of fluid which does not rotate is of interest. C. C. Lin (1966) in his monograph, The Theory of Hydrodynamic Stability, provides a comprehensive treatment on the stability of parallel shear flows. Of particular interest here is the stability of a sheared layer of fluid between two horizontal parallel plates. This problem is referred to as plane Couette flow (pCf hereafter). Hopf (1914), Wasow (1953), Zondek and Thomas (1953), and Grone (1954) have studied pCf by asymptotic methods as oR -+ o. Gallagher and Mercer (1962) studied pCf for finite Reynolds numbers using the Galerkin method. Deardorff (1963) studied pCf at finite Reynolds numbers using finite differences. The conclusions from the above studies are: as aR -+ (and considering the least stable disturbance), the linear theory predicts pCf to be stable with the wave speed approaching +1 and the wave number approaching 1. For finite R, the linear theory predicts pCf stable to R = 1430 with no indication that instability should be expected at larger R. When considering the effect of rotation, there are two problems of particular interest. These are, the flow induced by a rotating disk in a semi-infinite body of stationary fluid and, the socalled Ekman layer problem. One configuration which produces an Ekman layer is a rotating horizontal tank of fluid with sources distributed around the periphery and a sink at the center. Gregory, Stuart, and Walker (1955) (hereafter referred to as GSW) presented experimental and analytical results for the rotating disk. Experimentally, Gregory and Walker, in GSW, observed

10 instability as stationary roll vortices (stationary with respect to the disk) oriented at a small angle to the zonal direction on the disk. From this observation, Stuart, in GSW, used von-Karman's similarity solution to the free disk flow as a base flow and considered the stability of this base flow by using a local formulation which assumed perturbations as roll vortices. The analysis considered only the inviscid case, but led to successful prediction of certain features of the resulting unstable flow. Faller (1963) experimentally investigated the Ekman boundary layer flow. Two types of instability were observed. One was similar to that found by GSW. The other was a traveling wave oriented differently than the stationary wave. It should be noted that the method of observation (china clay) used by GSW precluded the observation of traveling waves. Lilly (1966) and Faller and Kaylor (1966) investigated the Ekman layer analytically using a local formulation, following Stuart, but included the effect of viscosity. The results were in excellent agreement with Faller's experimental results. 1.2.2 Stratified Flows A number of authors have investigated the Benard problem (a layer of fluid, heated from below) both experimentally and analytically. The Benard problem is comprehensively discussed by Chandrasekhar (1961) in his monograph. Chandra (1938), in an experimental study of plane Couette flow heated from below (denoted stratified pCf herein), observed that shear has no effect on the critical Rayleigh number, R.

11 That is, the W is 1708, the same as for the Benard problem. c However, Chandra did observe a tendency for instability to result in longitudinal roll vortices (axis of roll in direction of shear). This tendency was especially pronounced at large R (i.e., large shear). Brunt (1951) observed the same tendency although, at very small rates of shear, he did notice transverse rolls. Deardorff (1965), in a numerical study of stratified pCf, notes that the transverse rolls observed by Brunt were probably caused by finite amplitude phenomena at supercritical IR. Gallagher and Mercer (1965) applied the Galerkin method in a numerical study of stratified pCf. Using Squire's theorem, they study a twodimensional transverse wave and demonstrate that the solution indicates that the longitudinal waves are the most critical. Chandrasekhar (1953) and Chandrasekhar and Elbert (1955) studied the rotating Benard problem. The approach was to assume solid body rotation as the base flow by neglecting the "excess centrifugal force" term which is also neglected in the perturbation equations. These results, as well as the results of subsequent experimental studies, (Fultz, Nakagawa and Frenzen, 1954, 1955, Nakagawa and Frenzen, 1955, Dropkin and Globe, 1959, Goroff, 1960) are comprehensively discussed in Chandrasekhar (1961). The experimental works, cited above, are in excellent agreement with the theoretical predictions. Veronis (1959, 1966) applied the nonlinear theory to the rotating Benard problem and concluded that finite amplitude disturbances could be sustained at subcritical IR (iR lower than that

12 predicted by the linear theory) for Prandtl number, Pr < /2 and moderate T. Rossby (1969), in an experimental study, observed subcritical instability for the rotating Benard convection in mercury (Pr = 0.025) in the range 67 < T < 156 and in water for T > 112. Rossby attributed the subcritical results in mercury to finite amplitude phenomena, in partial agreement with Veronis. At present, the cause of the subcritical instability observed by Rossby in water is not known. It should be noted that they occurred at greater T than will be considered in the present study. Niiler and Bisshopp (1965) considered the horizontally unbounded rotating Benard problem as T + A. The approach was to consider the asymptotic form of the characteristic equations for the Rayleigh number which results from the normal mode analysis of Chandrasekhar (1953). They concluded an Ekman layer approach is appropriate with the main result that the critical ~R becomes independent of the boundary conditions, i.e., the results for rigid boundaries approaches the results for free boundaries as T -+. Homsey and Hudson (1971a), by asymptotic or boundary layer methods, considered the stability of a stratified, rotating, quiescent layer of fluid as T + X which included the effects of side walls. The quiescent base flow is a result of neglecting the "excess centrifugal force" term, as discussed above. They concluded that when centrifugal effects are ignored, the effect of side walls is negligible for all but very tall cylinders of fluid. Homsey and Hudson (1971b) also considered the case which includes the "excess centrifugal force" term. They concluded that away from the side wall, the effect of the base flow circulation is stabilizing but that near the wall a destabilization is possible.

13 References which include the effects of rotation, stratification, and independently controlled shear in a horizontal layer of fluid, as in the present study, have not been found.

CHAPTER 2 THEORETICAL ANALYSIS 2.1 BASIC EQUATIONS Consider the flow of a Newtonian fluid between two infinite, horizontal, coaxial disks rotating at slightly different rotation rates (discussed in paragraph 2.2.2), as shown below (Fig. 1). Both the stratified and the unstratified flow regimes will be considered. The invariant form of the basic equations in a rotating reference are given below for the stratified case. The following assumptions concerning the fluid are made: 1) Viscous dissipation and the work of compression are negligible. 2) The kinematic viscosity, v, and the thermal diffusivity, a, are constant. 3) The density, p, is a linear function of the temperature, e = T - T1, and is given by: C='c (1- is) (1) 4) The density variation is important only when considering the gravity and centrifugal forces. It will be shown, in the stratified case, that the "excess centrifugal force'" term can be neglected for the parameter range of interest. 14

15 T2 5 t~; 0r L^ T2. - I1 7...., i Figure 1. Horizontal, Infinite Disks - Rotating at Slightly Different Rotation Rates.

16 This leads to the following set of equations, written in a rotating coordinate frame with origin at mid-plane: (2) DV Dt l k AV v ve AC - + 0 ( 9 where: P - t o - 1 P0 is the hydrostatic pressure, and 2 is the rotation rate of the fluid at the midplane and must be determined. The basic equations for the unstratified flow regime result by setting 0 O0 in Eqs. (2). The stability of the system is investigated by considering the effect of a small perturbation on the base flow. Analytically, this is accomplished by substituting variables of the form (Kit): vA + X \V ("It) into the basic equations and collecting terms of like order in y', where y' is a small perturbation parameter. The lowest order, (y'), equations are the base flow equations. The (y')l order equations are the linearized perturbation or disturbance equations.

17 In modeling the unstratified and the stratified flow regimes, it is necessary to consider different restrictions or assumptions in each case. For this reason, the details of the development, outlined above, are considered separately. 2.2 GOVERNING EQUATIONS OF THE UNSTRATIFIED FLOW 2.2.1 Unstratified Base Flow The governing equations for the base flow in the unstratified case are developed in a rotating cylindrical coordinate system. The flow is assumed to be axisymmetric and Batchelor's similarity assumption is made. The restriction that the disks are rotating at slightly different rotation rates (i.e., small shear discussed in paragraph 2.2.1) is used to eliminate the nonlinear terms in the ordinary differential equations which result from the similarity assumption. The governing equations of the base flow are nondimensionalized by choosing, where dimensional variables are now denoted by ( ): C t} ) g - t rO t' ) L ~ i\ z Et v) )v, II,v ~ JVVqV (3) p- U0 p where U0 is chosen as the shear velocity given by Tuh= (L ro g ern (4) The nondimensionalized governing equations of the axisymmetric base flow are'

18 XV/+ vrv* +.V& -} (( V4 + () which must satisfy v _L _ \/;/r = t =&:l) \ = ) T ~0 SLt; _ -- V/r: =: = ~ V 0 - V~ 0 f~3 where j r rt 2 and T =, is called the Taylor number herein. Following Batchelor, it is assumed that Vz V (z). _ti,-iz z Equations(5) then require V /r = v(z) and V /r = ut(z). ro thefrsE From the first of Eqs. (5) =-L d.=- T I(

19 where W(z) is some unknown function of z and the T-1/2 is included by anticipating the solution. Equations 5 become, with primes denoting, d/dz, V V -u-T WZV — +T V tVi% &UT~~~ W l~~jZ~U —~ + V ~(7) t uV-T +W- T UL T(Lt- | viW + W = - -3 From the first and last of Eqs. (7), P - r FR (, where P1 is a constant. Eliminating the pressure and noting that, for small shear, /Q2 << 1, the first two of Eqs. (7) become V T+ ZTU -D u -2TV O E 1 DP Note that, the term -1 r ar was not dropped —it was mathematically eliminated. The reason for this is that the order of P (through P1) with respect to the velocities may not be obvious. It is shown below that P1i 0. Define Z = uL + *V (8) The differential equation to be solved is then

20 -l —) T Z - D (9) which must satisfy the boundary conditions; X(t) - CE ) Z- 1) I61 and from the first two of Equations 7; Si ( a L ZT(,, e ) The solution to Eq. (9) which satisfies the boundary conditions is rJ -T (I L) f -t jC,-,Tminlf ( I ia- I/r (10) where C 1 Q 1 From the first of Eqs. (5) and the boundary conditions on.V, the solution must satisfy, 3m Zdz + constant = 0 at z =+ From the requirement of no net radial flow'IL - 12The last two conditions result in P =- 21 + f2

21 The final solution is Z *(J) (12) Note that, since no restriction was made on the aspect ratio, L/ro, Eq. (12) is valid for all radii, and any limit on L/r0 would be a result of the restriction discussed in the next paragraph which concludes T = < 100. The above solution would also result by nondimensionalizing Stewartson's solution for small shear and transforming it to the present rotating reference with origin at the mid-plane. 2.2.2 Validity of the Unstratified Base Flow Solution In considering the validity of the base flow representation, two questions must be given attention. First, is the similarity assumption valid? Second, to what extent is the small shear solution, Eq. (12), valid? As discussed in paragraph 1.1, considerable experimental evidence is available which is in excellent agreement with the solutions to the similarity equations. In an effort to determine the range of validity of Eq. (12), numerical solutions were compared to available solutions of the nonlinear similarity equations (c.f., Eqs. (7)) provided by Lance and Rogers (1968). The results of Lance and Rogers for two corotatiag disks at a shear rate of e/Q = 1/3 and various T are compared to numerical solutions of Eq. (12) in Fig. 2. Even at this large shear rate, the linearized solution is a very good approximation

22 0.5- x 0.4 xx 0. 3 o z 0 x 0.2 x/ T = 18.7 0.1 /0 0 0.0,, I,, I. 0. 5._ 0.4 o 0.3 0 z0~~~~ O / tT = 48 0.21 o 0 o \ {/ 0.0 0.5 _ 0.4 O <0 < T = 126 0.3 0 0.2 e r/r 0.1 x 0x V0/r Lance & Rogers 0.1 e x O.l 0 nx x -Vr/r (E/I: 1/3) 0.0,. -0.2 0.0 0.2 0.4 0.6 0.8 1.0 Figure 2. Comparison of Linearized Base Flow Solution, Eq. 12, to Non-Linear Solution of Lance and Rogers (1962), /F:i = 2.

23 0.5 0.4 4 \ 0.3 0.2 ~~~~I T=I 0.5 - 0.4 0.3 I / T -5 0.2 O.1 I / 0.0,,'' i...', 0.5. 0.4 0.3 i )T = 10 0.2 _ 0.i -- r/r 0.0 I I. 0.0 0.2 0.4 0.6 0.8 1.0 Figure 3. Base Flow Velocity Profiles, Eq. 12, at Small T _QL U

24 to the nonlinear solution for this case. The only significant difference is the rotation rate of the fluid at the mid-plane. Notably, the shape of the velocity profiles are in excellent agreement at the lower T and in very good agreement even at T = 126. On the basis of this, it might be expected that Eq. (12) is accurate to shear rates as large as O(1)The boundary layer nature of the solutions is quite evident in Fig. 2 at T = 126. Since a boundary layer flow can be expected at large T, it follows that the parameter c/Q would no longer be a good measure of the intense shear occurring in the boundary layers. This leads to the conclusion that the parametric region of validity of Eq. (7) also includes an upper limit on the parameter T. From Fig. 2, one might expect a boundary layer analysis to be necessary for T > 100. However, as will be demonstrated in the next section, this will not create any difficulties in the unstratified case since the limiting form of Eq. (12) as T + oo is a boundary layer solution, namely, the Ekman layer solution. Fortunately, the stability of the Ekman layer has already been studied, both experimentally by Faller (1963) and analytically by Lilly (1966) and Faller and Kaylor (1966). This does, in fact, provide an excellent check for part of the present work. For later reference, the velocity given by Eq. (12) is shown in Fig. 3 for small T. In conclusion, the unstratified base flow solution, Eq. (12), is found to be valid for 0 < T < 100 and at all radii (c.f. paragraph 2.2.1).

25 2.2.3 Stability Equations of the Unstratified Flow The stability analysis is performed in a local coordinate system which is located at some radius, ro, much greater than the depth of the fluid, L. The justification of this approach has its basis in the following. The limiting form of the base flow, Eq. (12), as T + 0 is the plane Couette flow (pCf) velocity profile. Although pCf is suspected to be stable to all infinitesimal disturbances, the type of instabilities one would expect, should the flow be unstable, are those associated with parallel shear flows which are two-dimensional traveling waves. The formulation used below is quite similar to that used in studying the stability of parallel flows and does satisfy this limit. In the limit T + o, the base flow becomes the Ekman layer. Two experimental studies have been made which are of interest. GSW studied the stability of the flow resulting from a free-disk rotating in still air. They observed that instability occurred in the boundary layer as skewed, stationary vortex rings (skewed with respect to a zonal velocity and stationary with respect to the disk). Furthermore, and equally important, these instabilities occurred at radii much greater than the boundary layer depth. On the basis of this Stuart, in GSW, in an inviscid analysis, successfully predicted certain features of the secondary flow. Stuart used a local boundary layer viewpoint at a large radius which assumed a banded disturbance structure.

26 The other experiment was performed by Faller (1963) on the Ekman layer flow. Faller observed two instabilities. One was an almost stationary wave similar to that observed by GSW and the other was a skewed traveling wave, both occurring at radii much greater than the Ekman layer depth. In both of the experiments discussed above, the waves, except for their three dimensionality, are similar to the waves observed in unstable parallel shear flows. Lilly (1966) and Faller and Kaylor (1966), using the approach of Stuart, studied the viscous Ekman layer stability. Their results are in excellent agreement with Faller's experimental results. Since the present problem is similar to those discussed above, a formulation similar to Stuart's is taken. In the present case, however, divergence of boundary layers need not be considered. The local curvilinear coordinate system of Fig. 4 is chosen with origin at ro, where r0 >> L0. In the (x, y) coordinate system, X t -'0L Z (13) A differential length is given by d4a= h~drL ~ dX dt where ( _ I K*~(~ Y) ro The perturbation equations are developed from Eqs. (2) in this local coordinate system. The details of this development are presented in Appendix I for the stratified case; the unstratified

27 y + Y -A LLy Q2 4a Figure 4. Schematic Showing Local Coordinate System.

28 equations follow by setting 0e 0. A short discussion of the development of Appendix I, as it would apply to the unstratified case, follows. After expressing Eqs. (2) in the correct form for the local coordinate system, a perturbation to the base flow is considered in the form A h ~ A A v( X9,3,)) - V + V, )e (14) The form of the perturbation assumed is a product of a global description and a local description of the perturbation variable. This approach was first taken by Stuart in GSW. The form of the global description is not known and it is shown in Appendix I that by virtue of the instability occurring at a large radius (i.e. rO >> L) it need not be anticipated. Locally, the disturbance is assumed to be two-dimensional. This is based on the physical observation, in the experiments discussed above, that the disturbances are banded in structure. The elimination of the global description of the disturbance and the assumption of a twodimensional local disturbance is equivalent to neglecting the curvature of the band structure. It should be noted that in order to take a local viewpoint, which is valid for all radii of importance, the base flow must have a similarity representation, which it does. The resulting local perturbation equations are nondimensionalized by choosing

29 (15) [ "W- L,V, p = t doL L-~ L After dropping terms of O0rL), the equations which result are o 4V V 7, _; |( i/4 $-v~" _-TDU-O -T 0(16) 22 7,l-iXl( /-9)l' -Lk -i L - o -( -- \)<P - ~ Pt _ =+ /2 where R - -, a local Reynolds number T =, the Taylor number v D = d (perturbation quantities) dz 2 2 2 V D -2 d (base flow quantities) primes denote dz + is the local perturbation stream function, defined in Appendix I (Eqs. I-7)

30 The quantities u and v are the z-dependent similarity velocity profiles in the skewed x, y coordinate system and are given by L - cot' + V,4a (17) where ul and vl are the x1, yl components of the base flow obtained by transforming Eq. (12) from the cylindrical coordinate system to the local x1, Y1 coordinate system. Z becomes 7(18( 1 8(18) with components At + A - -nkcut4e -Z ~ A Z Equations (16) constitute an eigenvalue problem. Discussion of how the eigenvalue problem is treated numerically is included in a separate chapter.

31 Before continuing with the development of the governing equations for the stratified case, the limiting form of the governing equations for the unstratified case as T -* 0 and T -+ are developed. 2.2.4 Limiting Cases of the Unstratified Flow As T + 0, the base flow velocity profiles become (on applying L'Hopital's rule): thL= 2a (19) V,- 7 Setting T - 0 and with y = 900 (c.f., Eqs. (17)), the perturbation equations, Eqs. (16), reduce to the traditional Orr-Somerfeld equation for pCf. As discussed above, pCf is suspected to be stable to all infinitesimal disturbances. In the limit T +oo, both the governing equations of the base flow, Eqs. (7), and the disturbance, Eqs. (16), are singular. On transforming the base flow solutions, Eqs. (18), to the bottom disk, choosing Q1 as the reference rotation rate, taking the limit of the transformed equations as T + o and redefining the vertical space variable as the base flow velocity profiles become (20) which is recognized as the familiar Ekman spiral (c.f.,Ekman, 1905).

32 Considering the base flow continuity equation in the local coordinate frame (see Appendix I, Eq. (I-6)), after stretching the z coordinate, as above, the appropriate base flow velocity scaling appears as UO0Iv; thus define Similarly, from the disturbance continuity equation, define To obtain the correct form of Eqs. (16) as T -+ a, divide the first by T, the second by T, and expand ~ and u as L, = F - v T + The first-order inner expansion of Eqs. (16) becomes X -io,R(l V d; z6E(4 e\1(> 1 u,o= ~' ZU(21) V, Q0- LK(RIIY ytQo- 4 2 ~ ) with D) ()- -C)O At I- O where C E- - E A1The second boundary condition was obtained by observing that the limiting form of the base flow gives zero shear as E + oo The outer expansion is trivial, with the conclusion that the fluid at; + rotates in solid body rotation with zonal velocity rQ.

33 Equations (20) and (21) are identical to those used by Lilly (1966) in his study of the stability of the Ekman layer. For convenience, the present parameters are related to the Ekman parameters by Rj= by (22) The first relation in Eq. (22) is particularly useful in presenting the results. This completes the formulation for the unstratified flow regime. The formulation for the stratified flow regime follows. 2.3 GOVERNING EQUATIONS OFTHE STRATIFIED FLOW 2.3.1 Stratified Base Flow The nondimensional governing equations of the stratified base flow are written below in cylindrical coordinates. L.1 "rVr _V = C. (23) rKV2. =AC>P r ___L +, r i +am +T l\ STi) rA -'~~~ - - -- -- t - +

34 -_V V -L rl V b t b rt 8 with -_~ = Oa la in addition to the boundary conditions of Eqs. (5).....t as before. Whenever the "excess centrifugal force" term, ~/Q re, is negligible, the r and e momentum equations are coupled to the energy equation only through the pressure. Since the r and 0 momentum equations are now identical with those of the unstratified case, and since, in the unstratified case the r and 0 momentum equations are sufficient to determine the similarity velocity profiles, this suggests assuming the similarity velocity profiles and determining the appropriate form of the pressure and the temperature profile. This approach is taken and certain restrictions this places on the problem are discussed in a later paragraph. Neglecting the "excess centrifugal force" term and making the similarity assumption, viz, Vz = V (z), it is evident from Eqs. (23) that the pressure and temperature must take the form P r'P,,( (24) without regard to ordering.

35 The energy equation becomes, making use of Eq. (6), d 4uRPrrT Lv1- o or, from the parameter definitions - _Pr Y'0/ (25) it is immediately apparent that if the parameter L PrT << 1, the conductive temperature profile results. However, in the case of finite T, this condition is too restrictive and will be relaxed. Considering Eq. (25) and requiring that the gradient of -+ -1 as F PrT +1 0, the gradient becomes dO _ jfrfOT A solution to wdz must be of the form 1 W1( z)and, as can be seen from Eq. (6), must be odd since v is odd (c.f., Eq. (12)). Then, do - DPr,(FT3) ed $(26) and expressing the exponential in its series representation, 6 L I __ - I- [towI 4-LPv4A + -r WI — ] considering the parameter - Pr to be small, a two term expansion would be pa -C1-i _ ()l 1 (27) J3~~~~~

36 Thus, the first-order correction to the gradient is O ( P and the first-order correction to the temperature profile, at finite T, would be O(i Pr. On the basis of the above, only the conductive basic state is considered. The base flow temperature profile becomes @ 1, (1 Jr () 28) with the velocity profile, in cylindrical coordinates, given by Eq. (12). It should be noted that the temperature profile is valid at all radii since, on neglecting the "excess centrifugal force" term and assuming the similarity velocity profile, the temperature is not a function of the radius. The restriction the neglect of the "excess centrifugal force" term presents can be considered from two points of view. First, the "excess centrifugal force" is normally neglected if it is of smaller order than the "excess" gravitational force. Comparing these terms, OtAT/Ct 1 __ 2 0 where Fr = is the Froude number. As stated above, the traditional criterion given to justify the neglect of the "excess centrifugal force" term is that Fr << 1. However, this condition is too restrictive and, in the present problem, not entirely relevant.

37 An alternate point of view is to consider the "excess centrifugal force" parameter in the second of Eqs. (23) by itself. The AfT term called the "excess centrifugal force" parameter above is, in fact, a measure of the ratio of the centrifugal acceleration caused by stratification to the Coriolis acceleration. It is evident that whenever J AT K< < (29) the "excess centrifugal force" term can be neglected. The significance of Eq. (29) is that it is independent of the radial distance, ro, and therefore the "excess centrifugal force" term may be negligible even when the Froude number is greater than 0(1). Furthermore, consider the rotating Benard problem for which (even without shear, as discussed in paragraph 1.2.2) the results of the linearized stability theory which neglects the "excess centrifugal force" term are in excellent agreement with the experimental results. In some of these experiments, the Froude number was greater than one. Intuitively, it would seem that shear would only improve this situation since the shear (as imposed in the present problem) imparts a zonal motion. This motion, in turn, results in the Coriolis acceleration which is the acceleration with which the "excess centrifugal force" is compared to in Eq. (29).

38 The discussion of Greenspan concerning the time scales of convective and rotating phenomena (c.f., paragraph 1.2.2) is relevant to the above discussion. In this case, however, the relevant time scales are the time scales of motions due to centrifugal effects and to Coriolis effects. Since the ratio of their respective accelerations is considered to be small in the present problem, a similar result would be obtained for their time scales. The perturbation equations of the stratified flow regime are considered next. 2.3.2 Stability Equations of the Stratified Flow In formulating the governing equations for the stability of the stratified flow, the local viewpoint is again used. The discussion below concerns the motivation of this approach. The main questions raised regarding the local analysis are: 1) Is a local form of the disturbance valid in this case? 2) Can the curvilinear terms be eliminated in this case? Two interpretations of the stratified problem will be considered. First, the stratified problem is interpreted as a study of the effect of a small but increasing stratification, AT, on the unstratified flow regime. Consider the effect of a AT too small to cause thermal convection by itself. It is evident that, for such a AT, the discussion of paragraph 2.2.3 would be valid with the

39 conclusion that the local viewpoint would be valid in the immediate vicinity of the parameter space (R, T, ac, y, c) which causes the unstratified flow regime to go unstable. Similar reasoning can be applied to each of the parameters which represent a potential cause of instability. This is best summarized by first deriving the perturbation equations and then observing the limiting form of these equations as the relevant parameters become very small. This is done in the next paragraph. However, the reasoning above would seem sufficient motivation to pursue the local formulation. The other interpretation is to consider the effect of small shear, c/Q, on the rotating Benard problem. In both the Benard problem and the rotating Benard problem, the experiments result in cellular convection which indicates this is the most efficient means of the transfer of energy. This continues to be the case in the rotating Benard problem to Taylor number greater than the Taylor number for which the present similarity base flow solution is considered valid. There is some evidence that shear would not after this fact. Consider the effect of shear on the Benard problem or conversely, plane Couette flow heated from below (referred to as stratified pCf hereafter). The experiments of Chandra (1938) and Brunt (1951) on stratified pCf demonstrated that the only effect of shear is to orient the rolls with axes in the direction of shear. Extending the effect of shear on the stratified pCf to the rotating case, with regard to orientation only, one would expect the shear would impart a regular orientation for a given set of

40 parameters and, furthermore, limit the curvature of the rolls to the general direction of the inflowing or outflowing fluid. Rossby (1969) for T < 100 and Pr = 100 observed "buckled" circumferential rolls in the central region of the rotating Benard convection, randomly oriented "buckled" rolls at intermediate radii, and radially oriented rolls at larger radii, the radial orientation being caused by small centrifugal forces. Significantly, the small centrifugal force does not destroy the roll but serves to orient it. Koschmeider (1967), in an experimental study of rotating Benard convection in the range O < T < 16, found axisymmetric cellular convection which began at the outer rim. Homsey and Hudson (1971) explain this from results obtained in a boundary layer analysis as T + A. They conclude there are Ekman layers near the co-rotating disks which cause a shear layer and a radial temperature gradient near the outer rim, the latter causing the effect found by Koschmeider. The significant point is that the convective roll is not destroyed, but given a preferred, regular orientation. Based on the above discussion, it is concluded that, in the parameter range of interest, the effect of small shear on the rotating Benard problem will be to impart a preferred orientation to the convective rolls. This orientation would likely be determined by the directions of maximum shear which are, conveniently, anti-parallel in the rotating reference chosen. The wavelength would be expected to be O(L) and the curvature would not be expected

41 to be large. Side wall effects in a typical experiment would not be expected to be important because the base flow is limited to T < 100 where viscous effects permeate the vertical depth (c.f., Figs. 2 and 3). On the basis of the above discussion and in conjunction with the discussion on the base flow, the local formulation was used for the stratified flow regime limiting the parameter range to T < 100 and R < R for Gr = 0, the unstratified case. A separate restriction on the aspect ratio is not considered necessary since L is inherently restricted by considering the restriction 0 < T < 100. The governing equations for the stratified stability, in the local coordinate system, are developed in detail in Appendix I. A discussion of this development is not included here as it would only repeat the discussion of paragraph 2.2.3. The governing equations obtained are (30) \1, i t t 8 t(T- - = 0d - Z~~P-P

42 with boundary conditions = -u-G where u and v are again defined by Eqs. (17) and (18), and use has already been made of Eq. (28) for 0, the base flow temperature profile. Comparing Eqs. (30) with Eqs. (I-10) of Appendix I, the "excess centrifugal force" term has been dropped in Eqs. (30). This is consistent with the limitations on the stratified flow regime, as imposed by the base flow representation (c.f., paragraph 2.3.1). 2.3.3 Limiting Cases of the Stratified Flow Before the limiting cases of the governing equations of the stratified flow regime are pursued, it is convenient to summarize the limitations (either imposed or found) on the physical parameters, which are -- < 1 (i.e., small shear, c.f., paragraph 2.2.2) 0 < T < 100 (c.f., paragraph 2.2.2) ~aT <c 1 (c.f., paragraph 2.3.1) E/Q Pr -+ finite when T 0(1) or greater such that, C_ Premains small which establishes Q,T the conduction regime (c.f., paragraph 2.3.1) It should be noted that as either T + 0 or R + 0, the Pr restriction can be relaxed, although the above restriction is certainly not serious.

43 In addition to the above restrictions, as was noted in paragraph 2.3.2, the curvilinear terms in the local perturbation equations are dropped by considering L/r0 to be small. However, also discussed in paragraph 2.3.2, this restriction should not be extreme because of the strong tendency for convective rolls in thermal instability to be oriented in the direction of shear and the nature of the hydrodynamic instabilities in the Ekman layer, i.e., waves observed occurring at large radius. In other words, it is felt that the curvilinear terms can be neglected whenever 0 < T < 100. This conclusion is reinforced by observing the numerous correct limiting forms of Eqs. (30) given below: Pr = 0 or Gr = 0 With Pr = 0, the velocity and temperature profiles of the stratified base flow remain unchanged. However, the perturbation temperature, e, must now be identically zero. The stratified perturbation momentum equations are then uncoupled from the energy equation, becoming identically Eqs. (16). Similarly, setting Gr O0 implies AT O; Eqs. (16) again result. Pr 0 O i) T = R= 0 Equations (30) become the governing perturbation equations of the Benard problem. ii) R -= 0 Equations (30) become the governing perturbation equations of the rotating Benard problem.

44 iii) T = 0 The limiting form of the base flow velocity profile is that of plane Couette flow, Eq. (19), and the base flow temperature profile remains Eq. (28). Setting T - 0 in Eqs. (30) results in the governing equations in this limit which, on setting y = 90~ (c.f., Eqs. (17)), reduce to the modified Orr-Somerfeld equation coupled with the energy equation which governs plane Couette flow, heated from below. The limiting forms of Eqs. (30) are not pursued as T -+ o as they were for the unstratified flow regime. It is not expected that the results obtained from Eqs. (30) will approach the results obtained when T + A. This conclusion is based on the results of Homsey and Hudson (1971) who found that, in considering the rotating Benard problem, as T + A, the "excess centrifugal force" terms must be included and the effect of the side walls must also be considered. Equations (30) (or Eqs. (16), without stratification) represent an eigenvalue problem. The solution to the eigenvalue problem is considered next.

CHAPTER 3 NUMERICAL ANALYSIS —SOLUTION TO THE EIGENVALUE PROBLEM The method chosen to reduce the eigenvalue problem of the stratified flow, Eqs. (30), to an algebraic eigenvalue problem is discussed below. The same method was used for the unstratified flow. In solving the algebraic eigenvalue problem, two methods were used, one method for the case of stationary disturbances and another for the case of nonstationary disturbances. In the case of the unstratified flow, the critical disturbance was always nons tationary. Equations (30) represent an eigenvalue problem for which the eigenfunctions 4, u, and e are unknown. For fixed values of R, T, Gr, Pr, a, and y, solutions for P, u, and e are admitted only for certain values of c = c + ic i, the eigenvalue. In the solution method applied here, the unknown eigenfunctions are replaced by finite difference approximations and solutions, c, of the resulting algebraic eigenvalue problem are sought. The method is quite similar to that used by Lilly (1966) to study the stability of the Ekman layer. The finite difference representations of ~, u, and e are discussed first. The methods of solving the algebraic problem are discussed next, and finally the procedure followed in the parametric study and the convergence characteristics of the method are discussed. 45

46 3.1 FINITE DIFFERENCE REPRESENTATION The finite difference representation is presented in detail in Appendix II; however, a short discussion of this representation follows below. The eigenfunctions ~, u, and e are approximated by dividing the z domain into N intervals (N + 1 nodes) and expanding the variables in Taylor series around the interior nodes. The first and third of Eqs. (30) are expanded around N - 1 nodal points, the second around N half nodal points. The eigenfunctions 4 and 0 are defined only at nodal points and u is defined only at half nodes. In this manner, using a centered difference scheme, Eqs. (30) (including the boundary conditions) are represented by a system of 3N - 2 algebraic equations in 3N - 2 unknown eigenvalues (in the case of the unstratified flow there results a system of 2N - 1 equations). At a point of neutral stability, only the lowest eigenvalue is of interest. The Taylor series expansions take the form, after rearranging, 6 3 4 D = j 24 3 + -o(A )... (31) 2 2 A 4 4 D q = 6j - 6 - bj + -O(A )...+ j 12 6 D4,, _52<_ A <66i + ~O(A +... D4'p - 52j - 2- 4 where A N,, = j+j/2 6j a 2j = 6(6'j), etc.

47 To minimize truncation error, the terms of O(A ) (hereafter referred to as the remainder terms) are retained in all expansions except those of the highest order derivative in each equation. The remainder terms are not included in the approximate representation of the highest order derivatives because the difference equations would then be of higher order than the differential equations they are to represent. The approximate solution could not then be expected to behave in the same manner as the exact solution. The base flow velocity profile quantities appearing in Eqs. (30) are evaluated numerically at the nodal point (or half nodal point) in question (see Eqs. (II-10)). The resulting algebraic eigenvalue problem of order 3N - 2 can be represented in matrix form as: -1 [A + iB] - I = (32) where I = unit matrix X = the complex eigenvalue = -ia(c + ici) The elements of the submatrices A, B, and C are given in Appendix II by Eqs. (II-7) through (II-11). The solution of this eigenvalue problem is considered next. 3.2 ALGEBRAIC EIGENVALUE AND DETERMINANT SOLUTIONS The eigenvalues, X's, of the matrix defined in Eq. (32) are derived from the temporal description of the disturbance which is of the form e. The disturbances at neutral stability (r = aci = 0) fall into two categories; stationary (Ai = -iacr = O) and nonstationary (Xi # 0).

48 When the disturbance is nonstationary, the full eigenvalue problem, Eq. (32), must be solved. The method used was to find all the eigenvalues of the matrix of complex elements, Eq. (32), by converting it to an upper almost triangular form (Hessenberg form) by similarity transformations and then using the LR iteration, described in Wilkinson (1965), on the transformed matrix. Very efficient Fortran-IV versions (which use very little complex arithmatic) of the programs COMBAL, COMHES, COMLR originated in ALGOL by Parlett and Reinsch (1969) and Martin and Wilkinson (1968a, 1968b) are available. The Fortran IV versions were translated by a "task force" associated with the Argonne National Laboratories and are available in The University of Michigan computer system (MTS) under MATH:LIB. When the disturbance is stationary, then X - 0 and the determinant A + iBI - 0 (33) The zeros of this determinant are found by determining the sign change of the real part of the determinant. Consideration of the sign of the imaginary part of the determinant was found not to be necessary since the magnitude of the imaginary part was always much smaller than the real part. The method was frequently checked by computing the eigenvalues. This provides considerable computer economies as compared to the method described above which calculates all 3N - 2 eigenvalues of Eq. (32). The determinants were evaluated by writing a complex version (which uses only real arithmatic) of an available routine described by Moler (1971). Moler's routine uses Gaussian elimination and partial pivoting in an LU decomposition which completes the interchanges in U but only partially completes the interchanges in L.

49 The conventional routines using Gaussian elimination involve operations on rows of a matrix. Moler's routine involves operations on the columns which is the conventional method of storing twodimensional arrays in Fortran. This substantially reduces both the subscript calculations and repetitive addressing which involves locations separated from each other by a large increment. In a multiple user or virtual memory system, such as The University of Michigan's IBM System/360-67, reduction of repetitive addressing from widely separated locations results in a large reduction of paging operations. Moler reports a reduction of 32 percent in elapsed time on a 200 x 200 matrix of real elements compared to the elapsed time of the same routine which uses the traditional row operations. A version of the conventional routine (DECOMP) for matrices of real elements is found in Forsythe and Moler (1967). 3.3 CONVERGENCE AND NUMERICAL PROCEDURE Since the governing equations were reduced to a set of ordinary differential equations, the only error involved in the method is the truncation error associated with truncating the Taylor series expansions, Eqs. (31). It is evident from Eqs. (31) that as N + X, the finite difference approximation converges to the exact solution (assuming the computer round-off error is negligible). The method used to determine the convergence of the numerical results was to obtain the results at increasing values of N, the number of intervals, until the results obtained did not change appreciably. The magnitude of the acceptable change in the results with increasing N and the parameters involved will be discussed below.

50 As indicated in paragraph 3.1, the remainder terms (terms of 0(A )) were retained in the finite difference approximation of all terms except the highest derivatives in order to speed convergence. The contributions of the remainder terms to the elements of the matrix can be seen in Appendix II, Eqs. (II-8), as the underlined terms. The remainder terms were programmed to be included or excluded at will. Although a definitive study with and without remainder terms was not made, preliminary comparisons showed that, at the lower range of N used (15 < N < 25), inclusion of the remainder terms did effect computer efficiencies. The convergence of the value of a parameter of interest (such as Gr at neutral stability) was improved by as much as 3 percent. This improvement was usually more significant than the improvement obtained by increasing N to such an extent that the cost of operation was doubled. The improvement obtained by including the remainder terms is at virtually no cost. The remainder terms were therefore included in the study described below. Before the convergence of typical numerical results obtained is presented, it is useful to consider the nature of the parametric study. This requires a description of the numerical procedures followed. Two procedures were used, one for stationary disturbances and one for nonstationary disturbances. These procedures are discussed below for the specific case of stratified flow. The procedure followed in the numerical study of the unstratified flow was essentially the same as the procedure used when the disturbance of the stratified flow was nonstationary.

51 For both stationary and nonstationary disturbances, the parameters of the stratified flow are Pr, T, R, Gr, a, y, and X where I is the complex eigenvalue defined by Eq. (32). As stated previously, only neutral stability (r = 0) was investigated. If i = 0, the disturbance is stationary and if X # 0, the disturbance is nonstationary. In the discussion that follows, it is important to note that the parametric study was accomplished in conversational mode with the computer and with the aid of graph paper. When the disturbance was stationary (X = 0), the zero's of the determinant, Eq. (33), were found by searching for a sign change of the values of the determinant by programming the computer to iteratively calculate the values of the determinant using the method of secanting. The last two values of the determinant calculated (at constant Pr, T, R, a, and y) were required to bracket the zero (be of opposite sign) to within an interval of about 1 percent of the value of Gr. Thus, the Gr at neutral stability for constant Pr, T, R, a, and y was found. The initial values of a and y taken were determined by knowing the various limiting cases of the problem. This process would be repeated for various values of a (thus generating a partial neutral curve at constant Pr, T, R, and y) until the approximate value of a which resulted in the lowest Gr for these values of Pr, T, R, and y could be recognized. Using the new value of a, a partial neutral curve would be generated as a function of y (with Pr, T, R, and a constant) until the approximate value of y, which resulted in the lowest Gr at this Pr, T, R, and o, could be recognized. With these improved values of a and y, the above

52 procedure would then be repeated until the critical point was found (i.e., the values of a and y were found which resulted in the lowest Gr at constant Pr, T, and R). The above procedure is appropriately described as a double interpolation scheme. The double interpolation scheme was repeated for various values of Pr, T, and R, to generate critical curves. When the disturbance was nonstationary (i # 0), neutral curves could not be generated in a direct manner, such as described above for the case of stationary disturbances. The traditional approach for locating critical points has been to generate curves of constant growth rate, Ar. In the conversational mode with the computer, it was found that generating curves of constant growth rate was not necessary. The procedure followed was to roughly generate a trace of the growth rate, Xr, as a function of a at constant Pr, T, R, Gr, and y until the region of maximum growth rate (or minimum damping) could be recognized. In this region of the o plane, another growth rate curve would be roughly generated at a new value of R. From these two curves of Ar at constant R's, both the value of ca which resulted in the minimum R and the minimum value of R at neutral stability could be estimated. With these new values of a and R, a similar procedure would be initiated in the y plane. This whole procedure would be iterated until the critical point was found. That is, the values of a and y which resulted in the lowest R at neutral stability was found at constant Pr, T, and Gr. This procedure is appropriately described as a triple interpolation scheme. The triple interpolation scheme would be repeated for various values of Pr, T, and Gr to generate critical curves.

53 It should be noted that the triple interpolation scheme applies directly to the unstratified flow by observing Pr and Gr are not relevant parameters in that case. Frequently, when applying the triple interpolation scheme to the stratified flow, the roles of the Gr and R would be interchanged. Whether the growth rate curves would be generated at constant R or constant Gr depended upon whether the initial estimates of a and y were derived from the limiting case R + 0 or Gr -+ 0 with Gr used as R -+ 0 and R used as Gr + 0. It is to be emphasized that the above schemes were substantially aided by knowing the many limiting cases of the respective problems (c.f., paragraphs 2.2.4 and 2.3.2). From the above discussion, it is evident that convergence considerations must be given not only to the value of the critical parameter pursued (Rcrit or Gr crit) but also to the numerical values of a and y at which these critical values occur. In view of this, the convergence was studied by repeating the double or triple interpolation schemes at increasing N. In the convergence results presented below, it will be noticed that acceptable values of the a and y, at which the critical parameter (R or Gr) occurred, were obtained at lower values of N than was necessary to obtain acceptable values of the critical value itself. It is of interest to note that Korpella (1972) observed the same behavior using weighted residual methods on a study of natural convection in an inclined slot, heated from the lower inclined side. This result was used in the present study to effect computer efficiencies by using the double or triple

54 interpolation scheme at the lowest level of N (15 i N i 25) which resulted in accurate values of the a and Y at which the critical parameter occurred. The critical parameter would then be obtained by using these values of a and y and further increasing N (25 < N < 40). Observing that the numerical results converge to an apparent final value, on increasing N, does not in itself constitute a proof of convergence. The computer programs were used to. obtain selected numerical values for the known limiting cases of the flow. For convenience, the program which finds the zeros of the determinant, Eq. (33), used in the double interpolation scheme (described above) will be referred to as the determinant program. The program which finds all of the eigenvalues of Eq. (32) used in the triple interpolation scheme (also described above) will be referred to as the eigenvalue program. Table I contains the results of the study of the classical Benard problem with T = R = 0 used as input to the determinant program. The results of Table I were also obtained from the eigenvalue program with comparable convergence behavior. In Table I, the critical value of ~R obtained is 1700.016 and the a at which crit occurs is 3.115. Reid and Harris (1958) found crit = 1707.762 crit crit at a = 3.117 in their study of the Benard problem. The disagreement in cR is <0.45% and in a <0.1%. It is evident that 1Rhas not crit quite converged to its final value. The apparent convergence to within 0.5% of a final value was considered acceptable in the present study. It should be noted in Table I that the final value of c is already correctly indicated at N = 20. The R of Table I

55 were found by setting Pr = 1.0 and finding the Gr. It is known that the cit of the Benard problem is independent of Pr; this behavior (in the present solution method) was verified by finding results comparable to those in Table I at other values of Pr. Table II contains the results of a study of the rotating Benard problem at T = 11.8 and Pr = 5.0. These results were obtained by setting R = 0 in the determinant program; comparable results were again obtained with the eigenvalue program. To compare the present results with the known results (c.f., Chandrasekhar, 1961), the parameter defined herein as the Taylor number, T, must be converted to the traditional Taylor number, Ta. The conversion is Ta = 4T2; therefore, T = 11.18 corresponds to Ta = 500. The critical value of Gr and the value of a at which Grcrit occurred (Table II) are 386.250 and 3.30, respectively. This results in ~rit = Pr * Grcrit 1931.25. These values compare to 1crit 1940.0 and a = 3.30 given by Chandrasekhar. The disagreement in Rcrit is again <0.5%. The rit of the rotating Benard problem is not decrit pendent on Pr for Pr > 0.677. The fact that the above results were obtained by arbitrarily choosing Pr = 5.0 verifies the lack of Pr dependence in the present solution method for this case. The eigenvalue program was used to obtain results for the rotating Benard problem at Pr = 0.025 and T = 500 (Ta = 106). In this case, the critical disturbance is overstable (c.f., Chandrasekhar, 1961). The results obtained also compare to Chandrasekhar's to within 0.5%.

56 TABLE I CONVERGENCE, BENARD PROBLEM - RAYLEIGH NO., ~T, AS A FUNCTION OF WAVE NO., a, AND NO. OF INTERVALS, N 8 15 20 25 30 35 3.080 1657.973 13.090 1657.918 3.100 1657.923 3.105 1657.937 1679.317 1690.478 3.110 1657.960 1679.301 1690.464 3.115 1657.997 1679.290 1690.457 1696.210 1700.016 3.120 1658.062 1679.351 1690.479 TABLE II CONVERGENCE, ROTATING BENARD PROBLEM - GRASHOF NO., Gr, AS A FUNCTION OF WAVE NO., o, AND NO. OF INTERVALS, N [T = 11.18 (Ta = 500) and Pr = 5.0] N t 15 20 25 30 35 3.20 378.365 3.25 377.950 382.356 3.30 377.801 382.170 384.262 385.455 386.250 3.35 377.911 382.246 3.40 378.276

57 In addition to the results discussed above, some selected results were obtained for convective pCf. This was accomplished by two methods. First, the linear base flow velocity profile of convective pCf was programmed into the eigenvalue program and calculations were made with T E 0. Second, the eigenvalue program was used retaining the velocity profiles of the present problem, Eqs. (18), and obtaining results at very small values of T (T = 0.1). The second method was prompted by noting that the limit of the present base flow velocity profile, Eq. (18), is the linear velocity profile of pCf (c.f., Eq. (19)). The results obtained by both of these methods were compared to the results given by Gallagher and Mercer (1965) in their study of convective pCf. The selected results obtained with the linear velocity profile compared to the results of Gallagher and Mercer to within 0.5% and the results at T = 0.1 to within 2%. Typical convergence results concerning the stratified and the unstratified flows of the present study are now presented. Table III contains typical results of the double interpolation scheme (described above) for the stratified flow regime. The Grashof number at neutral stability is given as a function of a, y, and N in the vicinity of the critical point. Note that, in Table III, as N was increased from 20 to 25, the values of a and y, used to obtain results, were changed. This was necessary because the initial values of a and y used were not those which resulted in Grcrit. That is, a search of the y plane was initiated by finding values of Gr at various y using a = 5.58 at N = 20 and 25. The

58 minimum Gr occurred at Y = 74.00 for both N = 20 and N = 25. The search of the a plane, at N = 20, was calculated using Y = 76.00 and the minimum Gr occurred at a = 5.56. In the search of the a plane at N = 25, the information from the y plane search was available and, therefore, Y = 74.0~ was used. Again, the minimum occurred at x = 5.56. Therefore, a second search of the y plane, at N = 25, was calculated using a = 5.56 with the minimum Gr remaining at y = 74.0~. Therefore, the values of a and y at which the minimum Gr occurs are 5.56 and 74.0~, respectively. Furthermore, these values of a and y were indicated at N = 20. As is shown in Table III, the values of a = 5.56 and Y = 74.0~ were then used to calculate Gr at N = 30 and N = 35 to obtain Grcrit Considering the results at a = 5.56, Y = 74.0~, and starting at N = 25, it can be seen that the solution converges from above. This was not always the case. It is not known why the solution converged from above part of the time and from below the rest of the time; however, it is interesting to note that Korpela (1972), using weighted residual methods on a study of natural convection in an inclined slot observed the same behavior. Figures 5 and 6 present typical convergence behavior of the triple interpolation scheme (described above) using the eigenvalue program as applied to the unstratified flow. It should be noted that when the critical disturbance of the stratified flow was nonstationary, the convergence behavior therein was comparable to the convergence behavior of the unstratified flow. Therefore, the results discussed below are indicative of all the nonstationary results. In Figs. 5 and 6, the growth rate,?r, is plotted as a

59 TABLE III CONVERGENCE, STRATIFIED FLOW - GRASHOF NO., Gr, AS A FUNCTION OF WAVE NO., a, ORIENTATION, y, AND NO. OF INTERVALS, N (Pr = 1.0, T = 100, R = 100) N a 20 25 30 35 (y = 76.00) (Y = 74.00) (y = 74.00) (y = 74.00) 5.50 9974.992 5.53 9973.379 5.54 9973.035 9933.347 5.56 9972.688 9932.886 9914.460 9905.041 5.58 9972.759 9933.781 5.60 9973.286 5.62 9934.141 5.65 9976.354 rI~~~~ N * 0 I Y 20 -- 25 --- (a = 5.58) (a = 5.58) (a = 5.56) 73.0 9974.821 9935.034 9934.874 74.0 9972.275 9932.973 9932.886 75.0 9972.361 9933.033 9933.116 76.0 9972.759 9935.432 77.0 9940.277 TABLE IV CONVERGENCE, UNSTRATIFIED FLOW - CRITICAL R AS A FUNCTION OF NO. OF INTERVALS, N N | T a/yo g 20 25 30 35 40 100 3.24/-20.0 61.0 58.5 57.0 56.5 56.3 100 5.10/10.0 105.8 108.2 109.6 110.5 35 3.15/18.0 158.6 160.0 161.0 161.4

60 function of y (Fig. 5) and a (Fig. 6) for T = 100 and at the values of R = R//I and N given. The only curves shown from the triple interpolation are those for which the values of Y and a which result in the critical value of R can be recognized. That is, at N = 30, a neutral point (Ar = 0) is shown to occur at Y -20.0~ (Fig. 5) and a 3.24 (Fig. 6). At N = 30 and R = 57.0, any other values of Y and a resulted in a growth rate curve which remained in the negative Xr domain. Thus, the neutral point shown is the critical point at N = 30. Of particular interest in Figs. 5 and 6 is that at both N = 25 and N = 30, the point of minimum damping occurs at very nearly the same values of y and a.. In Fig. 5, the minimum point for N = 25 occurs at Y = -20.0~ and for N = 30 at y = -20.2~. In Fig. 6, the minimum point for N = 25 occurs at a = 3.23 and for N = 30 at a = 3.24. The shift in y should be considered with regard to the range of y, which is 180~. From this viewpoint, the shift in Y is <0.1% of the range. The shift in a from N = 25 to N = 30 is <0.3% of the value of a. In the present study, the values of Y and a at N = 25 were therefore considered accurate. The values of Y = -20.0~ and a = 3.24 (for T = 100) were then used to obtain results of R at N = 35 and N = 40. These scrit results are given as the first entry in Table IV. The final results for another instability at T = 100, as well as the final results for T = 35, are also given in Table IV. The reason for two sets of results at some of the values of the Taylor number, T, will be discussed in the next chapter. All of the results in Table IV are considered to have converged to within 0.5% of their final values.

61 It is evident in Table IV that the solution converges faster at T = 35 than at T = 100. In general, the convergence, on decreasing T from T = 100, improved with one exception. The exception was that, as T became increasingly smaller, the RE necessary to cause instability became increasingly larger, eventually causing the solution method to fail (in the sense that it became prohibitively expensive). Since this behavior was expected and is related to the nature of the instability at T -+ 0, further discussion of this behavior will be included in the next chapter (Results and Conclusions). It is sufficient to note, at this point, that this region of poor convergence caused by large R crit did not unduly restrict the parameter range which could be studied. This concludes the discussion on numerical procedure; the results and conclusions are discussed next.

62 -0.5 2 0 4 1 -0 3 Li -0.2 2 UI 2'r -0.1 3 0. 0 C3 +0.1 - 1-1 N=25, R =58.0, a=3.24 +0. 2 1 2-2 N=30, R =57.0, a=3.24 +0.3 _ -16.0 -18.0 -20.0 -22.0 -24.0 ORIENTATION, y, degrees Figure 5. Convergence, Unstratified Flow - Growth Rate, Xr, Near Critical Point as a Function of Orientation, y, and No. of Intervals, N, (T = 100). 1 -0.2 2 20,2 -0.1 - LUJ F0.0 o zC 11-1 N=25, R =58. 0, y=-20.0 2-2 N=30, R =57.0, y=-20.0O +0. 1 3. 00 3.10 3.20 3.30 3.40 WAVE No., a Figure 6. Convergence, Unstratified Flow - Growth Rate, Xr, Near Critical Point as a Function of Wave No., O, and No. of Intervals, N, (T=100).

CHAPTER 4 RESULTS AND CONCLUSIONS The numerical results and conclusions for the unstratified and stratified flow regimes are discussed separately. 4.1 UNSTRATIFIED FLOW Table V is a summary of the numerical results obtained from the parametric study of the unstratified flow regime, Eqs. (16). For all T > 48, two types of instability resulted. Following Greenspan's (1969) discussion of the stability of the Ekman layer, these are denoted "A" and "B" waves. Faller and Kaylor (1966) and Lilly (1966) refer to these waves as traveling and "almost" stationary waves, respectively. The latter waves (class B) are nearly stationary with respect to the rotating disk. In the rotating frame of the reference chosen here (rotating at the mean rotation rate at mid-point), both waves are traveling. The last of Eqs. (22) gives the transformation of the present wave speed to a rotating reference at a disk which corresponds to the Ekman studies. In the presentation of the results, the parameters R and a, also defined by Eqs. (22), have been used. This permits a direct comparison with the results of Lilly (1966) for the stability of the Ekman layer. Figure 7 presents the results for the critical (lowest) R as a function of Taylor number. Figure 8A presents the ac for which the critical R occurred as a function of T. ~~~ S~~~ Figure 8B presents the y for which the critical R occurred as a function of T. Lilly's results for the "A" and "B" waves are shown 63

64 TABLE V CRITICAL RESULTS, UNSTRATIFIED FLOW'Type! cr T Type R R y /R c'WaveEE 100.0 "A" 563 56.3 3.24 0.324 -20.0 0.226 0.568 100.0 "B" 1105 110.5 5.10 0.510 10.0 0.238 0.065 70.0 "A" 533 62.5 2.75 0.329 -16.5 0.226 0.472 70.0 "B" 931 111.3 4.40 0.527 10.5 0.239 0.057 50.0 "A" 774 109.0 1.83 0.259 - 9.0 0.236 0.394 50.0 "B" - 856 121.0 3.90 0.550 13.0 0.247 0.023 48.46, "A" t 855 122.8 1.70 0.244 - 8.0 0.238 NC* 48.461 "B" i 851 i 122.0 3.85 0.552 13.0 0.248 NC 35.0 "B" 956 1 161.4 3.15 0.531 18.0 0.267 NC 20.0 J "B" 4805 j1074.0 1.15 0.258 35.0 0.359 NC Not Calculated Not Calculated

65 1000 900 800 700 600 500 400 300 - ~r 200 / IIBm 100 90 80 70 60 50 40 0 20 40 60 80 100 TAYLOR No., T Figure 7. Critical Curve, Unstratified Flow - R as Function of T. crit

66 0.6 itBU O. 5.] 4 - _ _..Bu 0. 45 C- 0 A "A" d 0'3 1,. -- 0.2 0.1 0.0 0 20 40 60 80 100 Figure 8A TAYLOR No., T 90 80 - 70 60, 50 40 30 lo B 0 -10 -20 "A" -30 0 20 40 60 80 100 Figure 8B. TAYLOR No. T Figures 8A and 8B. Wave No., cx, and Orientation, y, at R as Function of T. crit

67 in these figures as dashed lines. It is evident that the present numerical solution does indeed converge to the Ekman layer results as T becomes large. At T = 100, the present solution already results in an R within 3% of Lilly's value for the "A" wave and within 4% for the "B" wave. As can be seen in Fig. 7, the "A" wave does not appear below a value of T of about 47. It is known from Lilly's work that on increasing R, the "A" wave will disappear again. That is, Lilly found that at R = 500, the "A" wave did not exist but at R = 150 it did. Thus, according to the linear theory, the domain of the "A" wave in the R vs T space is limited. There results a locus which has an R (shown qualitatively as a dashed line on Fig. 7) Emax above which the "A" waves do not grow and a R below which the ecrit "A" waves do not grow. A thorough search at T = 45 was performed. The vestige of the "A" wave was found but it was always damped. An example of this is shown in Fig. 9 where the damping rate, ac, is plotted against R, at T = 45, caE = 0.253 and y =-7.0~. The aE and y chosen are those which result in a minimum damping at R = 193.8 (R = 1300) and T = 45. The damping first decreases on increasing R and then increases. Since the minimum damping occurs at an a and Y which would be a continuous extension of the a and E: Y curves which resulted in R (Figs. 8A and 8B) for the "A" scrit wave, and since the growth rate could not be forced to become positive by increasing the Reynold's number, these results have been interpreted as the vestige of the "A" wave.

68 T = 45.0 -4.0 - = 1.7 = -7.0~ 0 2.0 rE. -1. O CD 0.0 I, I, I' I'! 150 170 190 210 230 250 R =RT- 2 Figure 9. Vestige of the "A" Wave.

69 The result that the "A" waves do not continue to be unstable as T becomes small is consistent with Lilly's interpretation of the "A" wave as a Coriolis-induced instability. This follows simply from the fact that T is a measure of the Coriolis acceleration. In considering the closed domain of the "A" wave, as well as the occurrence of the "B" wave at values of T > 47, it is logical to question the validity of the linear theory above the "A" wave R locus. However, Faller's (1963) experimental results show that both the "A" waves and the "B" waves occurred simultaneously in regions of overlapping radii in the Ekman layer (but started at different radii). Furthermore, the linear theory, as used by Lilly (1966), successfully predicts both the "A" and the "B" wave R.rt In considering the validity of the present results in this ccrit region of the R vs T domain, it is useful to first consider the nature of the results obtained. In the numerical work, the unstable modes always occurred in conjugate pairs with only one pair unstable for any given set of parameters. This conjugate pair is interpreted physically as indicating two waves with equal but opposite wave speeds, one in the positive z domain and one in the negative z domain. Thus, whenever an instability occurred, the solution always involved only one mode in each z domain. As is shown in Fig. 8, the "A" and "B" waves always occur at different wave numbers and orientations. This would inhibit interaction of the waves and, on that basis, one would expect nonlinear terms to be less important than if interaction were expected. It is also evident in Fig. 7 that, in general, the

70 "A" and "B" waves occur at different R. This can be interpreted (through U0) as indicating instabilities occurring at different physical radii. The instabilities observed by Faller appear to be quite localized at different positions in the z domain, in the vicinity of the disk. Therefore, the basic nature of the flow was not destroyed (after instability occurred) in a rather large range of R, i.e., radius. On the basis of the above, it is expected that in an experiment conforming to the configuration considered herein, both "A" and "B" waves would be observed for T > 47 and also that above some R,which would depend on T, the "A" wave would no longer exist. Referring to Fig. 7 at T! 48 and R 122, both the "A" and the "B" waves appear to be equally probable. At this Taylor number, one might then expect to see both waves begin at the same radius. This may not happen, however, if the growth rate of one of the waves increases faster than the other, as a function of R Figure 10 is a graph of the growth rates vs R for both the "A" and "B" waves in the immediate vicinity of the intersection of the two curves in the R vs T space. The growth rate of the "B" wave is shown to increase much faster than that of the "A" wave. Thus, the "B" wave would appear to be the more probable in the vicinity of T = 48. Also, the solution (c.f., Fig. 7) predicts the "A" wave to be damped out again at larger radii. From the above discussion, it is concluded that in an experiment, it is possible that one would not observe the "A" wave at values of T as low as predicted by the linear theory. This possible discrepancy is, however, minimized by the steepness of the slope of the REcrit vs T "A" wave locus in this

71 IIB LL 2.0o 1 T = 48.64 -1.0 IS.~ —- ~~=3.85 I:~ | " A " a =1 7 0 y 13.00 c= 850 855 860 REYNOLD'S No., R Figure 10. Growth Rates as Functions of Reynold's No. in Vicinity of "A" and "B" Wave Intersection.

72 vicinity. Interaction of the waves could alter this conclusion but, as discussed above, interaction is not expected to be serious. Referring again to Fig. 7, it is evident that as T decreases from T - 48, the Rit increases very rapidly. This is consistent with the limiting form of the governing equations as T > 0. It was demonstrated in paragraph 2.2.4 that as T -+ 0, Eqs. (16) reduce to the Orr-Somerfeld equation for plane Couette flow (pCf) and, as discussed in paragraph 1.2.1, the linear theory appears to predict that pCf is stable with respect to all infinitesimal disturbances. In terms of the present problem, it is therefore anticipated that at T - 0, the critical Reynold's number, Rcrit" should become very large. Referring to Table V, Rcrit has increased to 4805 at T = 20. Considering the parameter R since Rcrit must become large as Rcrit T + 0 R = +00 as T + 0. ccrit Hopf (1914) and Grone (1954) studied pCf by asymptotic methods and conclude that as aR + a, the wave speed, c /R -+ +1. Thus, as T - 0, the wave speed of the "B" wave should approach ~1, while Y approaches +900. Referring again to Table V, the numerical results for the "B" wave are consistent in this respect with the wave speed increasing rapidly from c r/R = 0.247 at T = 50, to cr/R = 0.359 at T = 20. An extrapolation of the y curve of Fig. 8B to T = 0 results in Y = 90~. Continuing the solution to small T (T < 20) becomes increasingly more difficult because numerical convergence becomes poor. This can be explained from two points of view.

73 i) The "B" wave is caused by vorticity distribution. Asymptotic analysis of this type of instability, valid at large aR, results in solutions of the form, (c.f., C. C. Lin, 1966): = e~- eIf (z) +...} where Q = /i(u c) dz. YC These solutions oscillate rapidly as z is varied. Thus, at large aR, a high spatial resolution of the eigenfunction is required. ii) The second point of view is similar to the point of view taken by Gallagher and Mercer (1962) using the Galerkin method in a study of pCf. In terms of the present eigenvalue problem, large aR manifests itself as follows. Reconsider the matrix representation of the algebraic eigenvalue problem, Eq. (32), which for the unstratified flow, can be represented as The highest derivatives are all represented in the submatrix [A] with dominant terms of O(N ). As aR increases, it is evident that N (the number of intervals) must also increase or the solution will be dominated by the [B] submatrix. Thus, the matrix solution with N constant behaves, on increasing aR, as if the governing equations were singular.

74 In the parameter range of the present analysis, the convergence was watched carefully and, as discussed in paragraph 3.3, the numerical results given appear to have converged to within 0.5% of their final value. In summary, it has been shown that the unstratified results converge to the results of Lilly (1966) for the Ekman flow as T becomes large and exhibit the proper behavior, that of pCf, as T + 0. The unstratified flow exhibits both of the instabilities observed in the Ekman layer, through a range of T, with the "A" (Coriolis induced) wave the critical disturbance above T = 48 and the "B" (vorticity induced) wave the critical disturbance below T - 48. Furthermore, it was shown that the "A" wave will not grow (become unstable) below T = 47. The stratified flow is considered next. 4.2 STRATIFIED FLOW Tables VI, VII, and VIII contain a summary of the critical results obtained from the parametric study of the stratified flow regime, Eqs. (30). Numerical results were obtained for Prandtl numbers, 6.0, 1.0, and 0.025 which closely correspond to water, air, and mercury, respectively. In the discussion that follows, the Grashof number, Gr, will consistantly be used as the critical parameter. That is, the results given in Tables VI, VII, and VIII (for Pr = 6.0, Pr = 1.0, and Pr = 0.025, respectively) are interpreted as follows. At the Pr, T, and R (or R ) considered, the numerical value of the Grashof number given is interpreted as Grcrit; i.e., if Gr is increased

75 TABLE VI CRITICAL RESULTS, STRATIFIED FLOW (Pr = 6.0) Pr = 6.0 jT;~~~~~~~~~~ r t T R i R Gr a c /R 100.0 100.0 10.0 2638.0 1 5.80 17.5 0 100.0 1 200.0 20.0 3854.5 1 5.65 - 25 100.0 300.0, 30.0 5154.8 6.10 -13.0 0 100.0 500.0 50.0 7435.5 6.50 -21.0 0 -- 4. —- ------ 50.0 50.0. 7.07 1242.7 4.65 68.0 0 50.0 j 200.0 28.28 3701.5 4.85 49.5 0 50.0 500.0' 70.70 8334.9 5.55 40.5 0 50.0 730.0,103.2 1 5000.0 1.90 f -10.0 0.231 ------------ ------------------ --— 1 —--------- 20.0 50.0 11.2 683.2 3.40 -61.5 0 20.0 200.0 44.8 2982.4 4.50 -71.0 0 20.0 500.0 111.8 7283.5 5.05 -75.0 0

76 TABLE VII CRITICAL RES1ULTS, STRATIFIED FLOW (Pr = 1.0) Pr = 1.0 T 1 T R R_ l Gr a y c /R 100.0 100.0 10.0 9905.0 5.56 i 74.0 0 100.0 200.0 20.0 12540.0 5.08 j 67.0 0 100.0 300.0 30.0 16199.0 5.33 50.0 0 100.0 400.0 40.0 19268.0 5.56 37.0 0 100.0 500.0 50.0 21998.0 5.68 31.0 0 100.0 559.0 55.9 15000.0 3.35 -20.0 0.211 70.0' 100.0 11.9 7208.0 4.80 86.0 0 70.0 300.0 35.8 14666.0 4.70 66.0 0 70.0 500.0 59.7 22230.0 5.10 51.0 0 ~. ~~....... —- --- - --- —. m. —.I —-- - - 50.0 50.0 7.07 4777.5 4.55 -82.0 0 50.0 200.0 28.28 9622.0 3.70 -87.0 0 50.0 500.0 70.70 25247.0 5.00 69.0 0 50.0 756.0 112.0 7000.0 1.88 -10.0 0.228 50.0 738.0 104.3 14000.0 1.92! -10.5 0.225 __________________________ ----------- ----------------— W —_-__ —__ 35.0 50.0 8.45 i 3660.0 4.05 -69.0 0 35.0 200.0 33.8 9370.0 2.90 -69.0 0 35.0 500.0 84.5 27663.0 5.30 87.0 0 35.0 947.0 160.0 5000.0 3.13 18.0 0.266 20.0 50.0 11.2 2610.0 3.50 -49.0 0 20.0 100.0 22.4 3709.0 3.30 -48.0 0 20.0 200.0 44.8 8675.0 2.50 -49.0 0 20.0 300.0 67.1 15176.0 5.00 -58.0 0 20.0 X 500.0 111.8 25415.0 5.85 -64.0 0

77 TABLE VIII CRITICAL RESULTS, STRATIFIED FLOW (Pr = 0.025) Pr = 0.025 C T R R 7 Gr 1 a y, c /R 50.0 50.0 7.07 173800.0 2.93 -41.5 1 0.332 50.0 175.0 1 24.77 186200.0 2.70 -39.0 i 0.141 50.0 240.0 i 34.3 197300.0 2.51 -36.0 0.135 50.0 300.0 42.4 204900.0 2.41 -31.5 0.144 50.0 500.0 70.7 144900.0 2.25 -15.0 0.206 50.0 600.0 84.8 92600.0 2.09 -12.5 0.223

78 above the value indicated, the base flow is predicted to be unstable. The values of a and y given in Tables VI, VII, and VIII are the values of the wave number, a, and the orientation, y, r at which Gr occurred. The value of c /R given is the wave crit speed at Grit In the graphical presentation of the critical results (from Figs. 11 through 21), the parameter RE, defined by Eqs. (22), has been used. The reason for using Rc in the presentation of the stratified results will be discussed below. 4.2.1 Critical Curves Figures 11 through 21 are graphical presentations of the critical results of Tables VI, VII, and VIII. Figures 11, 14, and 17 present the critical Grashof number, Grcrit, as a function of R at Pr = 6.0, Pr = 1.0, and Pr = 0.025, respectively. The Grcrit are presented as curves of constant T. Figures 12, 15, and 18 present the wave numbers, a, at which Grcrit occurred. Figures 13, 16, and 19 present the orientations, y, at which Grcrit occurred. In Figs. 12 and 15, it will be noticed that there are some regions which are denoted by broken lines in an otherwise continuous curve. For example, the wave number curve for Pr = 1.0 (Fig. 15) and T = 35 is presented as a broken line in the regions 40 < R < 55 and 130 < R < 145. By referring to Table VII, it can be seen that the critical results are not given in these regions. Sufficient numerical results were obtained in these regions with small numbers of intervals (15 < N < 25) to present the c vs R curves, but the values of Grcrit were not pursued c ci

79 9 Pr 6.0 8 6 CD) c 5L4.. 4 -- / T=O 0 20 40 60 80 100 120 Figure ii. Critical Grashof No., Grcrit as Function of R - RT-'1/. (Pr = 6.0).

80 Pr = 6.0 N" I 90 / \ / / / R RT Figure 12. Wave No., a, (at Gr crit) as Function of R = RT 1/2 (Pr = 6.0)

81 -30 Pr = 6.0 -40 -50 -60-70-80 90 m 807060o.,SO 50' 40 30 20 10 0 2' -0-I ------ -10 -20 -30 4.. 0 20 40 60 80 100 120 -1/2 R = RT Figure 13. Orientation, y, (at Grcrit) as Function of = RT1/2 (Pr = 6.0)

82 44 42 40 - Pr = 1.0 3836 - 34 / 32 / 30 / o, 2- /8 26/, 24 22 20 - v18 16 - 14 - 12 / 10 8 - A, 6 T 0 0. 0 20 40 60 80 100 120 140 160 R = RT Figure 14. Critical Grashof No., Grcrit, as fumction of R = RT-2 (Pr: 1.0).

83 7.0 Pr = 1.0 ~ T =20 6.0 -, 5.0 4 = m3.0 2.01.0 0.0 0 20 40 60 80 100 120 140 160 R RT-1 Figure 15. Wave No., c, (at Grcrit) as Function of R RT-2 (Pr = 1.0).

84 -30 - -40 - Pr = 1.0 T - 20 -50 - -60T = 35_ — -70 - 90 80- NN 70 60o 5040 30 - 20 - 0-10 -20-, -30.... 0 20 40 60 80 100 120 140 160 R - RT/2 Figure 16. Orientation, y, (at Grcit) as Function of R: RT1/2. (pr: 1 0)

85 since, as discussed in Chapter 3, this would require obtaining the results at a large value of N (and, consequently, increase computer expense). The additional calculations did not seem warranted since the Grcrit vs Rc results appeared to be well behaved, allowing extrapolation through these areas with confidence. Therefore, the Grcrit vs RE curves are not shown as broken curves in these regions (c.f., Figs. 11 and 14). The y vs RE curves also appear to be well behaved and are not shown as broken curves in these regions (c.f., Figs. 13 and 16). The broken curves referred to above are to be distinguished from the vertical dashed lines of Figs. 12, 13, 15, and 16 which are discontinuities in the curves. These discontinuities will be discussed further below. In considering the results, it is helpful to recall the limiting cases of the problem as described in paragraph 2.3.2 and to consider these limits as they apply to the Grcrit - R space (e.g., Figs. 11, 14, and 17). For T = 0 and any Reynold's number, R, Grcrit must take the value associated with the Benard problem for which the critical Rayleigh number is crit = 1708. Therefore, Grcrit = 1708/Pr whenever T = O. This is shown in Figs. 11, 14, and 17 where the curves labeled T = 0 are horizontal since, as discussed in paragraph 2.3.2, the stratified pCf always becomes unstable at crit = 1708. The critical disturbances of stratified pCf are stationary with respect to the velocity at the mid-plane. For the case of no shear, R = 0 (and, therefore, RE = 0), Grcrit must correspond to that of the rotating Benard problem at crit

86 the Pr and T of interest. Thus, the Grcrit curves, with T the parameter, must intersect the ordinate at the value associated with the rotating Benard problem. The critical disturbances of the rotating Benard problem are always stationary for all boundary conditions if Pr > 0.677 and are nonstationary (overstable) for free boundaries if Pr < 0.677 and T > 11. For Gr = 0, the results of the unstratified flow, presented in paragraph 4.1, are applicable with the conclusion that the Gr crit curves of constant T must intersect the abscissa, Rc, at the R ~ escrit (of Fig. 7) associated with the T of interest. These disturbances must be either "A" or "B" type traveling waves, depending on whether T is greater than or less than T = 48, respectively. The effect of shear and rotation (as represented by R and T) on Gr, y, and c/ will now be considered, first at T) on Grcrit a, y, and c /R will now be considered, first at Pr = 1.0 and 6.0 and then at Pr = 0.025. The effect of Pr with other parameters held constant will be considered separately. 4.2.1.1 Pr = 1.0 and 6.0 Both stationary and traveling-wave instabilities are predicted at Pr = 6.0 and 1.0 (Figs. 11 and 14). There are two types of traveling disturbances. Since, at these two Prandtl numbers, the results are quite similar, the remaining discussion is in terms of just the Pr = 1.0 results (Figs. 14, 15, and 16). The region of stationary instability is indicated in Fig. 14 by that region of the constant T curves which starts at R = 0 and continues through the region for which Grcrit increases with increasing Re. The traveling disturbances are again denoted as "A" or "B" waves,

87 the "A" wave occurring above T = 48 and the "B" wave occurring below T 48. It is evident on Fig. 14 that both shear and rotation (as measured by RE and T) have a stabilizing effect on the stationary disturbances, but have a destabilizing effect on the traveling waves. These results are not surprising in view of the apparent causes of these instabilities. The stationary disturbances are of the same type as those which occur in the rotating Benard problem, i.e., bouyancy driven convection cells. It is evident from Fig. 14 that increasing T has a stabilizing effect for all R as long as the critical disturbance is stationary. It is of interest to note that this conclusion cannot be discerned from Grcrit vs R curves. The Gr - R crit crit curves are, in fact, confusing in this respect. This becomes evident by referring to the tables of critical results (choosing Table VII for Pr = 1.0) and noting the results given for Gr rit at constant R. For example, at R = 200 and comparing Grcrit at T = 100 to Gr crit crit at T = 50, the Grcrit are 12540 and 9622, respectively. This indicates a stabilizing effect of T at constant R. However, at R = 500, again comparing at T = 100 and T = 50, the Grcrit are 21998 and 25247, respectively. This indicates a destabilizing effect of T at constant R which is, obviously, the reverse effect. This is the reason for choosing to use R instead of R in graphically presenting the results. As stated above, the effect of increasing R is also stabilizing as long as the critical disturbance is stationary. This effect was, in a sense, also noted by Gallager and Mercer (1965) (referred to as GM hereafter), in their study of stratified plane Couette flow (pCf). The conclusion concerning the results of GM is

88 is based on the following. The Grcrit were monotonically increasing functions of Reynold's number, R. Since Squire's theorem is valid for stratified pCf, a nonzero R can be interpreted as representing a disturbance having a transverse component. R = 0 then represents a longitudinal disturbance (disturbance with wave axis aligned in the direction of flow). Although GM concluded that the critical disturbance was longitudinally oriented, it is also clear that if Grcrit is a monotonically increasing function of R, then, on choosing some particular orientation, increasing R has a stabilizing effect on that particular disturbance. It is in this sense that the present results are in agreement with the results of Gallager and Mercer, although Squire's theorem is not valid in the present problem. As already implied by the labeling of the traveling disturbances as "A" and "B" waves, these disturbances are of the same nature as those found in the unstratified flow presented in paragraph 4.1. This point is discussed further below. Figures 15 and 16 represent the wave numbers, a, and the orientation, y, vs R for which the Grcrit occurred. The transition crit from a stationary to a traveling disturbance (on increasing R ) is very evident in these figures, with a sharp discontinuity occurring in both a and y when the new type of disturbance becomes the most critical. As was noted previously, the critical curves at Pr - 6.0 follow the same general trends as discussed above for Pr = 1.0. The Pr - 0.025 results follow.

89 4.2.1.2 Pr = 0*025 At Pr = 0.025, numerical results were obtained for T - 50. The results were limited to one value of T because the critical disturbance was always nonstationary. This required finding all the eigenvalues of Eq. (32) and, as discussed in paragraph 3.2, this involved considerable computer effort. That the critical disturbance is always nonstationary is consistent with two of the known limiting cases of the problem, R + 0 and Gr -+ 0. As R + 0, the results of the rotating Benard problem apply. Chandrasekhar (1953) and Chandrasekhar and Elbert (1955) found the critical disturbance of the rotating Benard problem with free boundaries to be nonstationary for Pr = 0.025 and T > 11. As Gr + 0, the results of the present unstratified flow apply and it was -found in this case that the critical disturbance was always nonstationary (c.f. paragraph 4.1). Thus, in terms of the parameter R, the limiting cases of the stratified flow are nonstationary at both ends of the range of RE (i.e., R = 0 to the R5 at Gr = 0 which is the largest R5 considered). It is, therefore, not inconsistent that the disturbance is nonstationary for all Rc considered. Considering Fig. 17, it can be seen that the effect of increasing R is, at first, slightly stabilizing (i.e., Grcrit increases) and then destabilizing. The transition between these two effects is, however, smooth in contrast to the results obtained at Pr = 6.0 and Pr = 1.0. The smooth transition is very evident in Figs. 18 and 19 which represent the a and y vs R at which Gr ocure crit occurred at Pr = 0.025, unlike the corresponding curves at Pr = 6.0

90 220 - 210 - 200 - Pr = 0.025 190 - 180 - 170 160 - 150 - 140 - 7 130 - x 120o0 100 LL\ 9080T 0 70 60 50 40 302010 0 20 40 60 80 100 120 R = RT Figure 17. Critical Grashof No., Grcrit, as Function of R =RT-1/2. (Pr = 0.025).

91 4.0 Pr = 0.025 3.0 - 2.0 ad 1.00.0 - ~ I, I' I...~ I' I I 0 20 40 60 80 100 120 R = RT 1/2 Figure 18. Wave No., a, (at Grcrit) as function of R= RT1/2. (Pr = 0.025).

92 0 -10 -20 -30 -40 -50 - Pr = 0.025 -60 - 70 - s-80 a90 0 * 80,j 70 0 60 50403020 10. - v...I.,,,. 0 20 40 60 80 100 120 R =RT-1/2 Figure 19. Orientation, y, (at Grcrit) as Functionof R = RT-1/2 (Pr = 0.025).

93 340 320 4.0 Pr = 0.025 - 300 3.8- - 280 3.6- - 260 3.4- 240 3,2- 220 3.0- 200 2.8- 180 2.6- 160 _2.4 / 140 2.2- 120 2.0- 100 1.8 80 1.6- 60 1.4, 40 1.2 20 1.0 l 0 0 20 40 60 80 100 R = RT-1/2 Figure 20. Wave Speed, cr/R, and Frequency X, at Grcrit as Functions of Re = RT /2. (Pr - 0.025, T - 50),

94 or 1.0 (Figs. 12, 13, 15, and 16). This smooth transition also occurred in the wave speed, c /R, or frequency, Xi Figure 20 presents these results. X is the imaginary part of the eigenvalue (c.f., Eq. (32)). The correspondence between Xi and c /R is cr/R = A /aR. Referring to the i curve of Fig. 20, it is quite evident that X varies smoothly as a function of R starting at the value of 44.5 at Re = 0, which is the value given by Chandrasekhar (1961) for the rotating Benard problem at T = 50 (Ta = 104), and increasing with increasing R to a final value of i = 333 at R = 109.0 ~ E which is the value obtained for the unstratified flow at T 50. Referring to the c r/R curve of Fig. 20, it is interesting to note that cr/R first decreases rapidly, goes through a minimum at R = 30, then increases going through an inflection point at R5 - 45 and continues to increase smoothly to a final value of 0.236 at R5 = 109.0 which is the c /R of the "A" wave of the unstratified flow at T = 50 (c.f., Table V). Referring now to Fig. 17, Grcrit is observed to increase with R5 going through an inflection point at R5 - 30, reaching a maximum at R 45 and then decreasing smoothly to Gr = 0 at R 109.0. The minimum of the c r/R curve occurs at crit e about the same Rc as the inflection point in the Grcrit curve and the inflection point in the c r/R curve occurs at about the same Rc as the maximum of the Grcit vs RE curve. 4.2.1.3 Prandtl Number Effect The discussion above has already pointed out that one of the Pr effects was that the critical disturbance was always nonstationary for Pr = 0.025 and T = 50. For this result to be consistent with the results of the rotating Benard problem (which it must be at RE = 0),

95 the critical disturbance of the present problem should be nonstationary for all Pr < 0.677 and T > 11 (where the criteria, T'> 11, has been taken from the free boundary case of the rotating Benard problem). Figure 21 is a comparison of the results obtained from the three Pr (6.0, 1.0, and 0.025) at T = 50. To reduce the results to a common bases, the axes chosen are ~crit vs R. It is clearly evident that increasing Pr has a stabilizing effect on both the stationary and the nonstationary modes, except at R - 0 where Pr has no effect on the stationary modes and at R = 0 where, as expected, Pr has no effect on the hydrodynamic or unstratified flow. It is of interest to note that the stabilizing effect of Pr on stationary modes (Pr = 1.0 and 6.0) was, in a sense, observed by Gallagher and Mercer (1965) in their study of the stratified pCf. That is, GM found that the Prandtl number had a strong stabilizing effect on disturbances with a transverse component. This effect appears to be quite strong in the present problem also. This concludes the discussion of the critical curves. 4.2.2 Neutral Curves and Eigenvalue Spectra As discussed in paragraph 3.3, the critical curve information was generated on the computer without generating detailed neutral curve results. However, even in this mode of operation, the neutral curve behavior could at least be partially recognized. The behavior of the neutral curves which seemed to be of interest is discussed below. The neutral curves considered below are all for Pr = 1.0. Comparing the results at Rc - 10, Fig. 22, and R = 40, Fig. 24, there is a very apparent effect of R. on both the shape of the Gr - Y neutral curve and the range of Gr. The range of Gr at

(OS = A)'ad''ON IlpueT jo oajj~g 2urzI-IqglS atUT'IZ arnold 3 z/L - = OZL 00 L 08 09 O0 OZ O'C 0 09 09 96

97 R = 40 is about three times as large as the range of Gr at R 10. At RE = 10, Fig. 22, the disturbance at neutral stability is stationary throughout the range of y. At R - 40, Fig. 24, the neutral disturbances are stationary in the range -54~ < y < 840 and outside of that range of y they are nonstationary with the stationary disturbances the most critical. This type of behavior was noticed at all of the values of T studied when R became sufficiently large. On further increasing Rc (R > 40 at T = 100), the range of Y, which resulted in stationary modes, would be reduced. Conversely, the two nonstationary loops shown on the extremeties of the y range in Fig. 24 would, on increasing R, encompass a larger range of Y. However, it was not this behavior which resulted in the "A" or "B" waves. The "A" and "B" waves were caused by a third nonstationary loop appearing in the Gr - Y neutral curves when RE became sufficiently large. At T = 100 (and Pr = 1.0), this loop occurred in the vicinity of y = -200. This additional loop is noet shown in any of the figures included in this thesis. The "A" or "B" wave loop was, in general, very distinct and therefore, quickly became the most critical portion of the Gr - Y neutral curve on slight increases of R. This is, of course, implied in the critical curves for Pr = 1.0 and T = 100 (e.g., Figs. 14 and 16) which show abrupt changes in both Gr it and the y at which Grcrit occurs when the "A" wave appears. Considering Figs. 23 and 25, it is evident that there is very little qualitative effect of R on the shape of the Gr - neutral curve at T = 100.

98 At lower Taylor numbers (T < 70) the Gr - a neutral curves were observed to have double noses in certain ranges of R. Figure 26 depicts such behavior at T - 20, R 67.1, and y = -50.0~. The E y chosen is close to the y at which Grcrit occurred for Pr = 1.0, T = 20, and R = 67.1 (y it= -58.0~, c.f., Table VII). This double nose structure was observed to occur in that region of the a - R space where a increases with increasing R. In Fig. 15, for T = 20 (Pr = 1.0), a increases with R for R 5 40. Below R - 40 ~ e ~ at T = 20, the double nose structure was not observed. As R increased past R = 40, the double nose structure became more evident. At T = 35 and T = 50, the double nose structure was evident but it was not as distinct as at T = 20. At T = 70 and T = 100, the double nose structure was not observed. The Gr - a neutral curves at T = 70 and T = 100 appeared to always be similar to the two shown in Figs. 23 and 25. It is to be emphasized that a definitive study of the neutral curve behavior was not made and, consequently, any of the observations stated above which are not shown by Figs. 22 through 26 must be considered only as probable typical behavior. In most cases when it was suspected that a critical point had been found during a preliminary search (which was guided by the limiting cases) in conversational mode with the computer, a thorough search of the a and y planes was made, using the eigenvalue program, to determine that no other mode could be found which would be more unstable. That is, the Gr ri was taken as constant input to the eigenvalue program and if, in a search of the a and y planes, the growth rate

99 was never positive, it was assumed that the critical value had indeed been found. This mode of operation worked quite well for it was through this method that the a's were first observed to decrease and then increase on increasing R (c.f., Fig. 15), and it was in this manner that the "A" and the "B" waves were first observed for the stratified flow. One particular characteristic of the eigenvalues was noticed which appears to be of interest. Figure 27 shows the lowest eigenvalue, X, in the spectrum, as a function of Gr, for T = 100, RE = 40, Pr = 1.0, a = 5.56, and y = 80~. It is noted that the 3 3 eigenvalue splits at Gr = 30.8x10. That is, at Gr < 30.8 X 10 X represents, in fact, the two lowest eigenvalues which are a complex conjugate pair and represents a nonstationary mode. On increasing Gr, this? splits into two stationary modes before becoming unstable, with one branch decreasing rapidly, the other increasing. If this splitting were to occur after the sign change, i.e., after X becomes positive, then the neutral disturbance was nonstationary. It was this behavior which was associated with the nonstationary regions of the neutral curve shown in Fig. 24. This was not, as far as could be determined, the same eigenvalue which was associated with the "A" and "B" waves. This completes the discussion of the results and conclusions concerning the stratified flow.

100 18 - T = 100 R = 10 16 ~ a= 5.42'o 14 8 12 10 -80 -60 -40 -20 0 20 40 60 80 Orientation, y Figure 22. Neutral Curve, Gr vs y, (Pr =1.0, T= 100 R= 10 O = 5.42). 38 - T = 100 R34 -= 10 34 y = 75.Q~ 302622 1816 12 0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 Wave No., a Figure 23. Neutral Curve, Gr vs a, (Pr = 1.0, T = 100, RE = 10, y = 75.0~).

101 50 T = 100 R= 40 a = 5.56 40 Co 30 20 stationary J non-stationary 10 I I -80 -60 -40 -20 0 20 40 60 80 Orientation, y Figure 24. Neutral Curve, Gr vs y (Pr = 1.0, T = 100, R = 40, a- = 5.56). 90 T = 100 80 - R = 40 70 - Y = 37.0~ 60 - L0 5- 5040 3020 - 10 0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 Wave No., a Figure 25. Neutral Curve, Gr vs a, (Pr = 1.0, T = 100, R = 40, y =37.0o).

102 25 2423 22 21 T = 20 20 R = 67.1 x \ I r = -50.0~ 19_ I/ 18 1716, 0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 Wave No., a Figure 26. Neutral Curve, Gr vs a, (Pr= 1.0, T= 20, R: 67.1, y= -50.0~). -20 -15 - *-10 - stable region -5 unstable 5 region 10 25 30 35 GRASHOF No., Gr x 103 Figure 27. Growth Rate as Function of Grashof No. (Pr = 1.0, T = 100, RE = 40, y = 80~, a = 5.56).

Appendix I STABILITY EQUATIONS OF THE STRATIFIED FLOW IN THE LOCAL COORDINATE SYSTEM The invariant form of the governing equations is given below. V.V= 0 V vvv +V t iv v V AX + C +V-7V = C V'G 6t These equations must be represented in the skewed local coordinate system of Fig. 4 where do= h&j, dL c at'.- ro - O ~- _ Define V +bv+ w where A -YLt +Xi 103

104 Equations (I-1) become, after eliminating the pressure between the. and k component, momentum equations LA 3u a~~~~~~~ aw'q~ (I-2) + I V V+} a DU I V 3l V /~ ) _ o Dt ax ra'a ql; X } D.a L v u 1 _ a _AO~qa, h 2. ose where hL at 8Vy 3L.1 a -A- The perturbation equations are obtained by assuming disturbances of the form given below, where dimensional variables are now represented by (^)

105 _^ _A iA _ A Q t) (M13) = VYj3) + -' V (rL45), pP - ea) t' p (,)ld) e. A A,, ^, A, ^ I A \ A A -C and nondimensionalizing with [ktA]- [hoj',o (1-4) A L. | to ay' + bj + t (I-5) The perturbation equations become j(+v+ W~Y -bC7) \ L L -l-3\d tt - h rt(Ve8 r tL 6, L LL (4 3 - L t a +v +VW V~b3 ro ro COr [?tybl >3 3'd3 -l( + L tL f A O\) Y ~ Av L (L \A / 32. TqA.1 A LV

106 U)L{ Du /.-Z) O V -V- t\N UV 4- qV) r2 UoL Pr| A + Pr w, Vu + - = D 2 where D and V are now defined as U L u 0 - t ( -+t L h i } a~a' -3L L + LIIII and use was made of o,-, 3 -, where w = O results by considering the conservation of mass in the local coordinate system, i.e., Li ( dth- Vw ) & =o (I-6)

107 thus correct to 0(t0 r~=0 and to satisfy the boundary conditions on w w O Dropping terms of 0(r-) in the perturbation equations and defining the stream function as w = -ia~ (I-7) v= ~_ ay and recalling the base flow solution (transformed to local coordinates), i.e., u = hai v= hv 0 = - -(1+ 2z) the perturbation equations become ) X\L/;L7 ( t(/}2 = IIVdI(I-8) olnV, (o-[V )V +Z — -3OIQ = _o _i( V- - U = t~~~oL~Q Prt( v UF)(3+t~ t,

108 where d (perturbation quantities) D =ddz 2 2 2 V = D - d (base flow quantities) Primes denote d To complete the separation of variables, recall that the approach has been to describe the flow in the vicinity of r. In this vicinity, h 0(1); therefore, define 0 U Lh U 0L R = ---- a local Reynolds number (I-9) 2 2 2r0h 2r0 Fr= 0 - a local Froude number g g also define T - V Cr = acTgL Gr =-7T2 2 The local perturbation equations become (I-10) V,-t - T o RiCV cdVIX -13 1 zT ou -Lo( Gr B + GrEyr D = O, - /- ~( Rv -:/R - \ + LT D~ -0 y1,'e - ci, P,-r -, / - c Pr 0 = O

109 which must satisfy 1 = D = u = 8 = 0 at z = _ where i and v are the z-dependent similarity base flow velocity profiles transformed to the local skewed coordinate system and are given by u = os y + v sin y (I-ll) v = -u1 sin y + vl cos y with system at r0 V, AIr'kC U) "JT C'k JT Since all of the curvilinear terms have been eliminated, Eqs. (I-10) are equivalent to assuming a local cartesian coordinate system at r0. The unstratified local perturbation equations result by setting 0 - O.

Appendix II FINITE DIFFERENCE SOLUTION TO THE EIGENVALUE PROBLEM OF THE STRATIFIED FLOW The finite difference solution technique (discussed in paragraph 3.1.1) applied to the stability equations for the stratified flow regime, Eqs. (30), is developed below in detail. This development also applies to the unstratified flow regime, Eqs. (16), by setting e _ 0. Equations (30) are written below after regrouping the real and imaginary operators separately. = (cIf) $ z TrD + t4-E.-+ i 4{RU ~ - RV - = q (11-2) tD e - - L { + ~ v - 3 R (II-3) (a )- xt o at _ = (11-4) where D. d dz Primes, ('), denote d (base flow quantities) dz Dots, (.), denote d — -isc 110

111 Equations (II-1) through (II-3) are approximated with finite differences by dividing the z domain into N intervals (N + 1 nodes) and expanding in Taylor series expansions around the interior nodes. Equations (II-1) and (II-3) are expanded around N - 1 nodal points, Eq. (II-2) around N half nodal points. The variables 4 and 0 are defined only at nodal points and u is defined only at half nodes. Choosing centered differences, define, -(j )'/3. (II-5) A 1 = where N and j denotes a nodal point. The same form of the above operators holds at the half nodes. Using the above difference representation, the Taylor series expansions take the form, after rearranging,

112 33 2 ~d A * z 6 - —'i * o+0(a) e+ - -- s+ 4-j- a +) - _ A The last two expressions in (II-6) are used to satisfy the conditions 4 and u be defined only at nodes and half nodes, respectively. Difficulty with this condition is met only for the ~ in Eq. (II-2), expanded around half nodes, and for the boundary conditions D(~+1/2) = u(+l/2) = 0. Equations (II-6) are substituted into Eqs. (II-1) through (II-3), including the terms of O(A 2) in all expressions, except the highest order derivatives in each of Eqs. (II-1) through (II-3). 2 The terms of O(A 2) will be referred to as the remainder terms. The remainder terms are included, where possible, to minimize the truncation error. The remainder terms are not included in the representation of the highest order derivatives because the finite difference equation would then be of higher order than the differential equation it is to represent. The approximate solution could not then be expected to behave in the same manner as the exact solution. Substituting (II-6) into Eqs. (II-1) through (II-3) and collecting like variables at indicated nodes, j, or half nodes, i,

113 the real LHS operators give -rc~v- Ec. I-14) (II-7) A 1(.1) + - 4(v u,.2) + A(.+A 1(5) A. A A2.() (Y- AZ(-2,: 4- AZ(3) bi,: + A Z(4) t),3 + AZ (5)qL_, + ())U; +AN(7))L,c1 frorAl (O]3) A3(1i e_ +A3(.' OG A3(3) G, Similarly, the imaginary LHS operators give +ron_. E (n:) BlUWS9j f _1(2)j 4B, B1(3i +BBl] 1(5)Ij q), j 4+ B(6) -l Brom E\ (_B-_.rom ~%..L t.-s) 8t(l)ul ti~ tt 3s L(2)u~., e U (3)U,~O., B;,, et Lt; tr +132(5)V; 4romrc iR (I-3) and the RHS operators (real) give troh ~q....... C1(1) Lw + 10(z)' C(3) fj 4) + cltits) uL where

114 (II-8) Al(1) = ( i + a Bi(1) = aR ~6 ~A2~ ~12A2 A14(2) + 2 R ( 3 Ai(2) = - 72 +' Al(3) =- + 2 (41) + 1(3) =+ = R 2 1) A 4 A2 A A1(4) = A1(2) BB1 = aR Al(5) = Al(1) B1(4) = B1(2) A1(6) = -12 B1(5) = B1(1) A1(7) = ( + 4 B(6) = -aGr A1(8) = -A1(7) A1(9) = -A1(6) A2 (1) = B2(1) = A2 2(2) = - R + B2 (2) = R A2(3) = -A2(2) B2(3) = B2(2) A2(4) = -A2(1) B2(4) = B2(1) A2(5) = 1 B2(5) = -oR A2 A2(6) - + 2) A2(7) = A2(5) A3(1) =2 Pr B3(1) = -a A3(~- Pr B3(2) = -aR A3(3) = A3(1)

115 Cl(1) - - 1 12A2 1+ 1/3 C1(2) - C1(3) = - +2+ cl(4) = C1(2) C1(5) = Cl(1) where the underlined terms are the remainder contributions. The boundary conditions (Eqs. (II-4)) become t -(.-o.U (II-9 L'4- Li -, = W ~,~ -gLoI= N-/ =3 E) = 0 On indexing j and i, in (11-7), as j = 1, 2,... N - 1 1 3 2N - 1 2' 2'''' 2 satisfying the boundary conditions (II-9) and evaluating the base flow quantities j.' i' and u'i at the nodes (or half nodes) thereby defined, there are obtained 3N -2 equations which represent an algebraic eigenvalue problem in 3N - 2 unknown eigenvalues. In general, for stability problems at the neutral point, only the lowest eigenvalue is of interest.

116 The equations necessary to evaluate the base flow velocity profiles and derivatives are — u, Amk cv,~c. U U, 2~. Z~A A + B (A+ + B V/ _ - u l ALAtE t \JI LQA a, - -io~h$A +T-t, V, - Lt?_FT_ d ~ 44lw >wnU

117 where, in the nodal notation, the argument z is evaluated as z = k-iN k = 1, 2,.. N (II-ll) 2N in the positive z domain, where k-odd is associated with half nodes and k-even nodes, when N is odd, the symmetry or the antisymmetry of (II-10) being used to obtain the values in the negative z domain. The eigenvalue problem defined by (II-1) through (II-4) can be represented in matrix form as [A + iB] u= [C] u on premultiplying by C-1 [C [A + iB] - XI] u = (II-12) where I = unit matrix X = -ia(cr + ici), the complex eigenvalue and the elements of A, B, and C are defined by (II-7) through (II-11). If u# O0, then the determinant C-1[A+ iB] - RI O The method actually used to find the eigenvalue is discussed in paragraph 3.1.2.

BIBLIOGRAPHY Barcilon, V. and Pedlosky, J., "Linear Theory of Rotating Stratified Fluid Motions." J. Fluid Mech., 29 (1967a), 1-16. Barcilon, V. and Pedlosky, J., "A Unified Linear Theory of Homogeneous and Stratified Rotating Fluids." J. Fluid Mech., 29 (1967b), 609-621. Barcilon, V. and Pedlosky, J., "On the Steady Motions Produced by a Stable Stratification in a Rapidly Rotating Fluid." J. Fluid Mech., 29 (1967c), 673-690. Batchelor, G. K., "Note on a Class of Solutions of the Navier-Stokes Equations Representing Steady Rotationally-Symmetric Flow." Quart. J. Mech. and Appl. Math., 4 (1951), 29-41. Bodewadt, U. T., "Die Drehstromung iiber festem Grunde." Z. angew. Math. Mech., 20 (1940), 241-253. Brunt, D., "Experimental Cloud Formation." Compendium of Meteorology American Meteorological Society, (1951), 1255-1262. Chandra, K., "Instability of Fluids Heated from Below." Proc. Roy. Soc., A 164 (1938), 231-242. Chandrasekhar, S., "The Instability of a Layer of Fluid Heated Below and Subject to Coriolis Forces." Proc. Roy. Soc., A 217 (1953), 306-3244 Chandrasekhar, S., Hydrodynamic and Hydromagnetic Stability. Oxford, 1961. Chandrasekhar, S. and Elbert, D. D., "The. Instability of a Layer of Fluid Heated Below and Subject to Coriolis Forces. II." Proc. Roy. Soc., A 231 (1955), 198-210. Deardorff, J. W., "On the Stability of Viscous Plane Couette Flow." J. Fluid Mech., 15 (1963), 623-631. Deardorff, J. W., "Gravitational Instability Between Horizontal Plates with Shear." Phys. Fluids, 8 (1965), 1027-1030. Dropkin, D. and Globe, S., "Effect of Spin on Natural Convection in Mercury Heated from Below." J. Appl. Phys., 30 (1959), 84-89. Duncan, I. B., "Axisymmetric Convection Between Two Rotating Disks." J. Fluid Mech., 24 (1966), 417-449. Ekman, V. W., "On the Influence of the Earth's Rotation on Ocean Currents." Ark. Mat. Astr. Fys., 2 (1905), 1-52. 118

119 Faller, A. J., "An Experimental Study of the Instability of the Laminar Ekman Boundary Layer." J. Fluid Mech., 15 (1963), 560-576. Faller, A. J. and Kaylor, R. E., "A Numerical Study of the Laminar Ekman Boundary Layer." J. Atmos. Sci., 23 (1966), 466-480. Forsythe, G. and Moler, C. B., Computer Solution of Linear Algebraic Systems. Prentice-Hall, 1967. Fultz, D., Nakagawa, Y. and Frenzen, P., "An Instance in Thermal Convection of Eddington's "Overstability"." Physical Rev., 94 (1954), 1471-1472. Fultz, D., Nakagawa, Y. and Frenzen, P., "Experiments on Overstable Thermal Convection in Mercury." Proc. Roy. Soc., A 231 (1955), 211-225. Gallagher, A. P. and Mercer, A. McD., "On the Behaviour of Small Disturbances in Plane Couette Flow." J. Fluid Mech., 13 (1962), 91-100. Gallagher, A. P. and Mercer, A. McD., "On the Behaviour of Small Disturbances in Plane Couette Flow with a Temperature Gradient." Proc. Roy. Soc., A 286 (1965), 117-128. Goroff, I. R., "An Experiment on Heat Transfer by Overstable and Ordinary Convection." Proc. Roy. Soc., A 254 (1960), 537-541. Greenspan, H. P., The Theory of Rotating Fluids. Cambridge, 1969. Gregory, N., Stuart, J. T., and Walker, W. S., "On the Stability of Three-Dimensional Boundary Layers with Application to the Flow Due to a Rotating Disk." Phil. Trans. Roy. Soc., A 248 (1955), 155-199. Grohne, D., "i'ber das Spektrum bei Eigenschwingungen. ebener Laminarstromungen." Z. angew. Math. Mech., 35 (1954), 344-357. Homsy, G. M. and Hudson, J. L., "Centrifugally Driven Convection in a Rotating Cylinder." J. Fluid Mech., 35 (1969), 33-52. Homsy, G. M. and Hudson, J. L., "The Asymptotic Stability of a Bounded Rotating Fluid Heated from Below: Conductive Basic State." J. Fluid Mech., 45 (1971a), 353-373. Homsy, G. M. and Hudson, J. L., "Centrifugal Convection and its Effect on the Asymptotic Stability of a Bounded Rotating Fluid Heated from Below." J. Fluid Mech., 48 (1971b), 605-624. Hopf, L., "Der Verlauf kleiner Schwingungen auf einer Stramung reibender Flussigkeit." Ann. Phys., Lpz., 44 (1914), 1-60.

120 Karman, T. von, "Laminare und Turbulente reibung." Z. Angew. Math. Mech., 1 (1921), 233-251. Korpela, S. A., "Stability of the Conduction Regime of Natural Convection 1. Flow in an Inclined Slot 2. Flow in a Narrow Rotating Annulus." Ph.D. Dissertation, The University of Michigan, Dept. of Mech. Eng., 1972. Koschmeider, E. L., "On Convection on a Uniformly Heated Rotating Plane." Beitr. Phys. Atmos., 40 (1967), 216-225. Lance, G. N. and Rogers, M. H., "The Axially Symmetric Flow of a Viscous Fluid Between Two Infinite Rotating Disks." Proc. Roy. Soc., A 266 (1962), 109-121. Lin, C. C., The Theory of Hydrodynamic Stability. Cambridge, 1966. Lilly, D. K., "On the Instability of Ekman Boundary Flow." J. Atmos. Sci., 23 (1966), 481-494. Martin, R. S. and Wilkinson, J. H., "Similarity Reduction of General Matrices to Hessenberg Form." Num. Math., 12 (1968a), 349-368. Martin, R. S. and Wilkinson, J. H., "The Modified LR Algorithm for Complex Hessenberg Matrices." Num. Math., 12 (1968b), 369-376. Mellor, G. L., Chapple, P. J. and Stokes, V. K., "On the Flow Between a Rotating and a Stationary Disk." J. Fluid Mech., 31 (1968), 95-112. Millsaps, K. and Pohlhausen, K., "Heat Transfer by Laminar Flow from a Rotating Plate." J. Aero. Sci., 19 (1952), 120-126. Moler, C. B., "Matrix Computations with Fortran and Paging." Standford Comp. Sci. Rept., (STAN-CS-71-196), (1971). Nakagawa, Y. and Frenzen, P., "A Theoretical and Experimental Study of Cellular Convection in Rotating Fluids." Tellus, 7 (1955), 1-21. Niler, P. P. and Bisshopp, F. E., "On the Influence of Coriolis Force on Onset of Thermal Convection." J. Fluid Mech., 22 (1965), 753-761. Parlett, B. N. and Reinsch, C., "Balancing a Matrix for Calculation of Eignevalues and Eigenvectors." Numer. Math., 13 (1969), 293-304. Pearson, C. E., "Numerical Solutions for the Time-Dependent Viscous Flow Between Two Rotating Coaxial Disks." J. Fluid Mech., 21 (1965), 623-633.

121 Picha, K. G. and Eckert, E. R. G., "Study of the Air Flow Between Coaxial Disks Rotating with Arbitrary Velocities in an Open or Enclosed Space." Proc. 3rd U.S. Nat'l. Cong. on Appl. Mech., (1958), 791-798. Reid, W. H. and Harris, D. L., "Some Further Results on the Benard Problem." Phys. Fluids, 1 (1958), 102-110. Riley, N., "The Heat Transfer from a Rotating Disk." Quart. J. Mech. and Appl. Math., 17 (1964), 331-349. Rogers, M. H. and Lance, G. N., "The Boundary Layer on a Disc of Finite Radius in a Rotating Fluid." Quart. J. Mech. and Appl. Math., 17 (1964), 319-330. Rossby, H. T., "A Study of Benard Convection with and without Rotation." J. Fluid Mech., 36 (1969), 309-335. Schultz-Grunow, F.,''Der Reibungswiderstrand Rotierenden Scheiben in Gehausen." Z. angew. Math. Mech., 15 (1935), 191-204. Stewartson, K., "On the Flow Between Two Rotating Coaxial Disks." Proc. Camb. Phil. Soc., 49 (1953), 333-341. Squire, H. B., "On the Stability for Three-Dimensional Disturbances of Viscous Fluid Flow Between Parallel Walls." Proc. Roy. Soc., A 142 (1933), 621-628. Veronis, G., "Cellular Convection with a Finite Amplitude in a Rotating Fluid." J. Fluid Mech., 5 (1959), 401-435. Veronis, G., "Motions at Subcritical Values of Rayleigh Number in a Rotating Fluid." J. Fluid Mech., 24 (1966), 545-554. Wasow, W., "One Small Disturbance of Plane Couette Flow." J. Res. Nat. Bur. Stand., 51 (1953), 195-202. Wilkinson, J. H., The Algebraic Eignevalue Problem. Oxford, 1965. Zeipel, H. V., "The Radiative Equilibrium of a Rotating System of Gaseous Masses." Monthly Notices of Roy. Ast. Soc., 84 (1924), 665, 684, 702. Zondek, B. and Thomas L. H., "Stability of Limiting Case of Plane Couette Flow." Phys. Rev., 90 (1953), 738-743.

UNIVERSITY OF MICHIGAN 9II0I JJIlllJ11111111115020 I86 I8710111 3 9015 02086 8710