THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING STABILITY OF A FLUID IN A LONG HORIZONTAL RECTANGULAR CYLINDER HEATED FROM BELOW Michael R. Samuels A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan Department of Chemical and Metallurgical Engineering 1966 January, 1966 IP-728

Doctoral Committee: Professor Stuart W. Churchill, Chairman Professor Vedat S. Arpaci Assistant Professor Joe D. Goddard Assistant Professor James 0. Wilkes Professor Chia-Shun Yih

ACKNOWLEDGEMENTS The author wishes to acknowledge the assistance of the members of his doctoral committee. Special thanks are due to Professor S. W. Churchill, Chairman of the Committee. The author also wishes to thank the National Science Foundation, whose financial assistance permitted the author to devote full time to this work. iii

TABLE OF CONTENTS Page ACKNOWLEDGEMENTS................................................... iii LIST OF TABLES.................. vi LIST OF FIGURES.................................................... vii NOMENCLATURE................................................... ix ABSTRACT............................................................. xiii I INTRODUCTION................................................. 1 II HISTORICAL BACKGROUND AND LITERATURE SURVEY.................. 3 2.1. Natural Convection Between Flat Plates of Infinite Extent............................................... 3 2.2. Natural Convection in Enclosed Finite Regions.......... 9 2.3. Other Numerical Solutions.............................. 12 III DEVELOPMENT AND DESCRIPTION OF THE MATHEMATICAL MODEL........ 13 3.1. Equations Governing the Problem........................ 16 3.2. Dimensionless Form of the Equations................... 19 3.3. Evaluation of UC....................................... 22 3.4. Elimination of Pressure from the Equations of Motion and Introduction of Stream Functions and Vorticities. 23 3.5. General Logic of-the Numerical Solution................ 26 IV FORMULATION OF FINITE DIFFERENCE APPROXIMATIONS AND THEIR STABILITY................................................... 27 4.1. Formulation of Finite Difference Approximations........ 27 4.2. Stability Considerations............................... 33 4.3. Definition of Stability................................ 34 4.4. Determination of the Stability of a Finite Different Approximation...................................... 35 4.5. The Method of Positive-Type Equations.................. 36 4.6. The von Neumann Stability Criterion.................... 38 4.7. Extension of Finite Difference Approximations to the Equations of Motion and Energy, and Their Stability.. 42 4.8. Boundary Conditions and Stability..................... 50 4.9. Advantages and Disadvantages of Explicit, Implicit, and Alternating Direction Implicit, Finite Difference Formulations......................................... 51 iv

TABLE OF CONTENTS (CONT'D) Page 4.10. Boundary Conditions.......................... 52 4.11. Determination of *n+l when ~n+l are Known....... 55 4.12. Evaluation of the Coefficients, U, V, and Nu..... 58 4.13. Sequence of Operations for Integration of the Vorticity and Energy Equations.................... 60 V EXPERIMENTAL WORK.......................................... 62 5.1. Apparatus.......................................... 62 5.2. Discussion of Experimental Observation.............. 68 5.3. Conclusions and Recommendations..................... 71 VI RESULTS.............................................. 73 6.1. Introduction..................... 73 6.2. Discussion of Results................ 76 a. Grid Size and Numerical Stability............... 76 b. Nusselt Number Variations with L/H, Gr, and Pr. 79 c. Determination of the Critical Ra, and Comparison with the Theoretical Solution of Kurzweg, and the Numerical Solution of Aziz............ 84 d. Streamline and Isotherm Plots.......... 88 e. Preferred Mode of Natural Convection............ 103 BIBLIOGRAPHY............................................... 105 APPENDIX A SOLUTION OF A SYSTEM OF LINEAR EQUATIONS HAVING A TRIDIAGONAL COEFFICIENT MATRIX................... 109 APPENDIX B DESCRIPTION OF THE PROGRAMS AND VARIABLES USED IN THE NUMERICAL STUDY................................. 111 APPENDIX C TYPICAL COMPUTER OUTPUT............................. 126

LIST OF TABLES Table Page I Summary of Cases Studied............................... 74 II Effect of Grid Size on Calculations in a Square Cavity for Gr = 3000., and Pr = 1.0.................. 76 III Central r and Nu as a Function of Pr at Ra = 10000., and L/H = 1.0............................ 83 IV Critical Ra as a Function of L/H and Pr......... 87 vi

LIST OF FIGURES Figure Page 1 The Physical Problem..................................... 15 2 The Numerical Grid........................................ 27 3 Front View of the Experimental Chamber.................... 63 4 The Cooling Plates...................................... 64 5 Schematic Diagram of the Heater Power Supply.............. 66 6 The Heating Plates........................................ 67 7 Sketch of Flow Patterns Seen in Liquid Glycerine at Ra 3000................................................... 69 8 Effect of Grid Size on Steady State Nu and *cent for a Square Cavity with Pr = 1, Gr = 3000.............. 78 9 Dependence of Nu on Ra and Pr for a Cavity of Square Cross Section..................................... 80 10 Comparison of Results with the Calculations of Aziz and the Data of Schmidt and Silveston..................... 82 11 Dependence of Nu on Gr and Pr for a Cavity of Square Cross Section...................................... 85 12 Effect of L/H and Ra on Nu for Pr = 1.0............ 86 13 Dependence of Strength of Convection on (Ra-Rac) for L/H = 1 0 and Pr = 1.0.................................. 89 14 Steady State Streamlines and Isotherms for L/H = 1.0, Pr = 1.0, and Gr = 1400................................. 90 15 Steady State Streamlines and Isotherms for L/H = 1.0, Pr = 1.0, and Gr = 2000................................. 91 16 Steady State Streamlines and Isotherms for L/H = 1.0, Pr = 1.0, and Gr = 5000................................. 92 17 Steady State Streamlines and Isotherms for L/H = 1.0, Pr = 1,0, and Gr = 10000................................ 9 18 Steady State Streamlines and Isotherms for L/H = 1.0, Pr = 1.0, and Gr = 20000................................ 94 vii

LIST OF FIGURES CONT'D Figure Page 19 Steady State Streamlines and Isotherms for L/H = 0.5, Pr = lo0, and Gr = 6000................................. 96 20 Steady State Streamlines and Isotherms for L/H = 0.5, Pr = 1, and Gr = 8000................................ 97 21 Steady State Streamlines and Isotherms for L/H = 0.5, Pr = 1.0, and Gr = 10000................................ 98 22 Steady State Streamlines and Isotherms for L/H = 0.5, Pr = 1.0, and Gr = 20000................................ 99 23 Steady State Streamlines and Isotherms for L/H = 2.0 Pr = 1, and Gr = 10000................................. 101 24 Steady State Streamlines and Isotherms for L/H = 3.0, Pr = 1, and Gr = 3000.................................. 102 viii

NOMENCLATURE Note: Certain symbols which have very limited scope and whose definition is clearly indicated in the text will not be mentioned here. Symbol Usage [A] Amplification matrix which describes growth rate of the a, P component of T and: a,b, c,d Components of the amplification matrix [A] cp Heat capacity at constant pressure, BTU/(lb.)(~F) f(x,y) Initial velocity conditions for integration of momentum equation F(X,Y) Initial conditions for integration of the dimensionless momentum or vorticity equation g Acceleration of gravity Gr Grashof number based on 1/2 temperature difference Gr = 1/2 Hg2 v2 h Local heat transfer coefficient, BTU/(hr)(ft2) (~F) H Cell height such that 0 < y < H k Thermal conductivity, BTU/(hr) (ft) (~F) K Thermal diffusivity, K = k/cpp L Cell length such that 0 < x < L Nu Nusselt number,. Nu = QH/LkAQ NX Number of vertical grid spaces NY Number of horizontal grid spaces p Pressure P Dimensionless pressure, P = pH2/poVo Pr Prandtl number, Pr = cp=/k q Heat flux density, BTU/(hr)(ft2 ) Q Total heat transferred/Unit of cavity depth, BTUJ/(hr)(ft) ix

Symbol Usage r AX2/AY2 Ra Rayleigh number based on 1/2 temperature differences = Gr x Pr =i Cp H3g-A 2 k v Rac Critical Rayleigh number Re Reynolds number t Time T Dimensionless temperature u x direction velocity U Dimenisionless X direction velocity, U = - v vy direction velocity vH _ V Dimei' sionless Y direction velocity, V = - x Horizontal coordinate x x/H y Vertical coordinate Y y/H Z a, component in the Fourier expansion for GREEK LETTERS ap X and Y component frequencies in a double Fourier expansion AX X direction grid spacing = L/NX/H AY Y direction grid spacing = 1/NY aL9 Temperature difference between upper and lower surfaces of cavity AT Dimersionless time step size 62 a2 V2 Laplacian operator, - + 2 5FDA to - x

Symbol Usage 2 b2 5X FDA to _ Dimensionless vorticity = - V2 G Temperature gcp ap component in the Fourier expansion for T xlX,2 Eigenvalues of the amplification matrix [A] General dependent variable —function of X, Y, and T v Kinematic viscosity ap Amplification factor of the a,p component of f p Density T Dimensionless time = tv/H2 a, Q ac, component of the Fourier expansion for ( Dimensionless stream function SUBSCRIPTS cent central h High i ith grid point.th ~~~j j grid point 1 Low o Evaluated at the mean temperature (Gh + Q1)/2 SUPERSCRIPTS h nth time step Variation from conditions at mean temperature xi

Symbol Usage R* Intermediate value Mean value ABBREVIATIONS ADI Alternating direction implicit FDA Finite difference approximation FDE Finite difference equation PDE Partial differential equation xii

STABILITY OF A FLUID IN A LONG HORIZONTAL RECTANGULAR CYLINDER HEATED FROM BELOW Michael Robert Samuels ABSTRACT If a fluid, initially at rest within an enclosure, is heated from below and cooled from above, it will remain at rest as long as it is not disturbed. If the fluid is sufficiently disturbed, it will either return to its quiescent state or undergo a complete transition to a convective regime, depending on the temperature difference, AQ, across the fluid, the physical properties of the fluid, and the geometry of the container. The critical temperature difference, A9c, is defined as that temperature difference above which convection will appear, and below which the quiescent fluid is stable. The objective of this study was to demonstrate the applicability of numerical techniques for the determination of AGc for various shaped containers, and different fluids. Numerical integration of the Navier-Stokes equations was used to compute the two-dimensional natural convection in a long, enclosed horizontal channel of rectangular cross section which is heated from below and cooled from above. Fields of temperature and stream function r (where r is defined such that = t, v= v _ - and v yy v ax u and v are the x and y velocity components and v is the kinematic viscosity of the fluid.) were calculated for a number of cases in which the length-height ratio, L/H, was varied from 1/2 to 3/1; the Prandtl number (Pr = C — where c is.K.~C i s

the specific heat, |j is the viscosity, and k is the thermal conductivity) was varied from 0.01 to 25.0; and the Grashof number (Gr 1 g2QH, where g is the acceleration of gravity, 3 is the temV perature coefficient of volumetric expansion and v is the kineri tic viscosity of the fluid) was varied from 60 to 4.0 x 106. The steady state mean Nusselt numbers (Nu = QH/kAQL where Q is the total heat flux), the streamlines and the isotherms were determined for each case studied. For a given L/H and Pr, the Rayleigh number, Ra = Gr x Pr, corresponding to AGc was determined by extrapolating a plot of Nu versus Ra to Nu = 1.0, The values of the critical Rayleigh number, Rac, so determined, were compared with the values predicted by Kurzweg from the simplified equations of motion which have been linearized by neglecting all second-order convection terms. The agreement was within 8 per cent for Pr > 10, Although Kurzweg's linear analysis predicted Rac to be independent of Pr, this was found to be true only when Pr > 10. For Pr < 10 a significant variation of Rac with Pr was observed, and taken to indicate that the assumptions made in linearizing the equations are not valid for Pr < 1.0. Y=H/2 The strength of convection, as measured by I UdY, was found o0 to vary in proportion to (Ra-Rac)l/2, confirming the theoretical predictions of Landau and Chandrasekhar. Rihen convection was studied in rectangles with L/H > 1.0, it was found that several different convective patterns may exist at the same L/H, Pr, and Gr, depending on the initial conditionso This multiplicity of solutions was investigated. It was found that for xiv

arbitrary unsymmetrical initial conditions the preferred mode of convection was that mode in which the number of cells was equal to L/H (for the integral values of L/H used in this work). Other symmetric patterns were forced to exist by using symmetrical initial conditions, but these patterns changed to the preferred mode when subjected to a non-symmetric disturbance of sufficient magnitude. In all cases where unsymmetrical starting conditions, or unsymmetrical disturbances were imposed on the same cavity and working fluid, identical Nu, streamlines and isotherms were found at steady state. To check the assumption of two-dimensional flow, an experimental program was undertaken to obtain a qualitative description of the natural convection in the channel. Unfortunately, it turned out that the choice of stable, or preferred, mode was extremely sensitive to minor perturbations in the boundary conditions such as the illumination, heat leakage, etc., and the equipment did not have adequate control for the desired experiments. Although interesting obser-ations were made, they did not provide the hoped-for critical test of the two-dimensional assumption. xv

I. INTRODUCTION Many situations are known to exist in nature where a small perturbation of a system under consideration will cause complete rearrangement of the system. These situations are termed unstable, and are of great interest to both scientists and engineers, since the presence of an unstable system may cause serious difficulties in the design and operation of process equipment. It is known that an unstable system may remain in an unstable condition for long periods of time, if it is not disturbed. However, a slight disturbance may cause the system to change completely into a more stable form. Thus, laminar flow may be maintained in smooth pipes up to Re = 40,00(46) - much greater than the critical value of 2100. However, if the laminar flow at Re = 40,000 is sufficiently disturbed (such as by a vibration of the pipe), the flow will undergo transition to a turbulent regime, and a considerable increase in the pressure drop will be necessary to maintain the same flow rate. Since the difference in pressure drop between laminar and turbulent flows is significant, the stable form must be known for design purposes. Many of the unstable situations which the engineer is most likely to encounter can be described by the equations of fluid mechanics -- i.e., the equations of motion, energy, and continuity. Since these equations are, in general, non-linear, coupled, partial differential equations, analytic solutions are difficult, it not impossible, except for very special cases. Therefore, other techniques have been developed -1

-2to solve these equations. In particular, the use of finite difference techniques has led to many extremely interesting and useful solutions (2,7,12,44,56) The purpose of this dissertation is to show how finite difference techniques may be utilized in the study and characterization of unstable systems. In particular, natural convection in a long enclosed horizontal channel of rectangular cross section which is heated from below and cooled from above was studied. It is known both experimentally and theoretically that pure conduction with no motion will be the stable condition for a region heated from below and cooled from above, provided the dimensionless temperature difference, or Rayleigh number (Ra), is below a critical value Rac. However, for Ra > Rac pure conduction is unstable and a disturbance may produce convection with a significant increase in heat transfer. The effect of the L/H ratio of the enclosed region, and the Prandtl number (Pr) of the enclosed fluid on the critical Rayleigh number are determined herein. In addition, heat transfer rates, isotherms, and streamline shapes are determined as a function of L/H, Pr, and Ra for Ra > Ra c

II. HISTORICAL BACKGROUND AND LITERATURE SURVEY This section is divided into three parts: 1) In the first part the Benard problem of natural convection between infinite horizontal plates which are heated from below and cooled from above is discussed. Although many of the references are only indirectly related to the problem of convection in an enclosed region, the Benard problem is the classical problem from which other studies have been extended. 2) In the second part, references on convection in an enclosed region are reviewed. 3) In the last part, other numerical solutions of the NavierStokes equations are reviewed. References dealing with numerical techniques and the problem of the preferred mode of convection will be discussed in a later section. 2.1. Natural Convection Between Flat Plates of Infinite Extent The first known reference to the existence of natural convection between horizontal plates of large extent which were heated from below, or cooled from above, appears to be a paper by J. Thomson, (55) who noted the presence of a cellular pattern in soapy water whose mean temperature was somewhat above the ambient temperature. Thomson suggested that these cellular patterns were caused by "convection circulation" -5

-4and thereby opened the whole field of convective instability. (There is modern evidence that the cells which Thomson observed were caused by surface tension gradients rather than thermal gradients.) The next (3) reference to cellular convection appears to be that of Benard, who in 1900 published photographs taken with a beam of parallel light that had passed through a layer of liquid paraffin which was heated from below. In 1916, Lord Rayleigh(42) set forth the basic approach to the analytical treatment for the problem of convective instability in a region heated from below. Rayleigh linearized the equations of motion and energy by assuming that all second order perturbations of the pure conduction solution could be discarded. A solution for the temperature and velocity perturbations was then expressed in terms of a Fourier Series whose components were given by eikx eily. In this manner, Rayleigh accounted for the cellular behavior observed by Benard, and was able to prove the existence of a critical Rayleigh number below which no motion could exist, and above which pure conduction would be unstable. For the unrealistic boundary condition of no slip along the upper and lower surfaces, Lord Rayleigh was able to predict a critical Ra of 654/2. However, he could not determine the critical Ra for the case of fixed, or no slip, boundary conditions. Several refinements of Lord Rayleigh's work appeared.(21'22'26'27) Finally, in 1940 Pellew and Southwell( developed a complete solution for the stability of the linearized equations

-5of motion with either free or fixed boundary conditions on either the upper or lower surface. As with Lord Rayleigh's work, Pellew and Southwell based their studies on an assumed cellular behavior which they characterized by the quantity "a". For a given "a", the critical Ra could be found. It was then a simple matter to find the value of "a" which corresponded to the lowest Ra. For any given cellular configuration, this value of "a" served to fix the cell diameter to height ratio, and therefore, for any given cellular pattern, a preferred cell size could be found. Owing to the linearization of the equations of motion, the hexagonal pattern photographed by Benard could not be explained. In 1935, Schmidt and Milverton(48) performed the first experimental verification of earlier theoretically predicted stability criterion. Agreement with the predictions of Jeffreys(21) was obtained within experimental uncertainty using liquid water as the working fluid. No attempts were made to visualize the flow patterns. In 1938, Schmidt and Saunders(49) used a modification of the equipment of Schmidt and Milverton(48) to study the critical Ra for both water and air. The critical Ra was determined by noting the discontinuity which occurred in the heat flux vs. AG curve as convection set in. Schmidt and Saunders also used a shadowgraph technique to determine the Ra at which laminar cellular motion broke down to turbulent motion.

-6In 1938, Chandra performed flow visualization as well as critical Ra number experiments. He found that for air between plates with a separation greater than ten millimeters, regular polygons seemed to be the preferred mode, while for spaces less than ten millimeters transverse rolls and irregular configurations predominated. Between 1940 and 1954, interest in regions heated from below was limited. In 1947, Stommell(53) presented a general summary of the theoretical work that had been performed to that time, and in 1950 Sutton(5 presented an attempt to explain the results of Chandra(5) for fluid depths less than ten millimeters. In 1946, Jakob(20) reviewed the data of Mull and Reiher, (34) who apparently were the first investigators to obtain precise and useful heat transfer data for air contained between infinite horizontal plates heated from below. In 1954, Malkus(28) extended the experimental work of Mull and Reiher to Ra = 1010, and found six discrete discontinuities in the plot of Nu vs. Ra in the region 1700 < Ra < 1010. In a later (29) paper Malkus attempted to explain the results of these turbulent experiments on the basis of a minimum eddy size which is effective in thermal transport in turbulent flow. He reported that experimental heat transfer rates were within ten percent of those calculated in this fashion, while the predicted transitions were within the experimental uncertainty.

-7In 1957, Ostrach(37) presented a thorough review of the literature dealing with natural convection in regions heated from below. An exhausting review of natural convection in general is provided by Schmidt. (47) The first attempt to predict the heat transfer rates in the region of the critical Ra was made by Malkus and Veronis.(30) The nonlinear equations of motion were expanded in a series of non-homogeneous linear equations which were dependent on the solutions of the linear stability problem. An infinite number of formal, finite-amplitude solutions which satisfy the basic equations were found to exist. The solution which maximized the heat flux was chosen as the most likely solution. Using these assumptions, the authors postulated the heat transfer rate as: Nu = 1 + k(l - a) (2.1-1) Ra where k was estimated to be: k = 1.51 In 1959, Schmidt and Silveston(50) published the results of an extensive experimental study of natural convection between two flat horizontal plates of great extent which were heated from below. Heat transfer rates for various Ra were presented for five different fluids and compared with the air data of Mull and Reiher. The heat transfer results for all six fluids appeared to fall on the same curve of

-8Nu vs. Ra. Shadowgraphs of the motion existing between the two plates were made by shining a bright light through the upper plate (lucite for the experiments) and reflecting it off the lower plate (a polished chrome surface). The reflected light was then photographed. Dark areas were interpreted as regions of downflow and light areas as regions of upflow. In this way the authors were able to distinguish between several different flow patterns that existed at different times and different Ra In 1960, Nakagawa(35) used the exact solution of Reid and Harris(43) for the velocity and temperature field at the point of instability to evaluate the k of Equations (2.1-1). He determined a value of 0.86, which did not compare well with the theoretical results of Malkus and Veronis(30) or the experimental data of Schmidt and Silveston.50) In 1965, K.G.T. Hollands(19) recalculated k from Reid and Harris' solution and obtained k = 1.44, which agreed well with the results of Malkus and Veronis, and of Schmidt and Saunders. (23) In 1961, Kuo3) presented a solution of the non-linear equations of motion for convection between infinite horizontal plates with free surface boundary conditions. The solutions were obtained in terms of orthogonal function series whose coefficients were determined by power series expansions. Once the expansions were obtained Kuo was able to predict temperature profiles, and hence heat transfer coefficients. More recently Herring(17,18) has attempted to predict heat transfer rates, and velocity and temperature fluctuations for turbulent convection between horizontal plates. In these papers Herring investigated only those non-linearities which had the form of an interaction

-9between a fluctuating component and the mean temperature field. Both free and fixed boundary conditions were considered, and fluctuating velocity and temperature fields were characterized by a series of horizontal wave numbers a. 2.2. Natural Convection in Enclosed Finite Regions In 1959, Globe and Dropkin(14) presented data for a circular region heated from below and cooled from above, but with such a high Ra that the motion must have been turbulent. The first investigation of two-dimensional convection in a region with H/L near unity appears to be the numerical work by Deardorff(7) in 1964. Deardorff used finite difference techniques to integrate the unsteady Navier-Stokes equations for an enclosed rectangular region heated below with fixed boundaries on all four walls. Length-height ratios of 1:1 and 2:1 were studied. The final solution was found to be independent of the initial conditions for the six cases studied. Steady state streamlines and isotherms were presented for the six cases investigated. Since Deardorff was attempting to describe turbulent convection, the lowest Ra he investigated was 6.75 x 105. The assumption of two-dimensional flow was found to suppress the appearance of random turbulent fluctuations, and in most cases a steady state solution was obtaind. In 1965 Kurzweg(24) developed a stability criterion for natural convection in an enclosed region, heated from below, with rigid sidewalls. The two-dimensional Boussinesq approximations to the

-10 - Navier-Stokes equations were linearized and a solution was assumed for the temperatures and stream functions in the form of two infinite series of orthogonal functions. A modified Fourier technique, as out(41) lined by Poots,(41) was used to evaluate the coefficients of the expansions. In order that these coefficients be non-trivial it was found that an infinite secular determinant (whose only unknown parameter was the Ra ) must equal zero. By terminating the series for T and at a finite number of terms, (N=3) the determinant reduces to an 18x18 matrix whose solutions give the various values of Ra for which a nontrivial solution exists. The lowest of these values is then taken to be Rac. Kurzweg presents values of the critical Ra as a function of the L/H for L/H between 0 and 4:1. Plots of the reduced stream function */*max and reduced temperature perturbation 0/Omax were presented for L/H = 1., 2., and 3. Absolute values were not available because of the linearization of the equations. (13) More recently, Fromm has presented the results of several numerical solutions of the two-dimensional, Boussinesq approximations to the Navier-Stokes equations for natural convection in an enclosed region heated from below, and subject to either free or rigid boundary conditions on the horizontal boundaries. Cyclic boundary conditions were used on the vertical boundaries in an effort to simulate the cellular behavior found experimentally between infinite flat plates. The results of Fromm's calculations for free horizontal surface boundary conditions agreed well with the theoretical predictions of Kuo.(23) Fromm also presents results for rigid boundaries in a cell seven

-11units long by three units high. Various cellular patterns were attempted, and results were given for two-, three-, and four-celled patterns. It was noted that a five-celled arrangement was unstable and broke down to a three-celled pattern. Approximately equal heat fluxes were found for the two-, three-, and four-celled patterns at steady state. Aziz(l) has recently presented the results of a numerical study of natural convection in a two-dimensional square region, and in a cube. Results are presented for Pr = 1.0 and Pr = 7.0 at various Ra > Rac. The three dimensional equations of motion were simplified by use of a three component vector potential stream function and a three component vorticity vector. The equation of motion and continuity then reduce to three equations —one for each component of the vorticity. The three equations were solved by standard numerical techniques. In 1961, Sorokin(52) presented the results of an experimental study of natural convection in a horizontal cylinder heated from below. Sorokin's experimental apparatus consisted of a long lucite block of square cross-section in which a circular hole had been drilled. The hole was filled with water, and some aluminum dust and sealed. Thermocouples were located in and around the hole. The whole block was heated from below and cooled from above, and heat transfer rates (as measured by dO/dR in the lucite block) were determined as a function of AG across the hole. The critical A0 was determined by plotting the heat transfer rate versus AG. Photographs of the dust patters show that the fluid may have either of two basic modes: 1) a two-dimensional circular rotation or; 2) a three-dimensional cellular motion.

-12 - 2.3. Other Numerical Solutions One of the first successful numerical solutions of the NavierStokes equations was by Hellums and Churchill5) who investigated free convection near an isothermal, vertical plate, and natural convection within a horizontal cylinder whose two vertical halves were at different temperatures. The computed results for the hollow cylinder agreed (31) satisfactorily with the experimental results of Churchill and Martini. (12) Fromm() investigated the unsteady wake behind a small rectangular body placed parallel to the flow between two moving flat plates. He was able to predict the fully developed vortex street behind the blunt obstacle, and by the use of an extremely fine mesh size was able to describe the behavior of the wake region in great detail. Richards(4) studied the effect of a point sink located at the center of an infinite rotating horizontal disc. Streamline patterns were calculated as a function both of sink strength and rotational speed of the disc. Wilkes )studied two-dimensional natural convection in an enclosed rectangle with side walls at different temperatures, and either a linear temperature gradient or no heat flux at the horizontal surfaces. Results were presented for a series of different L/H ratios and Prandtl numbers. Barakat and Clark(2) have studied transient natural convection of liquids in closed containers with a gaseous layer above the liquid. In particular the effect of various surface heat flux and temperature boundary conditions on the transient natural convection in the liquid were determined.

III. DEVELOPMENT AND DESCRIPTION OF THE MATHEMATICAL MODEL During his experiments with natural convection in a horizontal cylinder heated from below, Sorokin(52) discovered the presence of two distinct modes of natural convection. As the temperature difference was increased beyond the critical A Q, the first pattern, in which the fluid was observed to undergo a planar rotation in the Z. plane, developed. No component of velocity was noted in the third dimension. As the temperature difference was further increased, this planar motion broke down, and a three-dimensional cellular motion evolved. Similarly, in the shadow graph pictures presented by Schmidt and Silveston(50) of natural convection between infinite horizontal plates, a two-dimensional roll pattern appeared to be the preferred mode of convection for Ra slightly greater than Ra. The question of the preferred mode of natural convection has recently come under intense study. In 1960, Palm(39) presented an analysis of the non-linear equations of motion in an attempt to account for the effect of the variation of physical properties with temperature. He concluded that hexagonal cells should be the preferred mode of convection if the variation of kinematic viscosity was of a sufficient magnitude. In 1963, Segel and Stuart(51) extended and corrected the work of Palm. By assuming a particular form of viscosity variation with temperature, they demonstrated that hexagonal cells should be the preferred mode provided the kinematic viscosity variation is great enough. They also showed that two-dimensional rolls could be the preferred mode of convection if there was a sufficiently small variation of kinematic viscosity, or if the conditions during the onset of convection were just right. -13

-14It is apparent that the a priori prediction of the mode of convection for a given set of conditions is difficult, and not yet resolved. Fortunately, it appears that the various modes of convection all result in approximately the same heat fluxes. Thus, for the predictions of heat transfer rates almost any stable mode of convection may be assumed. (This point is well demonstrated in the results of Fromm(13), in which two-, three-, and four-celled two-dimensional patterns were found to yield approximately equal heat fluxes.) In the calculations presented herein it is assumed that the motion is two-dimensional. This mode is known to be physically stable, and therefore the calculated results should have physical significance in so far as the other assumptions are valid. Whether or not the twodimensional motion is more preferred than three-dimensional motion is unknown. Fortunately, the lumped parameters such as the total heat flux are essentially the same for the two- and three-dimensional modes near the critical Ra. Hence, the values of these lumped parameters for the two-dimensional case are good approximations for the three-dimensional case as well. It is indeed the similarity between the properties of the various modes which makes their relative stability so difficult to determine. Three-dimensional calculations are theoretically feasible by an extension of the methods presented herein, but will require several orders of magnitude more computer time. It should also be noted that in three dimensions the complication arises that the most preferred dimensions of the circulation are not known unless the length, height, and width of the region are approximately equal.

-15Under the assumption of two-dimensional flow, the physical problem may be stated as follows: A fluid of mean temperature o0 is at rest within the enclosed rectangular region O < x < L, and 0 < y _ H. (See Figure 1). The fluid initially has a linear temperature gradient with 0 = gh at the lower plate y = H, and Q = Q1 at the upper plate y = O. It is assumed that Qh > Q1. A small circulation is introduced into the central region and allowed to develop subject to the boundary conditions: Q = Qh on the lower surface; 0 = 1 on the upper surface; Insulated walls on both sides; and u = v = O on all walls. x 91 x=L y Oh + 1 0 28x 8x y=H With the initial conditions u(x,yO), v(x,y,O) = f(x,y) where f(x,y) represents the initial velocity fields. Figure 1. The Physical Problem

-16The temperatures, velocities, and heat transfer rates were calculated as a function of the time until steady state was approached. Since only the steady state solutions are of interest, the starting conditions may be changed to speed the approach to steady state. For the purpose of simplicity, the fluid properties will be assumed to be independent of temperature except for the buoyancy term. In particular, the fluid viscosity, thermal conductivity, and specific heat are assumed to have constant values 0, ko, and Cp these being the actual property at the mean temperature 0 = 2 (Gh + 91) Density variations due to pressure are neglected, and the density at any temperature ~ is related to its value at Q by the equation of state: p = po/(l + P[G-G0]) (3-1) where P is a constant. To a first approximation, P = V/V and has the value of l/Go for gases at low pressure. 3.1. Equations Governing The Problem For the geometry under consideration, the equations of motion, energy, and continuity may be written as: -+u + -+V -+- = - ^au. a v au -1 a a2U +(u; (3.1-1) t a x +v y P ax V x2 +y2 av av av 1 aP 2v d2v +- u - v = =g- - +v -2 + (3 1-2) at ax. y p dy ax dy

-17— + u -+ v K (31- =) - -.v^ =K ^+ 1i (3.1-3) at ax by ax2 by2 6+U6 a + v + p + 0 6t 6x 6y x dy / However, the equation of continuity, Equation (3.1-4), is simplified since constant density has been assumed except in the gravitational term of the equation of motion. Therefore: u +v = (3.1-5) ax 6y Now assume that: P = Po + pl (3.1-6) where p = actual pressure at any.point, and Po = the pressure which would exist at any point if there were no motion and all of the fluid was at temperature g0. Therefore, equation (3.1-2) yields: o = g -1 9 (3.1-7) Po ~y Equation (3.1-6) yields: 6p =8;6l (3.1-8) ax ax and therefore, g - g [1 + ( - ] ( pg + 1 (3.1-9) g P oy p y

-18g -(9 - ) - 0)) _ __ g ap =g( g ) ( + ) (3.1-1 0) If the assumption is made that: P(g-~o) < < 1 Equation (3.1-10) becomes: g 1 = gp(-) - (3.1-11) P ay Po ay Substitution of Equations (3.1-11) and (3.1-8) into (3.1-2) and (3.1-1) respectively yields the equations: u u au 1 up1 + v-2 + 2u -+u - + - ---- = -y (23.1-12) at ax ay PO 6x x2 ay2 v +ua+ v = g(Q-o) -1 1+ v ( -2v + (3.1-13) at ax a p y Pxy x y2 which, along with Equations (3.1-1) and (3.1-4) and subject to the initial and boundary conditions: u(x,y,0), v(x,y,0) = f(x,y); 6(x,y,0) = 1 + H (Q-) O(x,0,t) = 91 G(x,H,t) = 0h (3.1-14) u = v = 0 on all surfaces ao(o,y,t) o ax 6a(L,yt) ox

-19constitute a complete statement of the problem. 3.2. Dimensionless Form of the Equations Following the precedure developed by Hellums and Churchill(l6) dimensionless variables of the following form are postulated: T = t; T = 2(-Qo); X =; = t h-1 x x (3.2-3) 1 v P = 2; U = u; V _ v pt u v o o o where Q = (Gh+ @1) /2. These are substituted into Equations (3.1-3), (3.1-4), (3.1-12) and (3.1-13) to yield: ou auuo - - V + + V2U (3.2-2) to aT X0 ax aY poo ax x2 uo av u2 v av I gp(h-l91) P ap - T + 0 5- + V po=0- T 0 to aT xo x 6Y 2 poXo ax v u o + 00 V2V (3.2-3) (Gh-Gl) 6T u (Gh-Gl) T K( 2( - K(h-l)) V2T (3.2-4) 2 t OT 2 x I X OY 2 x2 o0 0 O -+t -v ) 0 (3.2-5) whrXo v -02 a 2 2 = where V2 =, + oX2 oY2

-20with dimensionless initial and boundary conditions: U(X,Y,0), V(XY,0) = F(X,Y) T(X,Y,0) = - 1 + - 2.Y H T(X,0,T) = - 1 T(X, H ) +1 (3.2-6) U = V = 0 on all walls aT(O,Y,T) = ax aT( L,Y,T) = _X X ax "o Equations (3.2-2) to (3.2-5) may be rearranged to: Xo 6U b U U pi P v x- — +U —+ V-= - 2 - + V2U (3.2-7) too bT ax bY pOu2 X xou x av av av xTgB( ) p- _, O+ U + V T- touo aT x aY 2 uO pOu oY + v V2V (3.2-8) xouO -T + U T + to V2T (3.2-9) aT xo aX bY xo 0 0 Therefore: U,V,T = o O xogp(Qgh-1) touo P o2 xou 2 u2 K0, (x L2 IV H, _o (3.2-10) H x o] 0 o

-21Which is simplified by setting certain groups equal to unity. let 1 x = H H o Vo = 1 V let = u~ = xuo H (3.2-11) let - = 1 to = H / touo let o = 1 p = V2 let 0 nPO 2 Puo H Substitution of these relations in Equation (3.2-10) yields the relation: U,V,T = ( [Pr,Gr,L/H] (3.2-12) where K = 1; H3g(Gh-91) = Gr v Pr 2 v2 Thus the dimensionless form of the equations of motion are written as: _U + u +v - + v2U (3.2-13) aT ax aY ax +U +V v Gr T - + (3.2-14) 8T 1 X Y /Y + U + V -= 1 V2T (3.2-15) aT ax aY Pr aU/aX + V/Ay = 0 (3.2-16)

- 22These together with the dimensionless initial and boundary conditions: U(X,Y,O), V(X,Y,O) = F(X,Y) T(X,Y,O) = —1 + 2.0Y T(X,O,T) = - 1 U =V = 0 (3.2-17) -T(O,,T) = 0 on all surfaces ax aT(LYT) = ax constitute a formal statement of the problem. 3.3. Evaluation of Nu The local heat flux is evaluated from: q = - k a |(3.3-1) while the local heat transfer coefficient is defined by: h = -q (3.3-2) (Gh-il) and the local Nusselt Number, Nu, is defined as: (3.3-3) Nu = - where xo = a convenient reference length. k Therefore: Nu = XO Io| (3.3-4) 4h-@l 6y Yo 0 in terms of the dimensionless variables Equation (3.3-4) becomes: Nu = 1 T (3.3-5) Th-Tl aY Y= but Th = + 1, T1 = - 1.

-23Therefore, the local Nu is: Nu 1 - (3.3-6) 2 IYY= 0 The overall Nu is defined as: Nu = - Nu dx L 0 which becomes: L Nu = - 1 H f at dX 2L a6Y Y=O 3.4. Elimination of Pressure From The Equations of Motion and Introduction of Stream Functions and Vorticites The pressure terms in the equations of motion may be eliminated as follows: 1) take a Equation (3.2-13) aY 2) take - Equation (3.2-14) ax 3) subtract 2) from 1) Rearrangement of the results of 3), followed by elimination of certain terms through use of the continuity equation, gives the following relation. F V 6U + u a [V - L+ V [ aV -u 6T l ax YJ ax l x YJ XY [lx Y J T ax 6 Lx - Y For simplicity, the term [ - U is defined as the vorticity and given the symbol 6 = _V - U.The equation of motion is then ax aY

-24written as: t + U X +V X= -Gr T + V2 (3.4-2) OT oX 3Y X (3.4-2) A new variable, the stream function r, is introduced so the continuity equation is always satisfied. r is defined such that: - = U aY and (3.4-3) a._ = v ax which always satisfies the continuity equation. The vorticity is expressed in terms of stream function as follows: aV a U - v2* (3.4-4) and the equations of motion and energy are now expressed as:.- + u a + v X = - GrL + v2~ (3.4-2) aT aX aY dX a+ U + a+v = 7 2T (3.4-5) OT 6X oY Pr The initial and boundary conditions for Equations (3.4-2) to (3.4-7) are given by: /(X,Y,O), S(X,Y,O) = F(X,Y) T(X,Y,O) = - 1. + 2. Y T(X,o,T) = - 1 (3.4-6) T(X,1,T) = + 1 aT(O,Y,)) = aT(,Y,T) = 0 ax ax H

-25The stream function boundary conditions on the four surfaces may be derived from the condition U = V = 0 on any surface. For example, examine the surface X = O. U = V = 0 implies X = - = 0. X = 0 can be inteaY aX aY grated to 4(O,Y) = constant. This constant may have any value as long as the same value is used on all four surfaces. For simplicity (O,Y) = 0 is chosen. Therefore, the stream function boundary conditions become: =|X= 0 (3.4-7) an on any surface, where n represents the direction normal to the surface. The boundary conditions (3.4-6) and (3.4-7), combined with Equations (3.4-2) to (3.4-5) constitutes the formal statement of the problem in terms of stream functions and vorticities. In Equations (3.4-2) and (3.4-5) the terms U and V have been included. These terms could just as well have been replaced by their equivalent expressions in derivatives of'. However, the computational effort will be slightly reduced if U and V are evaluated once for each set of 4, rather than evaluating the respective derivatives. Physical significance can be attributed to both vorticity and stream function as follows:(-2) The vorticity ~ represents twice the rate of rotation of a fluid element, while the stream function is usually represented by plots of constant $, or streamlines. The streamlines are tangent to the velocity vectors, and at steady state represent the path of a particle "floating along with the fluid."

-263.5. General Logic of the Numerical Solution A basic approach to the numerical solution of equations may be visualized as follows: 1) Initial conditions in the form of 4 and T are read into the computer. 2) Suitable finite difference approximations (FDA) to Equation (3.4-4) are used to calculate the vorticity field. 3) Suitable FDA to Equation (3.4-3) are used to calculate the U and V fields. 4) Using the U, and V fields previously calculated, a finite difference approximation to Equation (3.4-5) is used to advance the values of T. 5) With the new T, a FDA is used to advance the values of the vorticity. 6) With the new vorticities an iterative technique is used to solve the simultaneous FDA to Equation (3.4-4) and thereby advance the value of. 7) The program returns to step 3), and steps 3-7 are repeated until the desired number of time steps have been attained.

IV. FINITE DIFFERENCE APPROXIMATIONS AND THEIR STABILITY 4.1 Formulation of Finite Difference Approximations In the numerical solution of a differential equation the dependent variables are assumed to exist at a finite number of regularly spaced values of the independent variables known as grid points. The derivatives are then replaced by finite differences and the differential equations are replaced by a series of algebraic equations which are then solved by conventional means. Before attempting to derive FDA to the various derivatives that will be used, it is necessary to define the values of the independent variables at which the dependent variables will be specified. Therefore, the following grid is defined: where: i=0 1 2 3 NX-2 INX-1 i=NX O < Y < i.1 1 1 I 1 1, l < X < L - H 2 m- - — L —-- - - H NX 1.0 N -1 - - --. i ine 3v X = igX Y = jAY \X = space Ay NY-2 between vertical NY-1 grid lines AY = space +-~-j AX |*- between horizonFigure 2. The Numerical Grid. tal grid lines -27

-28 - The intersection of two grid lines serves to fix the locations of the grid points. The following notation will be used to represent the value of a dependent variable: inj (iX, jAY, n lT) (X, Y T) Once a grid system is established Taylor's Series Expansions about any grid point may be used to approximate derivatives with finite differences as follows: Suppose i j is known for all i,j and it is wished to approximate X 2 2 etc. The Taylor's series in $ may be written o X' X2 as:'i+lj = ij + a X AXX +... + o(~x3) (4.1-1) akX P 2.i 6X2 ^ -— ^j = i 1Mj - __L + _2i' - O(Z6X3) (4.1-2) whre oAx)1! 6X 2! 6X2 where 0(zX3) represents the remainder term of the series, and is read as "Order of 3". If SZX is small so that AX2 < < 1 we may terminate Equations (4.1-1) and (4.1-2) at the second terms to yield: oi+l,j = ~i,j + ( j + (X2) (4.1-3) i-l,:= -ij - AX ~^ + o(AX2) (4.1-4) which can be rearranged to yield:. _ +l. + o(ax) (4.1-5) ax ax

-29ij = 4i,-lj + o(Ax) (4.1-6) ax /x X'LAX Equation (4.1-5) is termed a forward difference approximation, while Equation (4.1-6) is termed a backward difference approximation. If greater accuracy then 0(AX) is desired, Equation (4.1-2) may be subtracted from Equation (4.1-1) to yield:,.. i+lj'i-lj + o(ALx2) (4.1-7) ax 2AX which is termed a central difference approximation. In a similar manner a FDA to - may be derived by addition aX2 of Equations (4.1-1) and (4.1-2) to yield: - _ X l-2J (4.1-8) 2a AX2 A forward difference approximation to - is given by: aT (=, j - Xj)/ (4.1-9) In a similar manner approximations of higher order accuracy could be derived by including more grid points, and more equations. However, in practice it is found desirable to avoid the use of higher order approximations for two reasons: 1) FDA which use only three adjacent grid points generate a special set of equations whose solution can be obtained by non-iterative techniques. 2) High order FDA are often plagued by numerical instability, a condition in which small errors become magnified to unreasonable proportions. Therefore, only those FDA which have already been resented will be used for the formulation of FDA to differential equations. However, when evaluating

-50coefficients for these equations, such as U and V, it is advantageous to use higher order approximations. Once suitable FDA to the various derivatives have been formulated is is possible to convert the partial differential equations (PDE) of Vorticity and Energy into a series of algebraic equations. Several different methods have been proposed for this conversion, and may be broadly classed as: 1) Explicit 2) Implicit, and 3) Alternating Direction Implicit (ADI), depending on whether the non-time derivatives are evaluated at T = nAT or T = (n+l)AT. Examples of each class can be demonstrated by use of the heat conduction equation: - = -. 2+ y2 (4.1-10) 6T ~X2 +Y2 Explicit: hn =2 Xn 2+ n 5T~ ij Xuiij -$ i~j where 52= i-l,j 2- i,j + Pi+l,j and n -n, j L\X 2 ( j \AT In the explicit formulation n+ can be calculated directly for each grid point once the n are known. Implicit:,Tj1 = 2 n+l+5 (4.1-11) ~~. ~~. +s ~ ~.1

-31In the implicit formulation it is not possible to calculate n+l n+ fore, a system of simultaneous equations must be solved for i.j. The solution of this set of simultaneous equations is usually performed by an iterative technique such as the Gauss-Seidel, or successive OverRelaxation iteration. These methods, however, usually require a significant number of iterations for satisfactory convergence, and therefore require a great amount of machine time. Alternating Direction Implicit: In order to overcome the necessity of an iterative procedure Douglas (8), Peaceman and Rachford (39), and Douglas and Peaceman (9) proposed the ADI procedure. This method requires the line-by-line solution of small sets of simultaneous equations which can be solved by non-iterative techniques, rather than the iterative solution required by the simple implicit formulations. In the ADI procedure the time step is broken into two equal segments. In the first segment all the derivatives in, say, the X direction are taken as implicit, while those in the Y direction are taken as explicit. In the second half step the process is reversed, and all X derivatives are explicit while all Y derivatives are implicit. In this manner the systems of simultaneous equations formed are of a special type known as "tri-diagonal", and have a solution which can be expressed in closed form (See Appendix 1).

-32The ADI procedure is illustrated by the following two steps: n A )/nT/2 =x + qf An (l - ij/ /2 2 j +2 n+l (4.1-12) (~)n * )/A//2 = 52 * 2 3n~l where $i,j is an intermediate whose value is meaningless without the second step. Unfortunately, this ADI scheme cannot be directly extended to three dimensions without encountering stability limitations. To overcome this weakness, Douglas and Rachford (10) developed a threedimensional version of the ADI procedure. In an effort to increase the accuracy (lower the truncation error) of Douglas and Rachford's three dimensional ADI procedure, P. L. T. Brian (4) developed a three-dimensional ADI version of the Crank-Nicholson procedure. It can be shown that Brian's three-dimensional procedure will reduce to Equations (4.1-12) when used in only two dimensions. Douglas and Rachford's three-dimensional scheme may also be used for two dimensions but involves a higher truncation error than (4.1-12). For completeness the three-dimensional procedures of Douglas and Rachford and P. L. T. Brian are given below: Rachford and Peaceman: 1 * n 2 2 n -- j_~n~j, b 52 *+ y + + T- [ -i]= -X Ij + 2Y l + bZ tlj (4.1-13) 1 n+ln 2 * 2 [itj~~i,j] = xi,j + +

-33P. L. T. Brian: 2 * n 2 2 2 n A-T [i,j- i,] = 5X 4i,j + by *ij + 5Z ij a) A i i, ij] = SX ijj + + Z i,j b) (4.1-14) i) j~, ^,^ X = Pi j + b il + 5 hj c) 1 [n+l n 2 2 2 A- i,- j] *= 5X i j + 5Y i,j + Z Pij d) iyj can be eliminated from Equations c) and d) so that only evaluation of *, ** and $n+l are necessary. 4.2. Stability Considerations At first glance it appears that explicit formulations would be preferred since they require the smallest number of computations per time step. However, it can be demonstrated that, in general, explicit FDA suffer from numerical instability if a critical time step is exceeded. This time step is often quite small and therefore a large number of time steps is required. Implicit methods, including ADI, on the other hand, require a greater number of calculations per time step, but usually do not have stability imposed limitations on the size of the time step. Therefore, a lesser number of time steps is required to reach the same value of time. A balance between the number of time steps, and the calculations per time step must be struck.

-54However, before this decision can be made, it is necessary to develop a rigorous definition of stability, and a criterion for determining the maximum allowable time step. 4.3. Definition of Stability Richtmeyer (52) introduces the concepts of convergence and stability as follows: a) convergence. Let $(X,Y,T) be the exact solution to a differential equation, and ~n j be the approximate numerical solution to the differential equation. Convergence of a numerical approximation is taken to mean,n j - (X,Y,T) as n- oo. A theorem by Lax (53) is introduced which guarantees the convergence of a FDA to the true solution provided the approximation is both consistent and stable. The consistency condition requires that the truncation error of a FDA must approach zero as AX, AY, AIT - 0. (A condition which is automatically satisfied by all FDA used in this work.) If the consistency condition is satisfied then stability is the necessary and sufficient requirement for convergence of the approximate solution to the true solution as AX, AY, AT - 0..b) stability. If the difference between ni j and ((X,Y,T) is called the error of the approximation two pertinent questions about the value of this error may be asked: 1) what is its behavior as n - oo for fixed AX, AY, Ar 2) what is its behavior at fixed T if AX,AY,AT -e 0.

-35In many cases that the answer to Question 2 depends on the relative rates with which AX, AY, AT -> 0. This question will be regarded as the more significant of the two since it is hoped, in the limit, to reduce this error to zero by letting AX, AY, AT -> 0. For both questions it is seen that an increasingly large number of time steps are required, and therefore the possibility exists for unlimited amplification of errors. If this occurs the final numerical approximation to 4(X,Y,T) will be completely meaningless, and therefore, Richtmeyer states the desired stability criterion as follows: A FDA is stable if, in the limit as n - oo, the error i (X,Y,T) remains bounded. Since, in general, $(X,Y,T) is bounded, the stability criterion can be rewritten as: A FDA is stable if in the limit as n - oo the approximation n pi, j remains bounded. 4.4. Determination of the Stability of a Finite Difference Approximation Two methods for determining the stability of a FDA will be established. It should be noted in particular that these stability criteria are sufficient, rather than necessary. That is, the FDA must be stable if the conditions are met*, however, if the conditions are not met one cannot definitely state that the FDA will be unstable. Specific examples can be sited of a stable FDA which does not satisfy * This analysis does not take into account the effect of boundary conditions which often introduce instabilities of their own.

-36one or the other of the stability criteria. In cases where the two procedures yield different results the less severe criterion is used since it guarantees stability. Experimental calculations have shown, however, that when neither of the criteria is satisfied, the FDA will generally be unstableo The stability criteria will first be applied to the simple heat conduction equation for explicit, implicit, and ADI formulations, and the results of the two tests will be compared. Then the concept of explicit, implicit, and ADI formulations will be extended for the vorticity and energy equations. The stability of the three formulations will be examined and the advantages and disadvantages of each formulation will be discussed. Next the FDA to the stream function equation, Equation (3.4-6) will be given, and finally the various FDA to U, V, and Nu will be derived. 4.5. The Method of Positive-Type Equations The simplest form of stability test makes use of the positivetype FDA as follows: An explicit finite difference equation of the form n-t-l ai~j~ ijn n n l j= + + a ij + aii+, + aii- l, + aij+li,l + ai-l^iJ (4.4-1) is termed positive-type if all the coefficients ai,j, ai+l j,.are nonnegative, that is, if ai, > 0 for all i,j ~? j

-37Similarly, an implicit finite difference equation of the form: n+1 n n+l j n+ 1 n+l + 1 inj = jbinjij + bi+l,j4i+l,j + bi-ljti-l,j + bi,j-lij-l + bij+i (4.4-2) is termed positive-type if bi,j > 0 for all i,j (any simple conduction or convection equation can be put in these forms where ai j or bij may, or may not, be > 0.) Barakat and Clark (2) have shown that explicit, positivetype FDA are stable. Similarly, Forsythe and Wasow (11) have demonstrated the stability of positive-type implicit FDA. Let us examine the stability of the explicit form of the heat conduction equation: x2" 5 Y2 - aT which can be put in the form of Equation (4.4-1) as: n+ 2zl n AT( n n n Qinj=[i- 2( 1+ m|)](P"' t i+1 j+< )+A + ) ZXX2 AY 2 CSX2 Y A)y2 i',j+l i -1) (4.4-3) For a positive-type equation: 2A< / AX2, X ( + A2) <1. (4.4-4) In a similar manner the implicit formulation can be put in the form:

-38n+l. n+l F 1 n+l 2X (Al+ 1 11+1 in, b j3 bAX2 ni+lj + 1 j + 1 i + 1ij -l) i,1j b 12 (4L4-5) where 2AT AX2 b = 1 + - (1 + AX2 AY2 Since all the coefficients are > 0 irrespective of AT, the implicit formulation is unconditionally stable. 4o(. The von Neumann Stability Criterion Unfortunately, the concept of positive-type FDA does not yield the most useful results when applied to multi-level approximations such as those of Crank-Nicholson, Dufort-Frankel, etc. For the study of these formulations a more general procedure than provided by the positive-type equations is needed. Therefore attention is turned to the von Neumann Stability test, as first reported by O'Brien, Hyman, and Kaplan (36). The method of von Neumann is strictly applicable only to linear equations with constant coefficients. Following von Neumann the difficulty of non-linearities, or non-constant coefficients can be circumvented by successively applying the procedure to a series of small overlapping regionso Each region is assumed so small that the coefficients can be assumed to be constant within ito Such a procedure, while it has worked well in practice, is, of course, open to serious question.

-39The method of von Neumann is applied by assuming that a solution of the form $(X,Y,T) = Z e7 ei(o+Y) exists for the finite difference equations. (y is a function of a and a and may be complex). One term of the solution is then expressed in the n n form p,(X,Y,T) = (, )eei(X+Y) where = eTT is the amplification factor for the a,P component, and may also be given by n+l ap, = Ba,/. The term is then substituted into the finite difference equations and a, is determined as a function of AX, AY, AT, a and P. Since it is required that $(X,Y,T) must be bounded as n - m, then ( (X,Y,T) must be bounded for all a,p as n -+ X, and the von Nuemann stability criterion is expressed as I, | _< 1 for all a,. If simultaneous PDE are being considered, the method is almost analogous. If the variables under consideration are Z and T n i (CYX+SY) for instance, general terms of the form Z a, = ia, e and T, a =,n ei(OaX+PY) are assumed and then substituted into the original simultaneous FDE and a relation of the form: = Ic j l =[Aj (4.5-1) is derived. The matrix [A] is called the amplification matrix. The von Nuemann stability criterion for this case may now be simply expressed as A < 1. This condition may be reduced to X1, X2 < 1 where X1 and X2 are the eigenvalues of the amplification matrix. For a proof see Wilkes (56).

-40The von Nuemann stability procedure will now be used to study the stability of the heat conduction equation: 4 S 6 + 6 ax2 + 2 Explicit Form: The explicit FDA to the heat conduction equation may be writtem as: $n+l 2n n n n ng Cj m nj.... (i+l^j-2$nj+1n 1 j+ r(in j+l-2i ~,j+inj-l)) (45-2) where b r = Assume that the general term in I, 2 - ZX= AY2A can be expressed as: n n, ~ e icoX t e (405=3) Substitution of Equation (. 5-3) into (4o5-2) and division by ei Xei1 Y yields n+l n ^aB - Sag' s ~E [culX_~i~`ic~R~i~~(~iasi-^^^^-Y~15^~B, (4.5-4) Therefore: i = [ 1n+l, = 1 + 25 [( —r) + (cos ax + r cos BaY)] (405-5) and the condition for stability becomes: -1 < 1 + 25 [l-r) + (cos oAX + r cos PAY)] < 1 The right hand inequality is always satisfied; the left hand inequality

- 41will be satisfied for all a,p if, and only if: ATX AY 2 1+ —Lh2 A\ y2 - < 1 Implicit Form: Substitution of the general term (n = n ioX i3Y,= - a3 ei e into n+il n,+n+ll n+l n+l - n+l fij - i+j = 5 i+lj 2i, + i-lj + r(<i,jl j + in +l) (4.5-7) yields: n+l1 n 1 ^PI/5a) J=', = l-26[cos(cAX) - 1 + r(cos(AY) - 1)] Since the largest positive value of (cos(oAg) - 1) + r(cos(AY) - 1) is 0, the greatest value of \ K a| < 1, for all A, and hence the implicit formulation is unconditionally stable. It should be noted that the stability criteria predicted by the von Neumann and positive-type methods are identical for these onelevel approximations, ADI Formulation ~i n [ * * * +n- n' )iJ ij = 25 [i+l,J 2- i + 4i-l,j + r($i,j+l-2iJ+i,:,J-1l) (4.5-8). -. - - 2i. + + r(. +1- 2.+.+) SC S SC_ in+,n~ljn+

-42Let i> be amplification of step 1, /; /, is the overall Lt,p Op amplification n++l/n.n The stability criterion is, < 1 Assume: *,j = ei((X+PY) a) $n+l =: n+l ei(OX+PY) b) Substitution of -) into step 1 and division by the common factor ei(X+PY) yields: */n *= - 1 + 45r(cos(fY) - 1) (45-10),fp 1 - 4-br (cos(AlX) - 1) Substitution of b) into step 2 and division by the common factor yields: n+ l/ = 1+ 48r(cos(c-xX)-l) (, 4 1 - 45r(cos(aY) - 1) The overall amplification factor is then: p = n+l/n = 1l + 4br(cos(_a X) -1)1 1 + 4br(cos(PAY) -1 Iaf h\1 - 4br(cos(fY) -l) 1 1 - 4&r(cos(aaX) -1)) which is always less than. for all AX, AY, AT, ca, and 5o Therefore the ADI formulation is unconditionally stable. 4.7. Finite Difference Approximations to Equation of Motion and Energy and Their Stability The following FDA may be applied to the vorticity and energy equations: 1) Explicit:

-43Vorticity: in+l i.n n+l 2 n 2 n n n AT - GrbXTij + 8X Sij + bY i, j -Ubx i j - Vby ti j Temperature: Tn+l Tn i, ij 1 2 n 1 2 n AT = p bX Tij + r 8y Ti j - UX Ti j -VSy Tij where 5X could represent either a forward, backward, or central difference derivative. As will be seen, 5X and by will have to be either a forward or backward difference, depending on the sign of U and V, if a stable formulation is to be obtained. 2) Implicit: Vorticity: n+l n iAT X l X X +ib J + - X - U n+l YV SIn (4.6-2) Energy: /XT - G 5X Tn+ + T Ux TX ~ TJ rpn+l n -Ll^J~ll^m 1 52 n+l1 2 n1 n+ l n+l AT + p - VbX T Where bX and by have the same significance as in the explicit case. ADI: (P.L.T. Brian Version) The vorticity equation will be discussed first:

-44lj = - Gr b Tn+l + 2 j + - U V a) 2p' ^ l nn+l 2 2 A7 ^,j-^j j = - Gr 5 T + 6x %i,j + by i,j - Ux ^ij - Vby ^, b) 1 [n+l n Gr 5X Tn+l + 52 * 2 3 V5 Equation a) is the first half step for all ADI methods. The * second half step is obtained by eliminating e from Equations a r3n n+l ]; ^ = _1 bY(i j^:+) - ~v^(gj ) (4~6-5) 2 for the second step. This formulation has the advantage of requiring slightly less calculations per time step than the normal ADIs The ADI equations are then: Step 1i ILJ 1,J2 * 2 * (46 6) Step 2: n+l n 2 iby jnijn+1i Vb, n _ n+l1 A T = 2 (,JJ - 2 1,J 1,J In an identical manner the temperature equations can be derived as:

-45Step 1: * n 2 2 -T, ji j 5X by. iT = Pr (T ) + (T ) - U Ti j - Vby T (4.6-7) Step 2: 2T? -T T n+Vl ij 1, j 1, j 1 2 T, Tn+ V) AT 2Pr bY 1(ij-ij - 2 (Tij1-Ti where 6X and by represent central difference derivatives. Since both the explicit and implicit formulations involve single level approximations to v2 and 5X, 6y, the method of positive type equations will be used to predict their stability criterion (the same results can be derived using the von Neumann method). Explicit: For purposes of example we will use a forward difference approximation to 6X ~, a backward difference for 5y 5, and central difference for 5X T. With these approximations the explicit form becomes: ijlj- ij Gr _ n l mnl n n n AT = X Ti+l,j [i+ + AY - [^ ij-T-2~iI] ++ l -,y21i-i+^il ^j] n+ 2^ n1J-l:in +l - AX in +l,j-i,j] AY [,j-i-l,j (4.6-8) therefore: n+l AT n AT- 2T 2AT UAT VAT n % ij = - a+ i-2ij + (1 y Y Insx L)AX2 AY2 AX AY wheA ri J+l (4.2 a 6-9) where a is a number.

-46For a positive-type equation the coefficients of each term in 5n must be > 0. The coefficient of the'i+lj term is examined first: [AT UAT 1.,I2 AXi If U < 0, this is > 0, and the coefficient is positive. If U > 0, then AT UAT 1 -7 > AX or >U (4.6-10) Since U > 103 is not unlikely, AX would have to be < 0.001 which would require too many grid points. Thus for a forward difference for sX, U must be < 0. i, j-1 term: + VY T4- V > 0 or > V, which will again require many grid points. A general statement can be made concerning the convective term: if U < 0 a forward difference must be used for 6X. If _ X U > 0 a backward difference must be used. (A central difference will not work in either case.) Now assuming U and V have the desired sign, the coefficient of the n -j term becomes: 1 - 2AT - -. UATI V AX2 AY2 AX AY V

-47which must then satisfy the condition: AT + + M + 1 <1(4.6-11) A2 Ay2 AX AY which serves to set the upper limit on AT. Identical conclusions can be obtained by studying the energy equation. Implicit: As before, using forward differences for 5X and backward differences for 6y, the implicit formulation may be rearranged as:.n+l [ 2T 2TAT 1 n GrAT [ n+l n+l AT n+l gi}j I1 + -'2 + --- = ~ij 2X Ti+l j-Ti-lj + i-l 51LX2j AY2J "LX I l i i^j + iXlj rAT T 1 n+1 FAT ATi xi1 n+ AT n+l LAX2 lX2 A+LX ] +1 [+ Ay2 + J 1J - 1+ 2 i +l (4.6-12) where U and V are assumed positive or negative depending on the form of difference used for the convective term. A check of this equation shows that all coefficients are positive, irrespective at AT and hence this formulation is stable for all time steps. (The von Neumann analysis may be used to show that central differences for 5X and 6y should yield a stable system with the implicit formulation.) Again, identical conclusions are obtained from considering the energy equation. ADI: The von Neumann procedure is applied to the ADI formulation to determine its stability. Solutions for n, *, n+l Tn, Tn, Tn+l of the form

-48t3. _,tn ideX iY * -* eiX eijY:,5 -':' e e, Tn n ei eiY T* e,_ iX ei na, n ~ o y Py T P ~ are assumed to exist, and an expression of the form: n+l b. n _e_ J p _c d_ L'ei will be sought. The stability criterion is that X1 and X2, the eigenvalues of the matrix [a b, must be < 1 c d The two steps will be considered separately: Step 1: Substitution of the relations for n t,, and T into the equations of motion and energy, followed by division with the constant factor eiOX eiY yields: Vorticity Equation: pn+l + q - + n - tn (4.6-13) - = pG + s where P = - G i sin OX A- (cos AY-1) q = - (cos wcX -1) tVAT i sin PAY AX 2AY r = T i sin oAX 2AX

-49Therefore: l+s-t ~n + P 1 n+l e = [q+r] [n+ + P 1 ~ (4.6-14) l-q+r 1-q+r Energy equation: 0* - Qn = 9 Q* _ rQ* + s ~n _ ton (4.6-15) Pr Pr or: s s = l-S-t n where Pr (4.6-16) l-q-r q = Pr Second half time step: A similar treatment for the second half of the vorticity equation yields: 2 -n n+l s(S - n+l)- t(n - nl) (4.6-17) Substitution of Equation (4.6-14) for * yields upon rearrangement: iD 2pQn+i [n+l = (l+4s-t)(l+q-r) ~n,4.6-18) (l1-:s+t)(l-q+r) (l-q+r)(l-c+t) Similarly, the second half step of the energy equation yields: 2x* n n+l - (n n+l l (4.6-19) 2 - = s( - ) - t( - ) (4.6-19) Substitution of Equation (4.6-16) for Q* yields: Qn+l (l+ -t l-q-r = /~_~ ~_~t. 8n (4.6-20)

-50Therefore, the overall system may be written as: On+"l a b n A Q [jnl = [a b~] [ - ]f = [A] f (4.6-21) L 9J tc d _e L e where a = (l+s-t)(l+q-r) b 2p (l-q+r)(l-s+t); (l-q+r)(l-s+t) d = (l+ —t)(l+q-r) c = (l-q+r)(1-3 +t) and the matrix [a b is the amplitude matrix whose eigenvalues must cd J be < 1 for a stable solution. Since c = 0 the two eigenvectors are identical to (a) and (d), and IXll = Ial; IX21 = Idl ~ Therefore, the stability criterion requires that lal < I; Idl < 1. By expanding the terms for a and b in the defining relations for p, q, r, s, and, it is possible to prove that both lal and Idl| 1 for all AT, AX, and AY. (See Wilkes (2) p. 67) Therefore it can be concluded that the ADI formulation for the equations of motion and energy should be unconditionally stable. 4.8. Boundary Conditions and Stability. The stability analysis just presented applies only to the interior region of the grid where no boundary points are used in the FDE. Unfortunately, it is often found that a boundary formulation is unstable, while the FDE for the interior points is stable. This problem occurs frequently when higher order approximations are used

- 51for the boundary conditions, or when the boundary conditions are non-linear. There appears to be no general procedure for determining the stability of a given boundary condition formulation. When the boundary condition formulation is unstable spurious values soon appear at the boundary grid points. As these errors become magnified they are carried into the central region and instability occurs. When this problem arises either a different formulation of the boundary conditions, or a shorter time step should be attempted. 4.9. Advantages and Disadvantages of Explicit, Implicit, and Alternating Direction Implicit Finite Difference Formulationso The explicit formulation has the advantage of requiring the fewest calculations per time step. However, stability considerations limit the maximum time step to AT < 10/( 2 2- + i+ -) (4'-1) X- X2 A2 AX AY If, for example, V = U = 100 and AX = AY = 0.1 the maximum AT is AT < 1/600. If, as seems likely, U and V > 100 then AT would have to be smaller yet. These exceedingly small values of AT require an unacceptably large number of time steps to reach steady state. The implicit formulation has none of the stability restrictions of the explicit formulation, but requires the solution of a large set of simultaneous equations in ~n+l and Tn+l for

-52each time step. The solution of these equations requires an iterative technique which can demand an unreasonably large number of iterations. Therefore, it is doubtful that the implicit formulation will offer any advantage over the explicit scheme. The ADI formulation requires neither a limitation on AT nor an unduly large number of calculations per time step, and combines the advantages of both the explicit and implicit formulations while suffering from none of their disadvantages. Therefore the ADI formulation was used. 4.10o Boundary Conditions. The thermal boundary conditions will be considered first: Two types of thermal boundary conditions were used: constant (and known) temperatures on the horizontal surfaces; and insulated surfSraes on the vertical walls. The temperatures are known at all times along the constant temperature boundaries, and hence no further consideration of the case is required, The insulated surfaces are treated as follows: A Taylor's series expansion is assumed for the point next to the boundary: AYX aToj AY2 2Toj Tij = Toj + A +.... (LX 2o- ) +1 aX 2 6X2 From the insulated surface condition, Equation (4 0 ) becomes:

-53a2T 2(Tij-T,j) ax"^4~~~~ = nx2~ ~~~~~~;(4.10-2) x2 ax2 Equation (4.10-2) is then substituted into the FDE of energy applied at the boundary. When these substitutions are made it is found that the FDE of the ADI procedure yield a series of simultaneous equations with tri-diagonal coefficient matrices. These equations can be solved by the recursion formulae presented in Appendix 1. The boundary conditions for the vorticity equation cannot be expressed in a simple form since the no-slip velocity requirements cannot be expressed in terms of vorticity, but must be expressed in terms of stream functions (i.e., i = - = 0 on all surfaces.) Therefore, no actual boundary conditions exist for the vorticity equation. Wilkes (56) and Fromm (12) have suggested the following scheme to circumvent this difficulty: 1) The initial vorticities are calculated throughout the grid using the initial stream functions. 2) The vorticity is assumed to remain constant at the old value Sn along the four boundaries while the internal vorticities are advanced. 3) Once +ij are known a new set of + can be calculated (see next section). 4) From the 1n+l boundary vorticities can be evaluated from Equation (4.10-5),

-545) With the new wall vorticities steps 2-5 may be repeated as often as desired. Although this procedure suffers from the drawback that the boundary vorticities are always one full step behind the interior vorticities, it is hoped that the boundary vorticities will not vary rapidly so this approximation will not introduce excessive errors. If only steady state values are desired, then this approximation should not affect the results since n+l _ n at steady state, The wall vorticity may be calculated from the stream function as follows: _ 72 (4.10-) For simplicity, let us consider the wall at X = 0. The boundary conditions are o = = 0. Therefore 82*=0 and 2 5XV may be calculated from a Taylor's series expansion:, =j = oj + x 2 - 2 (4 lo-4) but toj -~oj 0 Therefore, -2l _x22 (4 10-5) Similar relations can be derived for the other three boundaries, Wilkes (56) suggested the use of a higher order approximation to -2- which included 2j. However, when this author

-55used the higher order approximation the vorticity equation exhibited a boundary induced instability at relatively low AT, and therefore the higher order approximation to 72 was abandoned. 4.11. Determination of en+l When n1 Are Known. The stream function equation is written as: 92 * = _- (4.11-1) or V2*r + = 0 (4.11-2) and in finite difference form: X + t + i = (4.11- 3) with boundary conditions 0 = 0 on all boundaries. Equation (4.11-3) may be written for each grid point not on the boundaries. If the usual three point formulae are used for 2 2 6X and 6y a series of simultaneous equations in ij are obtained. Since there may be as many as four hundred i,j an iterative technique must be used. Several techniques have been proposed for the solution of these equations. In the Gauss-Seidel and Successive Over Relaxation (SOR) methods the FDE are expressed in the form: n+l n r n n n i,j = i,j + W[i-,j - 2Ti,j + i+l,j + AX/(Y ( j-l + 2j + ti j+l)+ AX2 (4.11-4) where W is a constant called the relaxation parameter. For W = 1, n SOR reduces to the Gauss-Seidell iterations, ij is the last approximation to ti j, while ~+1 is the new approximation.

-56The iteration scheme is extremely simple, but suffers from relatively slow convergence, especially if the initial values are not good approximations to *ij. To avoid this difficulty Peaceman and Rachford (39) suggested a different formulation. Equation (4.11-3) is rewritten as: 2X * + &y 2 + = | (4.11-5) That is, the problem is converted to an unsteady problem, and numerical integration is performed in time until the Vi,j field ceases to vary. The ADI formulation is chosen as the most desireable procedure for integrating Equation (4.11-5), and the two half time steps are represented as: 2 2 n AT2 x *i by i'i,j + iij AT 2 ~- tip + Y a + D (4711-6) n+l n b) 1'j-$ij-ij 1 2 b) _AT - 2 bY with boundary conditions r = 0 on all boundaries. Since the boundary values of 4 are known, Equations (4.11-6) need be applied only for the interior regiono Peaceman and Rachford (39) suggested using a varying time step with values of A = 0.25, 0.5, 1.0, 2.0. Each time step was used for a full time step in ascending order of AT o After using AT/AX2 = 2.0, AT was returned to its smallest value and the whole process repeated. A program was written to determine which sequence

-57of AT required the least number of iterations to converge for the problem at hand. For a grid size of AX = AY = 0.1 it was found that a constant value of AT = 0.04-0.06 produced the most rapid convergence. Similarly, for AX = AY = 0.05 a value of AT = 0.04-.06 produced the most rapid convergence. Therefore, AT = 0.04-0.06 gave satisfactory results irrespective of the grid size. It should be noted that this result is not general, and applies only to the problem and boundary conditions considered here. Once a time step size has been established a rapid check of the convergence of Equations (4.11-6) is most desireable. The most common test for the convergence of an iteration scheme involves calculations of the sum c= Z ( n+l-i )2. When e drops below a ~,J certain value the method is considered to have converged. This test,'holTever, has two drawbacks: the evaluation of c requires a rather large amount of useless calculation, and some values of r may vary slowly in a given range and the iteration might be ceased before convergence is attained, The convergence test used in this work avoids both of these disadvantages. During early calculations it was noted that the stream functions near the center approached their steady state values quite rapidly, while the corner stream functions required longer. Therefore, the convergence test was applied only at the grid point i = 1, j = 1 as follows:

-581) At the end of each iteration the vorticity Zl 1 is calculated from the current values of the stream function using Equation 2) The difference (Z 1-t1 1) is calculated. 3) If (Z1 1n+ 1) / 1 i is less than a predetermined constant EPS1, the iteration is said to have converged. If not, a new set of ii j are calculated and the test is applied againo Five to ten iterations were required for convergence of the stream function shortly after new initial conditions had been applied. However, the values of Til n+ canged less radically as the calculations proceeded, and for the major portion of the calculations one iteration was required when EPS1 = OOO1 4o12, Evaluation of the Coefficients U, V, Nu The vorticity and energy equations have been simplified by assuming that U and V are constants over any given time stepo The question -rises whether to evaluate U and V at the old values n n-l of *i 2j, the new values of jni, or some intermediate value of +/2 n n+l n+l n+l/2, Since only the i, j are known before + and T are determined we are forced to use the old values of i in the determination of U and V, unless iteration for the values of +1, Un+1 and Vn+l is performed around each time step. However, since only steady state values are desired, the lag of U and V one n+ l cycle behind: and T should cause no difficulty, and the extra machine time needed to iterate around each time step could not be justifiedo

-59During the solution of PDE it is desirable to use only twoand three-point approximations for 6X or by. However, when evaluating coefficients of the PDE there is no advantage in using these lower order approximations. Thus, in the evaluation of U, V, and Nu it is desirable to use higher order approximations to a ar - aT Ad 1 Ad and d Wilkes (56) studied the effect of using second- or fourthorder approximations on the numerical solution for natural convection between vertical surfaces at different temperatures. By comparing his numerical solutions which used second- and fourth-order approximations to U and V with an analytic solution to the problem Wilkes found significantly better results were obtained with the fourth-order approximation. Therefore, fourth-order central differences were used in determination of U and V. For points not adjacent to a boundary the fourth order approximation to'-X'j can be obtained by writing Taylor's series for the points i+l, i+2, i-l, i-2 in terms of the point at i. ~a2Terms in a 4t a 5 are eliminated from the equations to Ox2 dx5 dx7 8x5 yield: Aij i- 2,oi-r1, j i+lj-i j (4.FDA t ^X "~ 12AX (4.12-1) A third-order DA to may be derived for a point adjacent to the boundary in a similar fashion, and is given by:

-6o-a (-2= (20oj-531,j+6V2,j-3,j )/'6AX (4.12-2) Similar expressions were used for the other boundaries. The mean Nusselt Number Nu has been derived as: 1 H L/H - T --'Nu = 2 I dX 2 L 6y 0 y=O OT where - is evaluated by numerical means. Three different apY=O aT proximations to B of first-, second-, and third-order accuracy were attempted. For Nu < 2.5 the second-order approximation gave the highest values for Nu, while for Nu > 3, the third-order approximation gave the highest values. For Nu < 2.5 the second-and third-order approximations varied by less than 5 per cent, while at Nu > 4, the variation was less than ten per cent. The third-order approximation was chosen as probably yielding more accurate results at high Nu, while at low Nu the difference between the approximations was negligible. The derivative was expressed as: 0Ti o -11Ti o+l8Ti, -9Ti,2+2Ti (4,12-3) 1 Seence of O ons fr Intn of te Vy an 4o13o Sequence of Operations for Integration of the Vorticity and Energy Equations. A complete sketch of the logical sequence in the numerical integration of the vorticity and energy equations can now be presented:

-611) Parameters Pr, Gr, NX, NY, DTAU, and TAUMAX are read into the computer. 2) Initial conditions of T and r are read. 3) The initial velocities are calculated. 4) The initial vorticity field is calculated. 5) The temperature field is advanced over both half time steps. 6) The vorticity field is advanced over both half time steps 7) The new field of stream functions are calculated. 8) The new wall vorticities are calculated. 9) The mean heat transfer rate is determined. 10) If desired, U, V, T, *, and ( fields and Nu are printed. 11) If TAU > TAUMAX the program returns to step 1 and a new set of data are read, otherwise the program returns to step 5, and steps 5-11 repeated.

V. EXPERIMENTAL WORK 5ol. Apparatus In an attempt to prove the assumption of two-dimensional motion, an experimental program was undertaken to determine the mode of natural convection in a long rectangular channel which is heated from below, cooled from above, and sealed on both ends. The channel was formed by placing two parallel, vertical insulating strips between horizontal surfaces as illustrated in Figure 3. Glycerol and air were used as working fluids, with channel dimensions of 1" x 1" x 20". The twenty inch channel length was chosen to minimize the effect of the end walls, while the one inch channel height gave a Ra near Rac for a temperature difference of two degrees Fahrenheit with both glycerol and air at room temperature. During the experiments with air 1" x 1" balsa wood strips were used as side walls, while 1" x 1/4" strips of lucite were used for the glycerine experiments. The upper cooling surface consisted of two 1/2" x 6" x 20" copper plates which were bolted together. A 3/16" semi-circular channel was milled in the lower surface of the upper plate to allow passage of a stream of cooling water between the plates. Holes drilled at the ends of the channel served as an inlet and outlet for the cooling water stream. See Figures 3 and 4 for details. Thermocouples were located within 1/16" of the cooling surface at seven positions on the plate to provide a check on the constancy of the plate temperature. Cooling water at a constant temperature was provided by circulating water from a constant temperature bath. After passage through the - 62

-65TUBE FITTING TO 8 * -LEVELLING SCREWS FROM BATH 1 FOR UPPER PLATES UPPER' " - "\. -,' "\ "j ",.. PLATES' BALSA I WORKING BALSA OOD - REGION lWOOD LOWER PLATES PLYWOOD PLYWOOD LEVELLING PLYWOODLEVELLING SHIM ISHIM TABLE TOP ViEW BALSA WOOD SiDE WALLS LUCIT f SECTION BALSA WOOD i BALSA WOOD Figure 3. Front View of the Experimental Chamber

i__ 20" _2 r+ + |+ "+ WATER OUT 6' WATER IN HOLES FOR LEVELING SCREWS BOTTOM VIEW OF TOP PLATE 3/16" MILLED CHANNEL TUBE FITTINGS Tc TC -SIDE VIEW OF BOTH PLATES/ I''''s't'' J "'"'r`jrlI r I I SIDSOLD C P P E R

-65cooling plates the water was returned to the bath whose temperature was held constant to within + 0.01~F by a mercury temperature controller, and high sensitivity relay. The relay activated a 250W-120V infrared lamp which was operated at approximately 75Vto increase the bulb life. The lower heating surfaces consisted of two 1/2" x 3" x 20" copper plates between which a resistance heater was placed. Heat was supplied by passing an electric current through the resistance heater. Low voltage direct current for the heater was supplied by full wave rectification and filtering of a low voltage AC obtained from a "Powerstat" variable transformer operating on line voltage. A "Solar" voltage regulator was used to reduce major variations in line voltage. In this manner a low ripple, low drift source of power was available for heating of the lower plate. A schematic diagram of the heater power supply is given in Figure 5. Three Co-Cu thermocouples were located within 1/16 inch of the lower plate surface and used to measure the lower plate temperature which was raised or lowered by increasing or decreasing the power supplied to the heater. (See Figure 6.) Heater power was measured by measuring the voltage drop across a standard resistance placed in series with the heater. The thermocouple and resistor voltages were measured with an L and N type 8662 portable potentiometer. Levelling for the lower surface was provided by four adjustable shims which could raise or lower any corner of the lower surface. Once the lower surface was levelled four height adjusting srews were used to level the upper plate.

LINE VOLTAGE 105-120 VOLT 60 CYCLE AC. POWER STATY SOLA VOLTAGE REGULATORE 7120 MH. CHOKE AT 120 MA. 4GE TYPE IN1696 -0SEMI CONDUCTOR 20 /Lfd —' --- ^ —----- - == *- 20 /Xfd — RECTIFIERS RESISTANCE HEATER STD. RESISTOR TO POTENTIOMETER Figure 5. Schematic Diagram of the Heater Power Supply

SOLID COPPER - 20" O O TC O OTC HEATING WIRE NICHROME OR I 2" I' CHROMEL A OSIDE VIEW OF HEATING PLATES TOP VIEW OF BOTTOM PLATE, AND HeEATING WIRE COATS OF INSULATING VARNISH * = i i 1; i RESISTANCES — H1 ------.. I. HEATER SIDE VIEW OF HEATING PLATES Figure 6. The Heating Plates

- 685o2. Discussion of Experimental Observation In the earliest experiments, air was used as the working medium, and flow visualization was accomplished by blowing cigar or cigarette smoke into the working region. Lighting was provided by shining a 60 W incadescent bulb through a lucite window which had replaced a small portion of one side wall. The low intensity of the illumination, coupled with the high smoke density, made observation of the smoke particle paths difficult. To increase the contrast between the smoke particles and their surroundings, the side walls and copper surfaces of the channel were coated with flat black actylic laquer and the experiments repeated. However, the low intensity of illumination still prevented visual observation of the flow pattern. To increase the strength of illumination, the 60 W bulb was replaced by a 100 W and then 200 W bulb.- With the higher illumination produced by the 200 W bulb it was possible to distinguish the flow patterns. However, with the larger bulbs a distinct change in the convection pattern was observed to occur shortly after the bulb was turned on. At first the particles in a particular plane along the Z(20") axis appeared to be rising or falling. This necessitates a return flow at some other plane, and hence a threedimensional, cellular pattern. However, this pattern soon changed into a planar rotation with fluid always rising on the wall opposite the source of illumination. Apparently the radiant energy from the bulb had heated the black wall sufficiently to produce a two dimensional flow. The 100 W bulb produced the same effect.

-69A section of the black wall opposite the illuminating slit was replaced by a section of transparent lucite in order to reduce the absorption of radiant energy on that wall. Although absorption of radiation on the wall opposite the illumination was no longer a serious problem, the illumination slit itself absorbed enough radiant energy to cause a two-dimensional flow similar to that which had occured on the black wall. After these experiments with air, glycerine was used as the working fluid. Flow visualization was accomplished by dropping aluminum dust particles through a small hole in the cooling plates into the working fluid. By watching the traces formed by the Al dust particles it was possible to determine the flow patterns of the glycerine. To minimize the effect of conduction from the surroundings into or out of the working region, the cooling plate temperature was set at 77~F, and the whole apparatus covered with two layers of rock-wool insulation. The first experiments were run at temperature differences of 10-12~F (corresponding to Ra 5-6000) between the upper and lower plates. The flow patterns detected in these early experiments exhibited a three-dimensional pattern as illustrated in Figure 7. Al dust introducedthere QIY ) S / Front view Side view Figure 7. Sketch of Flow Patterns Seen in Liquid Glycerine at Ra e 3000.

-70In the front view of Figure 7 a series of darker and lighter regions were noticed. These regions apparently were the result of variations in the refractive index of the glycerine caused by temperature gradients in the liquid. The general shape of the refractive index lines appears to indicate the presence of two "cellular regions" in the X-Y plane. This conclusion was substantiated by the random presence of the Al dust particles only in the left "cell", indicating that no transfer of fluid occurred between the left and right "cells", but did occur within each "cell". These cellular regions appeared to be the result of secondary flows caused by a raising or lowering of the side wall temperature because of conduction into or out of the side walls. In the side view of Figure 7 a cellular pattern was observed in the Z direction. The vertical lines between the cells were well defined by Al dust particles that were trapped between the two cells, but could also be seen as a refractive index line when no Al dust was present. This three-dimensional pattern appeared to be similar in nature to the pattern observed with air in the short instant before radiation initiated the two-dimensional flow. In most cases the cells were approximately 1-1/4 to 1-1/2 times as long (Z direction) as they were high. As the temperature difference between the plates was lowered, the variation in room temperature (as much as 15~F between 9:00AM and 3:00PM) prevented the attainment of steady state since a heater power which caused a 2~F temperature difference in the morning might cause

-71an 8~F difference in the afternoon. At times a temperature difference was established by conduction from the room even though no power was supplied to the heater. Because of the extremely slow convective motion (on the order of one to two millimeters/minute) there was no chance to obtain any steady, or pseudo-steady, state results for the low temperature differences. 5.3. Conclusions and Recommendations The numerical calculations performed in this work have been based on the assumption of two-dimensional flow. If this assumption can be proven then there will be little doubt that the numerical calculations are a true representation of the natural convection within a long rectangular channel which is heated from below and cooled from above. Therefore, the preceding experiments were intended to yield a qualitative description of the natural convection which existed in the channel. The equipment was not expected to provide precise quantitative heat transfer rates, or temperature profiles (the numerical calculations would provide these), and therefore a minimum of temperature control and regulation was provided. It turned out that the stable, or preferred, mode was extremely sensitive to minor perturbations in the boundary conditions such as the illumination, heat leakage, etc. The equipment thus did not prove to have adequate control for the desired experiments. Although interesting observations were made, they did not provide the hoped for critical test of the two-dimensional assumption.

-72From the results of the experimental work the following recommendations are made for future experimentation: 1) The equipment should be built in a room with close temperature control to eliminate the problem of varying ambient temperature. 2) Good thermal control of all boundaries is necessary. 3) The working fluid should allow a critical temperature difference of about 5~F for a two inch cell height. The 5~F AG will prevent minor variations in the temperature of the heating or cooling plate from causing major secondary flows, while the two inch cell height would allow simpler and more precise optical visualization.

VI. RESULTS 6.1. Introduction A computer program for the numerical integration of the equations of vorticity and energy for natural convection in an enclosed rectangle which is heated from below, and cooled from above, was developed using the "MAD" computer language. The programmed calculations were performed on the IBM 7090 digital computer at the University of Michigan Computing Center. Approximately sixty different combinations of Gr, Pr, and L/H were studied. In each case a set of arbitrary initial conditions were fed to the computer and calculations were performed until a prechosen number of time steps was exceeded. The current values of T and 4 were then punched on cards and the calculations halted. A check of the behavior of the mean Nusselt number served to indicate whether or not the calculations had approached steady state. If the Nu was still changing rapidly with time the punched values of T and were used as input conditions for the next run and the procedure was repeated. In almost all cases two such returns were sufficient to yield nearly steady state behavior. (For high Pr near Rac more calculations were required unless the initial conditions were very close to the steady state values.) In the majority of the cases steady state was reached on the initial sequence of calculations. The condition of no motion and a linear temperature gradient is a valid solution to the equations of vorticity and energy for all -73

-74Gr, Pr, and L/H. Therefore, if a linear temperature gradient with no motion were used as initial conditions, no motion would be expected to develop until roundoff error had produced a sufficient disturbance to initiate a motion, which would then develop to its steady state. It seems plausible that sufficient roundoff error to initiate motion might not occur for a considerable amount of computer time, and therefore it was desirable to supply some form of disturbance in the initial conditions. By doing so, the unsteady calculations lose their significance except when applied to the specific disturbance used to initiate the flow. Since this investigation was primarily interested in the steady state flow patterns, this was not a serious drawback, and computational effort was minimized by using a distrubance which closely resembled the steady state convection pattern, or in most cases the final results of the calculation with the closest values of Pr, Gr, and L/H. A summary of the cases studied is presented below: TABLE I SUMMARY OF CASES STUDIED L/H Pr Gr 1/2 1.0 5000., 5000., 6000., 7000., 8000., 10000., 20000., 1.0 0.01 150000., 170000., 200000., 300000., 500000., 1000000., 2000000., 4000000 1.0 0.03 50000., 56670., 66670., 100000., 167000., 333000., 666700.

-75TABLE I CONT'D SUMMARY OF CASES STUDIED L/H Pr Gr 1/1 0.1 13000., 14000., 15000., 17000., 200000., 30000., 50000., 100000., 200000. 1/1 1.0 1300., 1400., 1500., 1700., 2000., 3000., 5000., 10000., 20000., 40000. 1/1 5.8 180., 230., 260., 360., 500., 860., 1720. 1/1 25.0 60., 75., 100., 200., 400. 2/1 1.0 1000., 1100., 1300., 1500., 2000., 3000., 5000., 10,000. 3/1 1.0 3000. For each case a plot of Nu versus time was prepared, and in those cases where steady state had been closely approached, but not attained, the Nu curve was extrapolated to obtain a reasonable estimate of the steady state Nu. In no cases did the uncertainty of this extrapolation exceed ~ 0.2 when Nu > 4.0. For Nu < 4.0 the uncertainty was in the order of 0.03 In addition to steady state Nu, plots of the streamlines and isotherms which existed after the last time step were prepared by an interpolating and plotting program. (See Wilkes (56) p. 93.) These plots were prepared with approximately equal increments of T and j to give a rapid estimate of the velocity and heat flux magnitudes. The program uses a linear interpolation to locate the position of an isotherm or streamline that does not fall directly

-76on a grid line. The coordinates of the various points on the streamline or isotherm are stored in column vectors, and plotted by the University of Michigan Executive System plotting subroutines. Although it is impractical to present every isotherm and streamline plot, an effort has been made to present those plots which best represent the various flow regimes encountered in this work. Plots of the steady state Nu versus Ra for various L/H and Pr are presented and compared with the numerical calculations of Aziz (1) and the experimental work of Schmidt and Silverton (50). 6.2. Discussion of Results a. Grid Size and Numerical Stability. To check the effect of grid size on the steady state results, grid sizes of AX = AY = 0.05, 0.10, and 0.20 were used with a square cavity for Gr = 3000., and Pr = 1.0. The results of these calculations are presented below: TABLE II EFFECT OF GRID SIZE ON COMPUTATIONS IN A SQUARE CAVITY FOR Gr 3000., AND Pr 1.0. X Nu I cent Computing time, min. 0.05 lo84 4.72 25. 0.10 1.92 4.92 6. 0.20 2,77 6. 1.

-77By extrapolating a plot of Nu and 4 cent versus AX to AX = 0 (see Figure 8) the values of Nu = 1.81 and 4 cent = 4.60 were determined and represent the estimated solution of the differential equations. The values of Nu and 4 cent for AX = 0.05 agreed within 2 per cent with the extrapolated values for AX = 0, while the values for AX = 0.1 agreed within 6 per cent with those extrapolated for AX = 0. It is noted that the computer time increases almost as the cube of number of grid spaces, and therefore AX = 0.10 was chosen as a reasonable compromise between accuracy and computer time. All figures and results are based on calculations for AX = 0.10. In the early calculations for Pr = 1.0, AT = 0.005 was used. However, at higher Pr the steady state was approached more slowly, and an attempt was made to find the maximum stable AT. It was found that AT > 0.01 produced a boundary-generated instability. Therefore, AT = 0.01 was used for the remainder of the high Pr cases. For lower Pr where steady state was approached more rapidly a lower AT was used in most cases. For those runs where Gr = 1000000. (Pr = 0.01) it was necessary to restrict AT to AT < 0.0002. For Gr = 4000000, AT was restricted to AT < 0.0001. In the high Gr, low Pr, cases the steady state was approached very rapidly so that 100-150 time steps were needed to reach steady state.

-783.0 7 2.8 6 2.6 -5 2.4 -4 Nu cent. 2.2- -3 2.0 / -2 Nu 1.8- -I I I I I I 0 0 0.05 0.10 0.15 0.20 0.25 0.3 GRID SIZE Figure 8. Effect of Grid Size on Steady State Nu and cent for a Square Cavity with Pr = 1, Gr = 3000

-79In general, 50 to 150 time steps were required to reach steady state in any given calculation. The program performed 100 time steps and the associated print out in approximately 2 to 2 1/2 minutes. For cases in which the starting values of T and 4 were poor approximations to their steady state values, a non-numerical instability in which the Nu oscillated in a bounded, but non-converging fashion was encountered. Decreasing the time step size had no effect on the oscillations. Conceivably the oscillations might have dampened out if the calculations had been continued. This was not attempted, however, since the difficulty could be circumvented by using better starting values for ~ and T b. Nusselt Number Variation With L/H, Gr, and Pr Both Kursweg (24), and Pellow and Southwell (40) found that the dedimensionalized, linearized equations of motion could be reduced to the form Nu = f(Ra). That is, Nu was found to be a function only of Ra, and not an independent function of Gr and Pr. Morgan and Warner (33) indicate that linearization is justifiable for all fluids with high Pr. To check the validity of Morgan and Warner's conclusion, and to define the range of Pr described as "large," a plot of Nu vs Ra for L/H = 1. and Pr = 0.01, 0.03, 0.1, 1,0, 5.8, and 25. is presented in Figure 9. It is apparent from Figure 9 that for Pr > 1. the assumption of linearization is well justified except possibly in the region of the critical Ra, where minute variations between Pr = 1 and Pr = 5.8 or 25.0 may be noted. Therefore,

go I t 1 I I I, 64z' 0 LEGEND 2-1 PR 5.8,25.0 0 PR - 1.0 0 PR a=0.1 A PR -0.03 0 PR a 0.01 10 2 4 6 8 10 2 4 6 8 1 CpRa gOAH3 Ro = fa yr — k 2v Figure 9. Dependence of N- on Ra and Pr for a Cavity of Square Cross Section

-81for all fluids with Pr > 1.0, Nu = Nu(Ra). Since almost all liquids and gases (except for liquid metals) meet the condition Pr > 1.0, this is a very widely applicable result. For comparison, the results of Aziz's (1) numerical calculations for two-dimensional natural convection within a square region heated from below and cooled from above have been presented in Figure 10. It should be noted that Aziz's calculations yielded significantly higher heat transfer rates, and a greater variation of Nu with Pr, than the present calculations. When Aziz's results for natural convection in an enclosed cavity were compared with the experimental data of Schmidt and Silveston (50), it was found that the critical Ra extrapolated from Aziz's calculations was lower than the experimental value for convection between two infinite plates. This situation does not appear to be physically realistic, since one would expect the introduction of sidewalls to increase rather than decrease the critical Ra because of the increased viscous drag. This apparent inconsistency in Aziz's calculations does not occur in the results present herein, or the analytical solution for Rac presented by Kurzweg. The dimensionless quantity (Pr *) as well as Nu, was found to be only a function of Ra when Pr > 1. Since * is the integral fUdY or fVdX, (Pr ~) can be expressed in terms of dimensional quantities as: Pr = c Py udy = $(Ra) k y=o

1 I |l 64<1 pr=?7. I,^ NUMERICAL CALCULATIONS 2 OF AZIZ FOR A SQUARE CAVITY NUMERICAL CALCULATIONS OF THIS WORK FOR A SQUARE CAVITY Pr = 1.0 XPERIMENTAL DATA OF SCHMIDT SILVESTON FOR INFINITE HORIZONTEL PLATE Pr2I.O 3 4 6 000 2 46 8 00 Cp/. 9g:^A H3 k 2v Figure 10. Comparison of Results with the Calculations of Aziz and the Data of Schmidt and Silveston

- 83That is, u is independent of pi for a given Ra, and the interesting conclusion is reached that all fluids with Pr > 1, and equal thermal diffusivities will have identical dimensional velocities independent of physical properties and cell size for the same Ra. For those cases with Pr < 0.1 the Nusselt number exhibits a noticeable dependence on Pr as well as on Ra. At low Ra, a decrease in Pr produces a decrease in Nu. However, as Ra increases, this monotonic behavior ceases and decreasing Pr first causes a decrease and then an increase in Nu. In order to check the consistency of this result against the total convective circulaH/2 tion in each case, the dimensionless circulation ( *t Pr = PrH/ UdY) cent o was calculated for various Pr at Ra = 10,000, and the results are presented below: TABLE III VARIATION OF CENTRAL (J*Pr) AND Nu WITH Pr at Ra = 10,000, AND L/H = 1.0. Pr Nu,*Pr 0.01 2.92 10.3 0.03 2.87 10.25 0.1 2.76 10.0 1.0 3.025 10.8 5.8 3.03 10.8 25.0 3.03 10.8

-84From these results it is seen that the minimum in heat transfer rate corresponds to a minimum in the circulation rate. When Nu is plotted against Gr at the various Pr studied, as shown in Figure 11, a series of monotonic functions which increase with increasing Gr is obtained. c. Determination of the Critical Ra, and Comparison with Theoretical Solution of Kurzweg (24) and Numerical Solution of Aziz (1) In keeping with commonly accepted practice, critical conditions for the onset of convection will be stated in terms of Ra, L/H and Pr. The critical Rayleigh number, Rac, is defined as that Ra above which convective motion is stable, and below which conduction is stable. The critical Rayleigh number was determined herein by extrapolating the curve of Nu vs Ra until it intersected the line of Nu = 1, representing pure conduction.(see Figure 12) At this point a decrease in Ra no longer causes a decrease in Nu. The results are presented in Table IVo In 1965, Kurzweg (24) published a study on the stability of fluid in an enclosed rectangle heated from below. He used the linearized equations of energy and motion to predict the critical Ra as a function of L/H, but could not predict the variation of Rac with Pr. However, a comparison of Kurzweg's results and those of this work indicates good agreement for Pr >1. For Pr < 1 a significant variation of Ra with Pr was found which indicated that

I I I I I I! 1 / I! I J X I I I I 1 1 6.0 5.0 C0 c 4.0c ~1 Prs 1.0 Pr'O.I — r-O —0 Pr-0.01 z 3.0 2.0 1000 2 4 6 8 104 2 4 6 8 10 2 4 68 10 2 4 6 Gr gAflHs 2U2 Figure 11. Dependence of Nu on Gr and Pr for a Cavity of Square Cross Section

8 -' 645 2 - 1 - - r LEGEND I~~~~~~~~~~ I~ ~ ~ ~ ~~~Je --— O —L/NO5,Pr:I.o 2^ ^^^^^ /^^~~~~~~~, -— O —L/HI.O,Pral.O ^y^^ ^^ - -.~ —.L/Hu2.OPrI.O 00' ^^^ ^^ —-- EXPERIMENTAL DATA OF.0e /SCHMIDT a SILVESTON FOR /^.^P~~~~~~/ ~CONVECTION BETWEEN IN/^' P'~~~~~~~~~~~~~~~~~ ^ /RFINITE PLATES HEATED 7' f FROM BELOW 8 9 10o 2 4 6 8 10 2 4 6 B 10 cp,~L 90,&8H Ro= gAH3 k 2v Figure 12. Effect of L/H and Ra on Nu for Pr = 1.0

-87the assumptions needed to linearize the equations of motion and energy are not valid. The variation in Rac does not become great until Pr < 0.1, and would probably become quite large for Pr < 0.001. TABLE IV CRITICAL Ra AS A FUNCTION OF L/H AND Pr L/H Pr Rac Kurzweg Aziz's Rac Rac 1/2 1,0 920 1008 1. 0.01 1500 1. 0.03 1380 1. 0.10 1270 1. 1.0 1240 1290 850 1. 5.8 1230 1. 7.0 750 1. 25.0 1230 2.0 1.0 5800 6060 For comparison, the values of Rac obtained by extrapolating Aziz's (1) results have been included in Table IV. The inconsistency with the values of Kurzweg and of this investigation are apparent. In 1961 Chandrasekhar (6) presented a work by Landau in which an attempt was made to predict the "strength of convection" as a function of (Ra-Rac). Landau's results were developed by combining the various solutions of the linearized stability problem in such a way as to satisfy the non-linearized equations of motion. From his

-88results, Landau predicted that the "strength of convection" should increase in proportion to (Ra-Rac)l/2. In an effort to check the validity of this result, (Nu-l), and total circulation *cent' (hopefully both measures of the "strength of circulation") were plotted vs (Ra-Rac) on logarithmic coordinates. A sample of these plots is presented in Figure 13. It was found that Vcent increases as (RaRac)l/2 over the whole range of the investigation, while (Nu-l) increased in proportion to (Ra-Rac)l for Ra; Rac, and then in proportion to (Ra-Rac)l/2 as Ra > Rac. Thus, the prediction of Landau that the strength of convection should increase in proportion to (Ra-Rac)1/2 appears to be valid, even in regions where the validity of the assumptions used in its derivation is doubtful. d. Streamline and Isotherm Plots: In Figures 14 to 24, streamlines and isotherms are plotted for various L/H, Pr, and Gr. In Figures 14 to 18 the build-up of convective strength can be observed as the Ra increases in a square cavity with Pr = 1. At low Ra the isotherms are almost straight, equally spaced, horizontal lines, while the streamlines are nearly circular. As Ra increases the isotherms become increasingly distorted as a nearly isothermal core with -0.2 < T > +0.2 developso In this core region conduction is almost negligible and convection produces the major portion of the heat transfer. As Ra increases, the size of the core region increases. The streamlines, on the other hand, have distorted slightly from their circular shape into

8~ 8 ~~~~~~6-.-i -S~~~~~~~~6 4 4 ~~~~~~~~~~~~^^2 2Slope=O.5 (Nu-I) - o,) 8- /1 } ^Slope=0.5 46- -Slope=1.0 2100 2 S8ope =O 100 2 4 6 8 1000 2 4 6 8 10,000 2 4 6 8 ( Ro -Rc) Figure 13. Dependence of Strength of Convection on (Ra-Rac) for L/H = 1.0 and Pr = 1.0

0.0 X 0.25 0.50 0.75 1.0 0.000 TX-I.0 y _s T00. T= -O. 4_4 0.a50nd Gr 00 T=-0. 2 T~0.4 0.750 =0.6 O.500 GT= 1.0..750 an G 10

-910oo X 0.25 0.50 0.75 1.0 0.00s-.0 0.2O50 — 7 —---—' O50T=- 0.60.,!/ / / 0.50. T.. i0O 25 --- ^ ^ \-\T=I.O 1.000 0.70 X 0.25 0. 50 0.75 1.0 * —0.5 i 1 =-I.07 0.2500 0Figure 15. Steady State Streamlines and Isotherms for L/H = 1.0, Pr = 1.0, and Gr = 2000 0.75

-920.0 X 025 0.50 0.75 1.0 0.250 T= —.6 -O.2 T:O. 4. 500 — __ Figure 16. Steady State Streamlines and Isotherms for L/H = 1.0 1.000 Pr = 1.0, and Gr = 5000 0.000 x 0.25 0.50 0075 1.0 i /=o Yl *^^ ^-2.4 00000 —^ ~ *=-3.6 ~= -7.2 0.500 [\ —., 0 750 il/=o 1.000 Figure 16. Steady State Streamlines and Isotherms for L/H = 1.0, Pr = 1.0, and Gr = 5000

-950.000 0200 050- -- — X 025. 050 075 1.0 o 500 T. _. \0^ 250T=04 T —002 \ 1.000 Pr = 1.0, and Gr 00 X 0 25 0 50 0 75 1.0 0.000 i=-1 8 q//=-3.6 =-7.2 0 500 \CI=o 1000 — -------— =0 —----- Figure 17. Steady State Streamlines and Isotherms for L/E = 1.0, Pr = 1.0, and Gr = 10000

-9)t0.o0 o 0 x 0.25 0.50 0 75 1.0 O.5000 ~~~~~~~~To.Q T-I. 00 — T=-O. 0.7 50 /0. I Pr = 1.0, and Gr = 20000 0.500 -~ T=-0.2*-2. 1 0050 -------- — ~ -------- w - 0.250^ ^ - 0 =-10.2 ~ ~ ~ ~ T=. 0.500 ii/=0 -^s^ v|/=0 I. OOO1 —----- — t! — ---------------- Figur l8 taySaeSralnsan stem o /. Pr = 1.0, and G = 2000

-95an oblong shape, but their main character has remained unchanged. It may be noted in several cases that a region of separated flow has developed at the two blunt corners of the oblong. Unfortunately, the coarseness of the numerical grid has prevented any quantitative description of this separated region which does not, however, appear to have any marked effect on the heat transfer rate. The separated region exists only in certain restricted ranges of Ra which vary with Pr. The existence or non-existence of a separated region was not a function of starting conditions, since non-separated initial conditions at times produced a separated final solution and vice versa. In Figures 19 to 22 streamlines and isotherms are plotted for a cavity with L/H = 0.5, and a fluid with Pr = 1.0. At low Ra the isotherms are nearly straight horizontal lines, while the streamlines are oblong and nearly vertical over a large region. As Ra increases, the isotherms become more distorted until at high Ra the isotherms in the central region are almost vertical and parallel; the streamlines in this region are also nearly vertical. This suggests a possible form of thermal boundary layer in which all resistance to heat transfer is confined to a thin layer along the upper and lower surfaces while convection transfers heat between the two layers. The appearance of this thermal boundary layer could be expected to decrease the rate of increase of Nu with Ra since increasing Nu can be accomplished only by decreasing the thickness of the thermal boundary layer. Figure 12 shows that Nu increases quite slowly with increasing Ra for Ra > 10,000. From this postulated

0.0 X 0125 0.250.375 050 0.0 X 0.125 0.250 0.375 0.50 0. 0.00( TI ^~.-X0.08 To-O. r~r~-0.160.250 0.250. Ton E =. =//* 0 Figure 19. Steady State Streamlines and Isotherms for L/H = 0.5, Pr = 1.0, nd Gr = 0.500' o0.50 ~~~~~~~~~~~~~~T0.70...0 ~~~T- 0.8 IO —601.000 ------------- T- 0.80 Figure 19. Steady State Streamlines and Isotherms for L/H = 0.5, Pr = 1.0, and Gr = c0

0. 0.125 0.250 O.375 0.500 0. 0 X 0.125 0.250 0.375 0.500 0. 0X 0 1 5020075000 0.000 - o~ooo - ~lO i I ~.,oT Too --- - -^ry-1 — 0 j,.o Ts~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - T-0.__^ ~~~__..___.___.. "-0.5761.. 0.125 0.25 T.-(~6j " — - - 0.250 k 0.3 0.15 /,/ 1 // 0,315 ^/ / ^^u /^ \\\\\ Tlr0.2 ~ ~ ~ ~~~~. ~0.2~~~~~.65 J::l^ / / 2i0'V // T 0._ TL 0.4,.~, -- -/-/ —T^/~ o..s U\v J/"~ 0. ------------- -----— 75 --- \\ — \ ^ y / / - O. ~875 —-" — - 0~ ~ O. ~-875 —= 0. ~ T=1.0 0I~,.000ooo000 Figure 20. Steady State Streamlines and Isotherms for L/H =0.5, Pr = 1.0, and Gr = 000

0.0 X 0.125 0.250 0.375 0.50 0.0 X 0.125 0.250 0.375 0.500 0.0 —T " ------------- 0.0 To- 1.0 Y I - I-~-~ -- 3~9 0.125 ----—. 1 1 _ ___ 0.125 0.250 --- ----- 0 250 9.000 __ __ __ *7 ___ _ Fgr20.375 0or37 L 0 P 00-11To20 8.0 O 0-500 ---------- 0~500 T-20 0.~250' --------------- / - --- / --- r --- 0. 25 --- 40 - 000 O.87S' - -5 _. - 0. 75 1.000 -T,.000 Figure 21. Steady State Streamlines and Isotherms for L/H =0.5,Pr = 1.0, and Gr = 10000

0.0 X 0 125 0 250 0.375 0.50 00 X 0.125 0 250 0.375 0.50 T. - 10______________ T o - - 0. 5 0 O 1.000 0 0 Figure 22. Steady State Streamlines and Isotherms for L/H 0.5 Pr 1.0. and r 20000-1.3 0.125 --- - -- r 0 —----- ~125 T-0.6 / f 0.2 ~.250 — 3 //.4 V \ \ ----- 12 -002\ 0.3. 0 \Oi 0.875 - 0.8 0. 00 -— 0 0 0 __0 Figsre 20 O.T50 O.T50 - Ts 0 T- 0. 8 Figcre 22. Steady State Streamlines and Isotherms for L/H 0.5 Pr 1. O., and Gr 00

-100thermal boundary layer one might expect the heat transfer rate to be independent of L(for L/H. 0.5) for those Ra at which the thermal boundary layer has developed, provided that L is much greater than the thickness of the viscous boundary layers on the side walls. However, no numerical, or experimental, data are available to confirm this hypoties is. In Figures 23 and 24 streamline and isotherm patterns are presented for L/H = 2.0 and 3.0, with Pr = 1.0. The existence of multi-celled patterns in these plots should be noted. The problem of determining the most stable number of cells for any L/H is rather complex, and will be discussed in the next section. It is sufficient to say that the most stable number of cells for a given cavity was found to be equal to L/H, for the integral values of L/H used in this study. The effect of Ra was studied in a cavity with L/H = 2.0 and Pr = 1.0. The Nu for these cases are presented in Figure 12, an examination of which indicates the similarity, in the critical region, between the two dimensional calculations for L/H = 2.0, and the experimental work of Schmidt and Silveston (50). This close agreement may indicate the presence of two-dimensional rolls in the experimental work for Ra ~ Rac. The divergence of the two curves for Ra ~ Rac may be traced to two factors: 1) Two-dimensional rolls have been replaced by cells in the experimental work, or 2) The effect of coarse grid spacing has become significant as the Ra is increased.

o0o X 025 0 50 0,75 1.00 1.25 1.50 1.75 2.0 0. T Fiue2j.SedySaeStemiesadIothrsfr/1.0 r n G 00 0.250 010 To~. To. 4.2 0 T-0.6 Two Too....07501 - ^ j__ -- T-0.8 Too Figur -- e ~-"-T-It 0 0 0 0X 0 25 0 5 0 075L 1.25 I 50 I 75 2.0 oo, ~-9 02 -- -— + —----- -.o_, —9' i 0 50 -. —--.o Figure 23. Steady State Streamlines and Isotherms for L/H = 2, Pr = 1, and G: = 10000

0.0 0.25 0.50 0.75 1.0 1.25 1.50 1.75 2.0 2.25 250 275 3. I.0 0 To T- 1.0 0~'""0- ~~- —.-.~_~0^=== ~ — ~~ ~ -..0.0 X 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 2.50 2.75 3.00 0.o I^-,.o 25] _o 0.50 0 "-.2 To-.4 To- TooO y-~o ayo i *o 1.0 0_-_____- 0_-______ _______ _ 1 Figure 24. Steady State Streamlines and. Isotherms for L/H == 2 Pr = 1 and. Gr = 10000 F igure 24. Steady State Streamlines and Isotherms for L/H 2 Pr 1 and Gr 10000

-103e. Preferred Mode of Natural Convection: In Figures 14 to 24, it is noted that the number of convective cells is equal to L/H, for the integral values of L/H used in this study. That these are not the only convective patterns which can exist is clearly indicated by the work of Fromm(13) who found with L/H = 7/3 that two-, three-, or four-celled patterns could exist, depending on the type of symmetrical disturbance introduced as initial conditions. This was a highly artificial initial condition, since any physical problem would probably be subjected to some form of unsymmetrical initial conditions, or unsymmetrical disturbance. Therefore, it was desirable to determine which cellular configuration was the preferred mode; that is, what configuration would result if the one-, two-, or three-celled patterns were subjected to non-symmetrical initial conditions or to non-symmetrical disturbances. To test the effect of initial conditions an unsymmetrical two-celled configuration was used as a starting pattern in a cavity with L/H = 1. In a short time the stronger of the two cells engulfed the weaker cell, and the single celled pattern proceeded to the same steady state conditions as were found when a onecelled pattern was used as a starting condition. In the second test, the same two celled unsymmetrical pattern was used in a rectangular cavity with L/H = 2.0. This time the two cells became symmetrical, and the final steady state was identical to that obtained when two symmetrical cells were used as a starting condition.

-104In the last test of the effect of initial conditions a small one-celled disturbance was placed in the upper left hand corner of a rectangular cavity with L/H = 3.0. As time progressed the one cell became two, three, and then four cells. After a time the fourth cell was engulfed and the final pattern consisted of three cells which were symmetrical about the center line of the cavity, as illustrated in Figure 24. Finally, one-, two-, and three-celled patterns were generated in a rectangular cavity with L/H = 2.0 by use of symmetrical initial conditions. The three patterns were then subjected to a non-symmetrica. disturbance of the temperature field. The one-celled pattern soon split into two cells. The two-celled pattern simply returned to its symmetrical condition, and the three-celled pattern was soon converted into a two celled arrangement. In all cases the final steady state conditions were identical to those which occurred when a symmetrical, or non-symmetrical, two-celled pattern was used as the starting condition. These tests are not absolutely conclusive, since little effort was made to study the effect of various disturbances, or to determine the minimum disturbance that will cause the shift from one configuration to another. However, they do suggest that, although many cellular arrangements can be forced to occur in any given cavity, only that arrangement in which the number of cells are equal to L/H (for an integral value of L/H) will be stable to a general non-symmetrical disturbance, and therefore is the preferred mode of convection.

BIBLIOGRAPHY 1. Aziz, K., A Numerical Study of Cellular Convection, Ph.D. Thesis, Rice University, Houston, Texas, 1965. 2. Barakat, H. Transient Natural Convection Flows in Closed Containers, Ph.D. Thesis, University of Michigan, Ann Arbor, Mich., 1965. 3. Benard, H., "Les Tourbillons Cellulaires," Rev. Gen. Sci. Pur. et Appl., 11, (1900), 1261-71 and 1309-28. 4. Brian, P.L.T., "A Finite Difference Method of High-Order Accuracy for the Solution of Three-Dimensional Transient Heat Conduction," A.I.Ch.E. Journal, 7, (1961), 367-370. 5. Chandra, K., "Instability of Fluids Heated from Below," Proc. Roy. Soc. London, A 164, (1938), 231-242. 6. Chandrasekhar, S., Hydrodynamic and Hydromagnetic Instability, Oxford University Press, (1961), 612-615. 7. Deardorff, J., "A Numerical Study of Two-Dimensional Parallel Plate Convection," J. of Atm. Sci., 21, (1964), 149. 8. Douglas, J., "On the Numerical Integration of d2u/ x2 + ~2u/y2 = au/6T by Implicit Methods," J. of the Soc. of Ind. and Appl. Math., 13, No. 1, (1955), 42-65. 9. Douglas, J., and Peaceman, D.W., "Numerical Solution of Two-Dimensional Heat Flow Problems," A.I.Ch.E. Journal, 1, (1955), 505-512. 10. Douglas, J., and Rachford, H.H., "On the Numerical Solution of Heat Conduction Problems in Two and Three Space Variables," Trans. Am. Math. Soc., 82, (1956) 421-439. 11. Forsythe, G. and Wasow, W., Finite Difference Methodes for Partial Differential Equations, John Wiley and Sons, Inc., New York, 1964. 12. Fromm, J., The Time Dependent Flow of an Incompressible Viscous Fluid, Los Alamos Scientific Lab. Report, No. La-2910, 1963. 13. Fromm, J., "Numerical Solutions of the Nonlinear Equations for a Heated Fluid Layer," Physics of Fluids, 8, (1965), 1757-1769. -105

-10614. Globe, S., and Dropkin, D., "Natural Convection Heat Transfer in Liquids Confined by Two Horizontal Plates and Heated from Below," J. of Heat Transfer, 81C, (1959), 24. 15. Hellums, J.D., Finite Difference Calculation of Natural Convection Heat Transfer, Ph.D. Thesis, University of Michigan, Ann Arbor, Mich., 1960. 16. Hellums, J.D., and Churchill, S.W., "Dimensional Analysis in Natural Convection," CEP Symposium Series, 57, No. 32, (1961), 75-80. 17. Herring, J.R., "Investigation of Problems in Thermal Convection — Rigid Boundaries," J. Atm. Sci., 22, (1964), 277-290. 18. Herring, J.R., "Investigation of Problems in Thermal Convection," ibid., 21,(1963), 325-338. 19. Hollands, K.G.T., "Convectional Heat Transfer Between Rigid Horizontal Boundaries after Instability," Phys. of Fluids, 8, (1965), 389-390. 20. Jakob, M., "Free Heat Convection Through Enclosed Plane Gas Layers," Trans. ASME, 68, (1964), 189-194. 21. Jeffreys, H., "The Stability of a Layer of Fluid Heated from Below," Phil. Mag., 2, (1926), 833. 22. Jeffreys, H., "Some Cases of Instability in Fluid Motion," Proc. Roy Soc. London, A 118, (1928), 195. 23. Kuo, H., "Solutions of the Nonlinear Equations of Cellular Convection and Heat Transport," J. of Fluid Mech., 10, (1961), 611-634. 24. Kurzweg, U., "Convective Instability of a Hydromagnetic Fluid within an Enclosed Rectangular Cavity," Int. J. of Heat and Mass Trans., 8, (1965), 35-41. 25. Lax, P., and Richtmeyer, R., "A Survey of the Stability of Linear Finite Difference Equations," Communications in Pure and Appl. Math., 9, (1956), 267. 26. Low, A., "On the Criterion for Stability of Viscous Fluids Heated from Below," Proc. Roy. Soc. London, A 125, (1929), 180. 27. Low, A., Proc. Third Int. Cong. for Appl. Mech., Stockholm, Sweden, 1930. 28. Malkus, W., "Discrete Transitions in Turbulent Convection," Proc. Roy. Soc. London, A 225, (1954), 185-195. 29. Malkus, W., "The Heat Transport and Thermal Spectrum of Thermal Convection," Proc. Roy. Soc. London, A 225, (1954), 196-212.

-10730. Malkus, W. And Veronis, G.,. "Finite Amplitude Cellular Convection," J. of Fluid Mech., 4, (1958), 225-260. 31. Martini, W., Natural Convection Inside a Horizontal Cylinder, Ph.D. Thesis, Univ. of Michigan, Ann Arbor, Mich., 1959. 32. Milne-Thompson, L.M., Theoretical Hydrodynamics,(4th ed.), MacMillan, London, 1960. 33. Morgan, G.E., and Warner, W., "On Heat Transfer in Laminar Boundary Layers at High Prandtl Number, " J. of Aeronautical Sciences, 23, (1956), 937. 34. Mull, W., and Reiher, H., "Der Warmeschutz von Luftschichten," Gesund.-Ing. Beihefte, 28, 1930. (As reported in Reference 20.) 35. Nakagawa, Y., "Heat Transport in Cellular Convection," Physics of Fluids, 3, (1960), 82. 36. O'Brien, G., Hyman, M., and Kaplan, S., "A Study of the Numerical Solutions of Partial Differential Equations," J. of Math. Phys., 29, (1951), 223-251. 37. Ostrach, S., "Convection Phenomena in Fluids Heated from Below," Trans. ASME, 9, (1957), 229. 38. Palm, E., "On the Tendency Towards Hexagonal Cells in Steady Convection," Journal of Fluid Mechanics, 8, (1960), 183-92. 39. Peaceman, D., and Rachford, J., "The Numerical Solution of a Parabolic and Elliptic Differential Equation," J. of the Soc. Ind. and Appl. Math., 3, No. 1, (1955), 28-41. 40. Pellew, A., and Southwell, R., "On Maintained Convective Motion in a Fluid Heated from Below," Proc. Roy. Soc. London, A 176, (1940), 312-343. 41. Poots, G., "Heat Transfer by Laminar Free Convection in Enclosed Plane Gas layers," Quart. Mech. and Appl. Math., 11, (1958), 257-273. 42. Rayleigh, J.W.S., "On Convection Currents in a Horizontal Layer of Fluid when the Higher Temperature is on the Under Side," Phil. Mag., 32, (1916), 529-546. 43. Reid, W.H., and Harris, D.L., "Streamlines in Benard Convection Cells," Phys. of Fluids, 1, (1958), 102.

-10844. Richards, C.G., Flow Due to a Rotating Disk with Center Sink, Ph.D. Thesis, University of Michigan, Ann Arbor, Mich., 1964. 45. Richtmeyer, R.D., Difference Methods for Initial Value Problems, Interscience Publishers, Inc., New York, 1957. 46. Schlichting, H., Boundary Layer Theory, McGraw-Hill Book Co., New York, 1962. 47. Schmidt, E., Heat Transfer by Natural Convection, presented as a Lecture at the International Heat Conference at the University of Colorado, Boulder, Colorado, 1961. 48. Schmidt, R.J., and Milverton, S.W., "On the Stability of a Fluid when Heated from Below," Proc. Roy. Soc. London, A 152, (1935), 586-594. 49. Schmidt, R.J., and Saunders, O.A, "On the Motion of a Fluid Heated from Below," Proc. Roy. Soc. London, A 165, (1938), 216-228. 50. Schmidt, E., and Silveston, P., "Natural Convection in Horizontal Liquid Layers," C.E.P. Symp. Series No. 29, 55, (1959), 163. 51. Segel, L.A., and Stuart, J.T., "On the Question of the Preferred Mode in Cellular Thermal Convection," J. of Fluid Mech., 13, (1962), 289-306. 52. Sorokin, M.P., "Thermal Instability of a Liquid in a Horizontal Cylinder," Inzh. Fiz. Zh., 4, No. 8, (1961), 107-110. 53. Stommel, H., "A Summary of the Theory of Convection Cells," Annals of the N.Y. Acc. of Sci., 48, (1947), 715-726. 54. Sutton, O.G., "On the Stability of a Fluid Heated from Below," Proc. Roy. Soc. London, A 204, (1951), 297-309. 55. Thomson, J., "On a Changing Tesselated Structure in Certain Liquids," Proc. Glascow Phil. Soc., 13, (1882), 464-468. 56. Wilkes, J.O., The Finite Difference Computation of Natural Convection in an Enclosed Rectangle, Ph.D. Thesis, Univ. of Michigan, Ann Arbor, Mich., 1963.

APPENDIX A SOLUTION OF A SYSTEM OF LINEAR EQUATIONS HAVING A TRI-DIAGONAL COEFFICIENT MATRIX Application of the implicit formulation to one-dimensional problems, or the ADI formulation to two- and three-dimensional problems produces a system of linear equations which have a tri-diagonal coefficient matrix as follows: b v+CoV'l = do alVo+blvl+clv2 = dl a2 v l+b2v2+c2v = d anv +b v + v = d n- lVn-2 +bn-l v n-l n-lVn n-l anvnl + bnvn = dn n n where the ab, c, and d's are known constants, and the v's are the -109 - may~~~~~~~~~~~ be qucl eemndfo* h euso omle *~ ~~~-0

-110v = i vi = c- civi+l i = n-l, n-2,...., 0 where the y and P are determined from do P = bo o =o i= b. aci- di-aiyil i = 1i,2...., n m i-1

APPENDIX B DESCRIPTION OF THE PROGRAMS AND VARIABLES USED IN THE NUMERICAL STUDY 1. Programs The program used for the numerical integration of the vorticity and energy equations consisted of one main program and ten subroutines. The function of each of these was as follows: MAIN - reads in the initial conditions, sets the necessary constants, and calls on the various subroutines to advance t, T, *, U, V, and Nu one time step. After each time step the program decides whether or not to print t, T, 4, U, V, and Nu, and then either advances another time step or punches the T and t fields and reads another set of initial conditions. PRTOUT - does the actual printing and labeling of the fields. INSUL - advances the temperature field through both halves of the time step. VORTIC - advances the vorticity field through both half time steps. STRMFN - calculate the new field of stream function from the new vorticity field. ZWALL - calculates the new wall vorticity from the new stream functions. VEL - evaluates U and V velocity fields from the stream functions. - 111

-112THOM and TOM - solve the tridiagonal matrices generated by VORTIC, INSUL, and STRMFN. DELY - evaluates 6, and 6 of T and ~ as used by INSUL, y y VORTIC, and STRMFN. GRADY - evaluates Nu from the temperature field. In addition to the main program and various subroutines used for the numerical integration, various other programs were used to determine the optimum time step for the stream function iteration; plot the streamlines and isotherms; and prepare initial conditions for nonsquare cavities. These programs are discussed below: CHCK - used to determine the optimum sequence of time step sizes for the stream function iteration. PLOT - locates and plots the streamlines and isotherms from the T and * fields punched at the end of each calculation. SPLIT - used to convert an l1xll field into an 11x21, 21x11, or 21x21 field to simplify punching the initial conditions for L/H = 2.0 and 0.5 cases. SETUP - used to simplify punching the initial conditions for the L/H = 3.0 case. Several stream functions at specified grid points were read as data. Zeros were then placed in the remainder of the grid points. A linear temperature gradient is assumed for the temperature field, and the T and 4 were punched on cards to provide the desired initial conditions.

-1132. Variables A description of the variables used in the numerical integration of the vorticity and energy equations is given below: Name Physical Description Interpretation CARD whenever CARD=lB the final T and fields are punched on cards. DT AT time step size for ~ and T fields DX AX X increment between vertical grid lines. DY AY Y increment between horizontal grid lines. EPS1 used in testing the convergence of the stream function iteration FPCH f and T fields are punched after every FPCH time steps FREQ U, V, T, ~, x, and Nu are printed after every FREQ time steps. GR Gr Grashoff number. HEIGHT Height Height of cell —usually 1.0. INTPRT whenever INTPRT=1B the intermediate values TSTAR, ZSTAR, and PHI are printed. ITER number of time steps already performed ITMAX maximum number of iterations allowed in determining new stream functions. When ITMAX is not sufficient, program assumes error has occurred and calculations are ceased.

-114Name Physical Description Interpretation L At/(Ax) dimensionless times step size for stream function iteration. LENGTH Length length of cell. NDATA when NDATA=1B, new T and ~ fields are read into the computer. Otherwise the last set of T and V are used as initial conditions again. NX NX number of X grid spaces. NY NY number of Y grid spaces. PHI intermediate values of ~ for the stream function iteration. PR Pr Prandtl number. PRT whenever ITER is less than PRT,, T, U, V, and Nu are printed. PSI stream functions T T temperature TAUMAX T maximum values of time the max calculations are to be carried. TAUZER T value of time before the first time step. zero TSTAR T* intermediate values of T between time steps. U U X velocity V V Y velocity Z vorticity ZSTAR * intermediate value of ~ between time steps.

- 1153. Program Listings The main program and various subroutines associated with the numerical integration of the vorticity and energy equations are presented in the following pages.

$ COMPILE MAD, PRINT OBJECT* PUNCH OBJECT MAINOO0 R PROGRAM TO NUMERICALLY SOLVE. THE EQUATIONS OF MOTION AND R ENERGY AS APPLYED TO A TWO DIMENSIONAL CAVITY WALL HEATED BR FROM BELOW AND COOLED FROM ABOVE R BOOLEAN CARD, LOOK. NDATA, INTPRT INTEGER FREQ. ITER, ITMAX, SWITCH, NX. NY * It J 1. PRT, FPCH, W, WI, WIJ 2, ITERMX, SWITCH. ITE, NN * IT PROGRAM COMMON WU.VPSI. Z. T * NX, NY, ITER9 DXDY, LX9 1LY. D2Y, L2X9 L2Y, L2X2, L2Y2, R, R2, D2X, DSQX9 L2XPR9 L2YPR 29L2X2PR, L2Y2PR9 LY2GR, LX2PR, LY2PR, LY4,DIMLOOK. DX29 DY2 3, A, C. D, AJo CJ. DJ. LY2 i LX2, EPS1i L 4, INTPRT 9 DY6, DX69 DY129 DX12 DIMENSION U(500~DIM)> V(500,DIM)t PSI(5009 DIM), Z(5009DIM) 1, T(500,DIM)*A(50), C(50) D(50), AJ(50)9 CJ(50)t DJ(50) 1, W(30) t TD(500.DIM). NPSI(500. DIM) VECTOR VALUES DIM = 29 0 0 EXECUTE FTRAP. ITMAX = 0 L=O R.**.*READ AND PRINT INPUT VARIABLES START PRINT COMMENT $1$ READ AND PRINT DATA HEIGHT, LENGTH, NY, NX, GR, PR, DT.ITMAX 1 ITERMX, FREQ 9 TAUZER, TAUMAX9 FPCH READ AND PRINT DATA9 EPS1i L, PRT9 NDATA, 1 INTPRT 9 CARD DIM(1) = NY + 2 DIM(2) = NY + 1 THROUGH SETW. FOR I=0,1 IeG. NX SETW W(I) = ( NY +1 )*I EXECUTE ZERO. ( U(O)U**U(500)t V(O).**V(500)9 PSI(O)0..PSI( 1500), Z(0)...Z(500') T(O)*..T(500). ITER, TAU ) WHENEVER *NOT. NDATA9 TRANSFER TO SET READ FORMAT DATA, (J=0,1,J.G.NY,(I=0,1,I.G.NXTD(W(I)+J))) READ FORMAT DATA9 (J=O,1,J.G.NY ( I=01,I.*G.NX NPSI[W(I) 1 +J))) VECTOR VALUES DATA =$5E16.8 *$ SET THROUGH SET2, FOR I = 01i I.G.NX WI = W( I ) THROUGH SET29 FOR J = 09,1 J.G. NY WIJ = WI + J T(WIJ) = TD(WIJ) SET2 PSI(WIJ) = NPSI(WIJ) PRINT COMMENT $0 AT TAU ZERO THE FOLLOWING CONDITIONS 1EXIST $ RA = PR*GR PRINT RESULTS TAUZER. RA R,...SET CONSTANTS TO BE USED LATER

-117NN = NY + 1 DX = LENGTH/ NX DY = HEIGHT / NY NU = GRADY. ( T, NX, DY ) DX2 = 0.5/ DX DY2 = 0.5/ DY DX6 = 6.0*DX DX12 = 12.0*DX DY12 = 12.0*DY DY6 = 6.0*DY LX = DT/DX LY = DT/DY LX2 = LX/ 2.0 LY2 = LY/2. L2X = DT/DX/DX L2Y = DT/DY/DY L2X2 = L2X/ 2. L2Y2 = L2Y/ 2. R = DX/ DY R2 = R*R D2X = 1.O/DX/DX DSQX = 2.0*( 1.0 + R2) D2Y = 1.0/DY/DY L2XPR = L2X/PR L2YPR = L2Y/PR L2X2PR = L2X2/ PR L2Y2PR = L2Y2/PR LY2GR = LY2*GR LX2PR = LX2/PR LY2PR = LY2/PR LY4 = LY/4. R.....CALCULATE INITIAL VORTICITIES AND VELOCITIES EXECUTE VELOC. THROUGH ONE, FOR J=1,1,J.G.NY-1 THROUGH ONE, FOR I = 1,1,IG. NX-1 WIJ = W(I) + J ONE Z(WIJ) = - D2X*((PSI (WIJ+NN) + PSI( WIJ-NN))+ R2*(PSI (WIJ+) 1 + PSI (WIJ - 1))- DSQX* PSI( WIJ)) EXECUTE ZWALL. EXECUTE PRTOUT. ( TZ, PSI, U, V, W, NX, NY R...**ADVANCE VORTICITIES AND TEMPERATURES TAUZER = TAUZER + DT THROUGH OUT1, FOR TAU = TAUZER, DT, TAU.G. TAUMAX LOOK = OB ITE = ITER +1 WHENEVERITER.L.PRT.OR.((ITE /FREQ)*FREQ-ITE ).E.O, LOOK = 18 EXECUTE INSUL. EXECUTE VORTIC. R.....CALCULATE NEW STREAM FUNCTIONS EXECUTE STRMFN.( ITERMX. IT ) R.....CALCULATE NEW WALL VORTICITIES

-118EXECUTE ZWALL. R*..**CALCULATE NEW VELOCITIES EXECUTE VELOC. ITER = ITER + 1 R PUNCH RESULTS WHEN DESIRED WHENEVER (( ITER/FPCH)*FPCH - ITER) *E.O PUNCH FORMAT DATA,(J=O,1,J*G.NY.(I=01.1,IG.NXT(I.J))) PUNCH FORMAT DATA,(J=O,1,J.GeNY,(I=O,1,I*.GNX,PSI(IJ))) END OF CONDITIONAL R.....PRINT RESULTS WHEN DESIRED WHENEVER.NOT. LOOK. TRANSFER TO OUT1 PRINT COMMENT $4$ PRINT RESULTS TAU, GR, PR * IT EXECUTE PRTOUT. ( T,Z, PSI, U. V, W. NX, NY ) NU = GRADY. ( T.NX, DY ) PRINT RESULTS NU OUT1 CONTINUE WHENEVER. NOT. CARD TRANSFER TO END PUNCH FORMAT DATA,(J=O,1,J*G.NY,(I=O01,I.G.NXT(IJ))) PUNCH FORMAT DATA,(J=O,1,J.G.NY*(I=,01,I.G*NXPSI(IJ))) END TRANSFER TO START END OF PROGRAM $ COMPILE MAD, PRINT OBJECT, PUNCH OBJECT PRTOUTOOO EXTERNALFUNCTION (T, Z,PSI. U, V, W, NX, NY ) ENTRY TO PRTOUT. INTEGER N INTEGER NX, NY, W, I, J WHENEVER NX.G. 10, TRANSFER TO BIG PRINT COMMENT $0 VORTICITY FIELD IS GIVEN BY $ PRINT FORMAT OUT, (J=O,1,J*GeNY,(I=O,1,I.G.NX, Z(W(I)+J))) PRINT COMMENT $0 NEW FIELD OF TEMPERATURE IS GIVEN BY $ PRINT FORMAT OUT, (J=O,1J.G.NY,(I=O,1,I.G.NX, T(W(I)+J))) PRINT COMMENT $0 NEW FIELD OF STREAM FUNCTION IS GIVEN BY $ PRINT FORMAT OUT. (J=O,1,J.G.NY,(I=O,1,I.G.NXPSI(W(I)+J))) PRINT COMMENT $0 THE U COMPONENT IS $ PRINT FORMAT OUT. (J=O,1,J.G.NY.(I=O01I.G.NX, U(W(I)+J))) PRINT COMMENT $ THE V COMPONENT IS $ PRINT FORMAT OUT, (J=01J.GeNNY,(I=,01,I*G.NX, V(W(I)+J))) VECTOR VALUES OUT = $ 1HO/( 1H, llE114)*$ VECTOR VALUES OUT1= $ 1HO/( 1H, 10E11.4)*$ FUNCTION RETURN BIG PRINT COMMENT $0 VORTICITY FIELD IS GIVEN BY $ PRINT FORMAT OUT, (J=Ol0J.G.NY,(I=O01,I.G.10, Z(W(I)+J))) THROUGH P2, FOR N=11, 10, N.G. NX P1 PRINT FORMAT OUT1,(J=O,1,J.G.NY,(I=N.1iI.G.N+9,Z(W(I)+J))) PRINT COMMENT $0 NEW FIELD OF TEMPERATURE IS GIVEN BY $ PRINT FORMAT OUT. (J=O,1,J.G.NY,(I=O01lI.G.10' T(W(I)+J))) THROUGH P2, FOR N=11. 10. N.G. NX

-119 - P2 PRINT FORMAT OUT1,(J=O,1,J.G.NY,(I=N,1,IoG.N+9.T(W(I)+J))) PRINT COMMENT $0 NEW FIELD OF STREAM FUNCTION IS GIVEN BY $ PRINT FORMAT OUT, (J=O,1,J.G.NY,(I=0,1,I.G*10PSI(W(I)+J))) THROUGH P3, FOR N=11 10, N.G. NX P3 PRINT FORMAT OUT1,(J=0,1J.G.NY,(I=NtiI*G.N+9, PSI(I,J))) FUNCTION RETURN END OF FUNCTION $ COMPILE MAD, PRINT OBJECT, PUNCH OBJECT INSULOO0 EXTERNAL FUNCTION ENTRY TO INSUL. BOOLEAN LOOK, INTPRT INTEGER I,J,NX, NY * ITER, WI WI WJ PROGRAM COMMON W,U,V,PSI, Z, T, NX, NY, ITER, DXDY, LX, iLY, D2Y, L2X, L2Y, L2X2, L2Y29 R, R2, D2X, DSQX, L2XPR, L2YPR 2,L2X2PR, L2Y2PR, LY2GR, LX2PR, LY2PR, LY4,DIMLOOK, DX2, DY2 3, A, C, D, AJ, CJ, DJ, LY2, LX2, EPS1, L 4, INTPRT, DY6, DX6, DY12, DX12 DIMENSION U(500,DIM), V(500,DIM), PSI(500, DIM), Z(500,DIM) 1, T(500,DIM),A(50), C(50), D(50), AJ(50). CJ(50), DJ(50) 2,TSTAR(5009 DIM), B(50), BJ(50), DIM(2), W(30) WHENEVER ITER.NE. 0O TRANSFER TO GO CX = -2.0* ( L2XPR + 1.0 ) CY = - ( 1.0 + L2YPR ) EXECUTE ZERO. ( TSTAR(0)...TSTAR(500)) THROUGH NEXT, FOR I = 0, 1, I.G. NX TSTAR( I,0) = T(IO) NEXT TSTAR(INY) = T(I.NY) R.....SET COEFF TO ADVANCE TEMPS, X DIRECTION A(O) = 0. B(0) = -2.*( L2XPR + 1 ) C ( 0 ) = 2. * L2XPR A(NX) = C(0) B(NX) = B(0) C(NX) = 0.0 THROUGH SET, FOR I - 1 1, I.G. NX-1 SET B(I) = CX GO THROUGH HALF, FOR J = 11, J.G.NY-1 THROUGH SET1, FOR I = 1,,1I.G. NX-1 WIJ = W(I) +J A(I) = L2XPR + LX2*U(WIJ) C(I) = L2XPR - LX2*U(WIJ) SET1 D(I) = V(WIJ)*LY2*DELY.(TI,J) -L2YPR*DELSQY.( T,I9J) - 2* 1T(WIJ) D(0) = -L2YPR* DELSQY. ( T0O,J) - 2.*T( J) D(NX) = -L2YPR*DELSQY.( TNXJ) - 2.0*T(W(NX)+J) R.... CALCULATE TSTAR BY THOMAS ELIMINATION

-120HALF EXECUTE THOM. ( TSTAR, 1B, A,B, C, D, NX, O,J, W WHENEVER.NOT. LOOK, TRANSFER TO NO WHENEVER.NOT. INTPRTt TRANSFER TO NO PRINT COMMENT $0 INTERMEDIATE TEMPS GIVEN BY $ PRINT FORMAT OUT, (J=O,1,J.G.NY,(I=O,1,I.G.NX,TSTAR(IJ))) VECTOR VALUES OUT = $ 1HO/( 1H, 11El1.'4)*$ NO CONTINUE WHENEVER ITER.NE. O0 TRANSFER TO G02 R,....NOW BEGIN SECOND STEP IN TEMP. THROUGH SET3, FOR J=1l1, J.G. NY SET3 BJ(J) = CY G02 THROUGH DONE, FOR I = 0, 1, I.G.NX WI = W(I) Ro....SET COEFF MATRICES SECOND STEP THROUGH SET2, FOR J = 1, 1, J.G. NY-1 WIJ = WI + J AJ(J) =L2Y2PR + LY4*V(WIJ) CJ(J)= L2Y2PR - LY4*V(WIJ) SET2 DJ(J)= T(WIJ) - 2.0*TSTAR(WIJ) +L2Y2PR*DELSQY.(T,IJ) - V(WIJ 1)*LY4 *DELY.(TI,J) DJ(1) = DJ(1) - AJ(1)*T(WI ) DJ(NY-1) = DJ(NY-1) - CJ(NY-1)*T(WI + NY ) DONE EXECUTE THOM.( T, OB, AJ, BJ, CJ, DJ, NY-1, 1, I W ) WHENEVER.NOT. LOOK, TRANSFER TO N02 WHENEVER.NOT. INTPRT9 TRANSFER TO N02 PRINT COMMENT $0 FINAL TEMPS ARE GINEN BY $ PRINT FORMAT OUT9 (J=0,1,J.G.NY,(I=01, I.G.NX,T ( IJ))) N02 FUNCTION RETURN END OF FUNCTION $ COMPILE MAD, PRINT OBJECT, PUNCH OBJECT VORTIOOO EXTERNAL FUNCTION ENTRY TO VORTIC. BOOLEAN LOOK, INTPRT INTEGER I, J, NX, NY, ITER, W, WI, WIJNN PROGRAM COMMON W,U,VtPSI, Z, T, NX, NY, ITER, DX,DY9 LX, 1LY, D2Y, L2X, L2Y, L2X2, L2Y2, R, R2, D2X, DSQX, L2XPR, L2YPR 2,L2X2PR, L2Y2PR, LY2GR, LX2PR, LY2PR, LY4,DIMtLOOK, DX2, DY2 3, A, C, D, AJ, CJ, DJ, LY2, LX2, EPS1, L 4, INTPRT, DY6, DX6, DY12, DX12 DIMENSION U(500,DIM), V(500,DIM), PSI(500, DIM), Z(500DIM) 1, T(500,DIM),A(50), C(50), D(50), AJ(50), CJ(50), DJ(50) 2, DIM(2), ZSTAR(500, DIM), B(50), BJ(50), W(30) WHENEVER ITER.NE. 0, TRANSFER TO GO

-121BX = - 2.0* ( 1. +L2X) BY = - ( 1. +L2Y ) R....eSET COEFF TO ADVANCE VORTICITIES ONE STEP THROUGH SET, FOR I=1,1, I.G. NX-1 SET B(I) = BX GO THROUGH HALF, FOR J=1ii J.G.NY-1 THROUGH SET1, FOR I = 1,t, I.G.NX-1 WIJ W(I) +J A(t) =L2X + U(WIJ)*LX2 C(M) =L2X - U(WIJ)*LX2 SET1 D(I) = V(WIJ)*LY2*DELY.(ZI,J) -L2Y*DELSQY.(ZI,J) - 2.0*Z( 1WIJ) + LY2GR*DELX.( TIJ) D(l) = D(1) - A(1)*Z( J) D(NX-1)= D(NX-1) - C(NX-1)*Z(W(NX) + J) R......ADVANCE FIRST VORTICITY STEP HALF EXECUTE THOM. ( ZSTARt 1B, A.BtC.D~NX-1, 1 J. W ) WHENEVER.NOT. LOOK9 TRANSFER TO NO WHENEVER.NOT. INTPRTt TRANSFER TO NO PRINT COMMENT $0 INTERMEDIATE VORTICITIES ARE GIVEN BY $ PRINT FORMAT OUT. (J=O,1,J.G.NY,(I=01OI.G.NXZSTAR(I.J))) VECTOR VALUES OUT = $ 1HO/( 1H 9 1lE11.4)*$ R.....SET COEFFS FOR SECOND STEP NO CONTINUE WHENEVER ITER.NE. O, TRANSFER TO G02 THROUGH SET39 FOR J = 1919 J.G. NY-1 SET3 BJ(J) = BY G02 THROUGH DONE, FOR I = 1 1 I.G. NX-1 WI = W( I ) THROUGH SET2, FOR J=ll J.G. NY-1 WIJ = WI + J AJ(J)=L2Y2 + V(WIJ)*LY4 CJ(J) = L2Y2 - V(WIJ) * IY4 SET2 DJ(J)= Z(WIJ) - 2.0*ZSTAR(WIJ) +L2Y2*DELSQY.(Z,I,J) - V(WIJ) 1*LY4 *DELY.( ZI,J) DJ(1) = DJ(1) -AJ(1)*Z(WI) DJ(NY-1) = DJ(NY-1) - CJ(NY-1)*Z(WI + NY) R......ADVANCE SECOND STEP IN VORTICITY DONE EXECUTE THOM. ( Z. OB, AJ, BJ, CJ, DJ, NY-1i 1* I 9 W ) WHENEVER.NOT. LOOK9 TRANSFER TO N02 WHENEVER.NOT. INTPRT, TRANSFER TO N02 PRINT COMMENT $0 FINAL VORTICITIES ARE $ PRINT FORMAT OUT. (J=0,1,J.G.NY,(I=O,1,I.G.NXZ (I,J))) N02 FUNCTION RETURN END OF FUNCTION $ COMPILE MAD, PRINT OBJECT. PUNCH OBJECT STRMFOO0 EXTERNAL FUNCTION (ITERMX, IT )

-122ENTRY TO STRMFN. BOOLEAN LOOK t BOOL, INTPRT INTEGER IT, ITER, I, J, M, NX, NY, W, WI, WIJX, Y. ITERMX 1,PTMX, PTFQ PROGRAM COMMON W,U,VPSI9 Z, T, NX, NY, ITER, DXDY, LX, 1QY, D2Y, L2X, L2Y, L2X2, L2Y2, R, R2, D2X, DSQX, L2XPR, L2YPR 2,L2X2PR, L2Y2PR, LY2GR, LX2PR, LY2PR, LY4t DIM, LOOK. DX2,DY2 3, A, C, D, AJ, CJ, DJ, LY2, LX29 EPS1, L 4, INTPRT, DY6, DX6, DY12, DX12 DIMENSION U(500,DIM), V(500,DIM), PSI(500, DIM), Z(500,DIM) 1, T(500,DIM),A(50)* C(50) D(50), AJ(50), CJ(50), DJ(50) 2,DIM(2), PHI(500,DIM),W(30), D(30), DD(30), 3 G(30), GG(30), BETA(50), BETB(50) BOOL = OB WHENEVER ITER *G. 0, BOOL = 1B R L = DT/ (DX)2 PTMX = 2 PTFQ = 10 IT = 0 EXECUTE ZERO. ( PHI (0)...PHI(500)) WHENEVER BOOL, TRANSFER TO OUTER LY = L * R2 DT = L *DX * DX A = L A2 = A*A BETA(1) = -2.0*( 1.0 + L ) THROUGH SETB, FOR I=2,1,I.G. NX-1 SETB BETA( I ) = BETA(i) - A2 / BETA ( I-1) OUTER CONTINUE THROUGH HALF, FOR J=1,1,J.G.NY-1 THROUGH SETD, FOR I=1,1,I.G.NX-1 WIJ = W( I) +J SETD D(I) =-2.0*PSI(WIJ) - LY *DELSQY.(PSI,IJ) 1 - Z(WIJ) * DT G(1) = D(1) / BETA( 1) THROUGH SETG, FOR I:2,1, I.G.NX-1 SETG G(I) = ( D(I) - A *G(I-1))/ BETA( I) HALF EXECUTE TOM.( PHI, AG,BETA, 1BNX-1lJ, W ) R.*...PERFORM SECOND STEP WHENEVER BOOL, TRANSFER TO OUTB AA = L*R2/ 2.0 AA2 = AA*AA BETB(1) - -(R2*L + 1.0) THROUGH SETBB, FOR J= 2, 1, J.G. NY-1 SETBB BETB( J) = BETB( 1) - AA2 / BETB( J-1) OUTB THROUGH DONE, FOR I = 1,1, I.G.NX-1 WI = W(I) THROUGH SETDD, FOR J = 1~1, J.G. NY-1 WIJ = WI + J SETDD DD(J) = PSI(WIJ) + AA *DELSQY.(PSI,IJ)- 2.*PHI(WIJ)

-123GG(1) = DD(1) / BETB(1) THROUGH SETGG, FOR J - 2,1. J.G.NY-1 SETGG GG(J) =( DD(J) - AA *GG(J-1)) /BETB(J) DONE EXECUTE TOM. ( PSI, AA, GG. BETB OBNY-1,I9 W ) IT = IT + 1 WHENEVER *NOT. INTPRT, TRANSFER TO RESET WHENEVER *NOT. LOOK, TRANSFER TO RESET WHENEVER IT *G*PTMX *AND. ((IT/PTFQ)*PTFQ - IT) *NE.O0 1 TRANSFER TO RESET PRINT RESULTS IT PRINT COMMENT $ INTERMEDIATE STREAM FUNCTIONS ARE GIVEN BY $ PRINT FORMAT OUT, (J=0,1,J.GeNY,(I=0,1,I.G*NX, PHI(IJ))) PRINT COMMENT SONEW STREAM FUNCTIONS ARE GIVEN BY $ PRINT FORMAT OUT, (J=0,1J.G.NY,(I=0,1,I.G.NX, PSI(IJ))) VECTOR VALUES OUT = $ 1HO/( 1H, 11E114)*$ RESET CONTINUE BOOL = 1B INTERNAL FUNCTION VOR.(X,Y) =-(PSI(X+1,Y) + PSI(X-1,Y) +R2*( 1PSI(X,Y+1) + PSI( X, Y-1)) - DSQX*PSI(XY))*D2X WHENEVER.ABS. (VOR.(1,2) - Z(1,2)).L. *ABS.( EPS1*Z(1,2)), 1 FUNCTION RETURN WHENEVER IT *L. ITERMX, TRANSFER TO OUTER EXECUTE ERROR. FUNCTION RETURN END OF FUNCTION $ COMPILE MAD, PRINT OBJECT, PUNCH OBJECT VELOCOO0 EXTERNAL FUNCTION ENTRY TO VELOC. INTEGER NX. NY, I, J, ITER, W. WI, WIJ, NN PROGRAM COMMON WUV,PSI, Z, T, NX, NY, ITER, DX.DY, LX, 1LY, D2Y, L2X, L2Y, L2X2, L2Y2, R, R2, D2X, DSQX, L2XPR, L2YPR 2,L2X2PR, L2Y2PR, LY2GR, LX2PR, LY2PR, LY4,DIM,LOOK, DX2, DY2 3, A, C, D, AJ, CJ, DJ, LY2, LX2, EPS1, L 4, INTPRT, DY6, DX6, DY12, DX12 DIMENSION U(500,DIM), V(500,DIM), PSI(500. DIM), Z(500,DIM) 1, T(500,DIM),A(50). C(50), D(50)9 AJ(50), CJ(50), DJ(50) 2, DIM(2), W(30) R.....CALCULATE U VELOCITIES.... THROUGH ONE, FOR I = 1,1, I.G. NX-1 WI = W(I) U(WI+1) = (-3.0*PSI(WI+1)+ 6.0*PSI(WI+2)- PSI(WI+3))/DY6 U(WI+NY-1)=(3.*PSI(WI+NY-1)-6.*PSI(WI+NY-2)+PSI(WI+NY-3))/DY6 THROUGH ONE, FOR J=2,1, J.G.NY-2 WIJ = WI +J ONE U(WIJ) = ( 8.*(PSI(WIJ+1) - PSI(WI J-l)) + PSI( WI J-2) 1 - PSI( WIJ+2))/ DY12 R..... CALCULATE V VELOCITIES.... THROUGH TWO, FOR J= 1,1, J.G. NY-1 V(W(1)+J)=(3.*PSI(W(1)+J) -6.*PSI(W(2)+J) + PSI(W(3)+J))/DX6

-124V(W(NX-1)+J)= (-PSI(W(NX-3)+J)+6.*PSI(W(NX-2)+J)-3.*PSI 1 (W(NX-1)+J))/DX6 THROUGH TWO, FOR I=2,1, I.G. NX-2 WIJ = W(I) +J TWO V(WIJ) =(8.*(PSI(W( I-1)+J)-PSI(W( I+1+J))-PSI (W(I-2)+J) + 1PSI(W(I+2)+J))/DX12 FUNCTION RETURN END OF FUNCTION $ COMPILE MAD, PRINT OBJECT, PUNCH OBJECT DELY.000 EXTERNAL FUNCTION ( A, It J) PROGRAM COMMON W DIMENSION W(30) INTEGER IJ,W, WIJ ENTRY TO DELY. WIJ = W(I) +J FUNCTION RETURN( A(WIJ+1) - A(WIJ-1)) ENTRY TO DELSQY. WIJ = W(I) +J FUNCTION RETURN ( A(WIJ+1) + A(WIJ-1) - 2.*A(WIJ)) ENTRY TO DELX. FUNCTION RETURN ( A(W(I+1)+J) - A(W(I-1)+J)) END OF FUNCTION $ COMPILE MAD, PRINT OBJECT, PUNCH OBJECT TOM.0000 EXTERNAL FUNCTION(PHIA, G, BETA, BOOL, N. K, W) ENTRY TO TOM. BOOLEAN BOOL INTEGER IJ, N, M,K, W 9 WI, WIJ WHENEVER BOOL PHI(W(N)+K) = G(N) THROUGH NEXT1, FOR I=N-1, -1, I.L.1 WIJ = W(I) + K NEXT1 PHI(WIJ) = G(I ) - A * PHI(W(I+1)+K)/ BETA(I) OTHERWISE WI = W(K) PHI(WI +N) = G(N) THROUGH NEXT2, FOR J = N-l1 -1, J.L. 1 WIJ = WI + J NEXT2 PHI( WIJ ) = G(J) - A*PHI( WIJ + 1) / BETA(J) END OF CONDITIONAL FUNCTION RETURN END OF FUNCTION

-125$ COMPILE MAD, PRINT OBJECT, PUNCH OBJECT ZWALLOO0 EXTERNAL FUNCTION ENTRY TO ZWALL. PROGRAM COMMON WUVPSI, Z, T, NX, NY, ITER, DXDY, LX, 1LY, D2Y, L2X, L2Y, L2X2, L2Y2, R. R2, D2X, DSQX, L2XPR, L2YPR 2,L2X2PR, L2Y2PR, LY2GR, LX2PR, LY2PR, LY4,DIMLOOK, DX2, DY2 3, A, C, D, AJ, CJ, DJ, LY2 DIMENSION U(500,DIM), V(500,DIM), PSI(5009 DIM), Z(500,DIM) 1, T(500,DIM),A(50), C(50), D(50), AJ(50), CJ(50), DJ(50) 2, W(30), DIM(2) INTEGER I,JNX,NY, ITER, W, WI, WIJ, NN NN = NY+1 THROUGH SETX, FOR I=1,1,I *G. NX- 1 WI = W(I) Z(WI) - -2*0*PSI( WI+1)*D2Y WIJ = WI + NY SETX Z(WIJ) = -2.0* PSI(WIJ -1) * D2Y THROUGH SETY, FOR J=11,J.G.NY-1 Z(J) = -2.0* PSI( NN+J) * D2X WIJ = W(NX) +J SETY Z( WIJ ) = -2.0*PSI( WIJ - NN) * D2X FUNCTION RETURN END OF FUNCTION $ COMPILE MAD, PRINT OBJECT, PUNCH OBJECT THOM.000 EXTERNAL FUNCTION ( P, BOOL, ABCDt N,M, L, W ENTRY TO THOM. BOOLEAN BOOL INTEGER I, N, M, J,L,WWI, WIJ DIMENSION BETA(50), GAM(50), V(50) BETA(M) = B(M) GAM(M) = D(M) / B(M) THROUGH SET1, FOR I = M+1,1, I *G. N BETA( I ) = B(I) - I) ) * C(I-1)/ BETA(I-1) SET1 GAM(I) = ( D(I) - A(I)*GAM(I-1))/ BETA(I) WHENEVER BOOL P(W(N) + L) = GAM(N) THROUGH SET2, FOR I= N-l, -1, I.L.M SET2 P(W(I)+L) = GAM(I) - C(I)* P ( W(I+1) + L) / BETA(I) OTHERWISE P(W(L) + N) = GAM(N) THROUGH SET3, FOR J= N-1, -1, J.L.M WIJ = W(L) + J SET3 P(WIJ) = GAM(J) - C(J) *P(WIJ+1)/ BETA(J) END OF CONDITIONAL FUNCTION RETURN END OF FUNCTION

APPENDIX C TYPICAL COMPUTER OUTPUT Two sets of typical computer output are presented below. In the first set new initial conditions have just been introduced and the computer has reproduced these initial conditions. In the second set typical T, Z, PSI, U, and V fields and Nu are presented. -126

-127HEIGHT - 1., LENGTH a 1., NX~10, NY - 10, GR * 1.OE 05, PR20.03, OT- 0.002,ITERMX ~ 10, FREQ - 5, TAUZER.C000, TALMAX = C.200, FPCH - 350 * EPS1 z 1.0E-03, L = 6.0, PRT = 0 NDATA = IB, INTPRT0OB, CARD = 1 8 AT TAU ZERO THE FOLLOWING CONDITIONS EXIST TAUZER z.000000, RA = 2SS9.999969 STARTING TEMPERATURES ARE -.1000E 01 -.1OOOE 01 -.1000E 01 -.1000E 01 -.1000E 01 -.1000E 01 -.10O0E 01 -.1000E 01 -.10001 01 -.1O00GE 01 -.10001 01 -.8620E 00 -.8558E 00 -.8248E 00 -.7671E CO -.6871E 00 -.5981E 00 -.5226E 00 -.4847E 00 -.4946E 00 -.53341 00 -.5586E 00 -.7366E 00 -.7274E 00 -.66391 00 -.5486E 00 -.4000E 00 -.2506E OC -.1325E 00 -.68591-01 -.6953E-01 -.12081 00 -.16741 00 -.6297E 00 -.6216E 00 -.5383E 00 -.3911E 00 -.2132E 00 -.4436E-01.9070E-01.1800E 00.2116E 00.1820E 00.1307E 00 -.5388E GO -.5387E 00 -.4595E 00 -.3150E 00 -.1441E 00.1763E-01.1585E 00.2752E 00.3503E 00.3627E 00.3261E 00 -.4481E 00 -.4639E 00 -.4109E 00 -.2947E 00 -.1491E CO -.1525E-07.1491E 00.2947E 00.4109E 00.4639E 00.4481L GO -.32616E 00 -.3627E 00 -.3503E 00 -.2752E 00 -.1585E 00 -.1763E-01.1441E 00.3150E 00.4595E 00.5387E GO.5388L 00 -.1307E 00 -.1820E 00 -.2116E 00 -.1800E 00 -.9070E-01.4436E-C1.2132E 00.3911E 00.5383E 00.6216E 00.6297E 00.1674E 00.1208E 00.6953E-01.6859E-01.1325E 00.2506E 00.4000E 00.5486E 00.6639E 00.7274E 00.7366E 00.5586E 00.5334E 00.4946E 00.4847E 00.5226E 00.5981E OC.6871E 00.7671E 00.8248E 00.8558E 00.8620L 00.1OCOE 01.1000E 01.1000E 01.1000E 01.1000E 01.1000E 01.1000E 01.1300E 01.1000OE 01.1O00E 01.1OO0F 01 STARTING STREAM FUNCTIONS ARE.OO00E 00.OOOE 00 000 00.OOO00E 00.OCCOE 00.OOO0000E OC U000E 00.GOOOE 00.0000 00.0.00E 00.GOOOE 0.OOO00OE 00 -.5054E 00 -.1796E 01 -.4024E 01 -.6493E 01 -.7853E 01 -.7174E 01 -.4670E 01 -.1737E 01 -.1939E-01.0000 00.000OE 00 -.2907E 01 -.7257E 01 -.1291E 02 -.1762E 02 -.1929E 02 -.1721E 02 -.1200E 02 -.5672E 01 -.1026E 01.COOOL CC.0000E 00 -.6150E 01 -.1387E 02 -.2277E 02 -.2946E 02 -.3173E 02 -.2880E 02 -.2132E 02 -.1158E 02 -.3235E 01.COOOE 00.0000E 00 -.8355E 01 -.1878E 02 -.3001E 02 -.3815E 02 -.4100E 02 -.3772E 02 -.2906E 02 -.1728E 02 -.6142E 01.0000t 00.0000E 00 -.8292E 01 -.2006E 02 -.3228E 02 -.4114E 02 -.4436E 02 -.4114E 02 -.3228E 02 -.2006E 02 -.8292E 01.0000E 00.OOCOE 00 -.6142E 01 -.1728E 02 -.2906E 02 -.3772E C2 -.4100E C2 -.3815L 02 -.3001E 02 -.1878E 02 -.8355E 01.0000E 00.0000E 00 -.3235E 01 -.1158E 02 -.2132E 02 -.2880E 02 -.3173E C2 -.2946E 02 -.2277E 02 -.1387E 02 -.6150E 01.COOOE 00.O000E 00 -.1026E 01 -.5672E 01 -.1200E 02 -.1721E 02 -.1929E 02 -.1762E 02 -.1291E 02 -.7257E 01 -.2907E 01.0000E 00.000OE 00 -.1939E-01 -.17371 01 -.46701 01 -.7174E 01 -.78531 01 -.64931 01 -.40241 01 -.17961 01 -.5054E oO -.0000C 0.00006 00.OOOE 00 uOOOE 00.0000 00.OCOOE 00.GCCE OC.00E 00.0OOOE 00 OOOOE 00.0000O 00.COOOE CC NU = 1.792045 THE U COMPONENT IS.OOOOE 00.OOOCE 00.0000E 00.0000E 00.OOOOE 00.0000C E OC 00.00OOE 00 OOOE 00.00OE CO.0OOO 00.0000E 00 -.1629E 02 -.4048E 02 -.7098E 02 -.9460E 02 -.1008E 03 -.8821E 02 -.6111E 02 -.2874E 02 -.4769E 01.COOO 00).0000E 00 -.3067E 02 -.6482E 02 -.9999E 02 -.1213E 03 -.1250E 03 -.1127E 03 -.8677E 02 -.5123E 02 -.1632E 02.0OOOt 00.0000E 00 -.2983E 02 -.61571 02 -.9049E 02 -.1080E 03 -.1143E 03 -.1085E 03 -.9073E 02 -.6209E 02 -.2721E 02.0000E 0.OOOOE 00 -.1158E 02 -.3292E 02 -.4991E 02 -.6107E 02 -.66066 02 -.6481E 02 -.5807E 02 -.4558E 02 -.2761E 02.0000E 00.OOOOE 00.1232E 02.8094E 01.5124E 01.23C5E 01.7749E-05 -.2305E 01 -.5124E 01 -.8094E 01 -.1232E 02.COOO 00.OOOOE 00.2761E 02.4558E 02.5807E 02.6481E 02.6606E 02.6107E 02.4991E 02.3292E 02.1158E 02.COOOF 00.0000E 00.2721E 02.62091 02.9073E 02.1085E 03.1143E 03.1080E 03.9049E 02.6157E 02.2983E 02.COOOE 00.OOOOE 00.1632E 02.5123E 02.8677E 02.1127E 03.1250E 03.1213E 03.9999E 02.6482E 02.3067E 02.OO0000E 00.OOOOE 00.4769E 01.2874E 02.6111E 02.8821E 02.1008E 03.9460E 02.7098E 02.4048E 02.1629E 02.COOO 00.OOOOE 00.OOOE 00.00OOOOE 00 000E 00.0000E 00.0000E O.0000E 00.0000E 00.0000 00.0OOOO 00.0000E 00 THE V COMPONENT IS.OOOE 00.0000E 00.OOOE 00.0000E 00.0000E 00.0 6 OC 0000E 00.OOOOE 00 0000E 00 0000 OE 00.0000E 00.OOOOE 00.8723E 01.1805E 02.2519E 02.2104E 02.4004E 01 -.1725E 02 -.2972E 02 -.2503E 02 -.9490E 01.OOOO 00.OO00OE 00.3653E 02.5198E 02.5541E 02.3429E 02 -.1971E 01 -.3867E 02 -.6168E 02 -.5882E 02 -.3160E 02.000OOOOE GO.OOOOE 00.6995E 02.8628E 02.8266E 02.4729E 02 -.3218E 01 -.5453E 02 -.9103E 02 -.9656E 02 -.6411E 02.OO0000E 00.OO0OE 00.9596E 02.1126E 03.1020E 03.5746E 02 -.2067E 01 -.6219E 02 -.1073E 03 -.1214E 03 -.9361E 02.OOOOE 00.0000E 00.1053E 03.1256E 03.1105E 03.6294E 02.000OO0E 00 -.6294E 02 -.1105E 03 -.1256E 03 -.1053E 03.00001 OC.000OE 00.9361E 02.1214E 03.1073E 03.6219E 02.2067E 01 -.5746E 02 -.1020E 03 -.1126E 03 -.9596E 02.0OOOE 00.0000E 00.6411E 02.9656E 02.9103E 02.5453E 02.3218E 01 -.4729E 02 -.8266E 02 -.8628E 02 -.6995E 02.0000 00.OOOOE 00.3160E 02.5882E 02.6168E 02.3867E 02.1971E 01 -.3429E 02 -.5541E 02 -.5198E 02 -.3653E 02.000OE 00.0000E 00.9490E 01.2503E 02.2972E 02.1725E 02 -.4004E 01 -.2104E 02 -.2519E 02 -.1805E 02 -.8723E 01.0000F 00.OUOOE 00.OOOOE 00.OO0OE 00.0000 00.OOOOE 00.0000E 00.00OE 00.0OOE 00.OOO0E 00.0000E 00.COOO 00 THE VORTICITY FIELD IS GIVEN BY.000OE 00.1011E 03.3591E 03.8048E 03.1299E 04.1571E 04.1435E 04.9340E 03.3474E 03.3877E 01.COOOE CO.1011E 03.2681E 03.4604E 03.5099E 03.3522E 03.1549E 03.1034E 03.2230E 03.3414E 03.2685E 03.3877E 01.5813E 03.2286E 03.2445E 03.4933E 01 -.2313E 03 -.2760E 03 -.1567E 03.8701E 02.3654E 03.4824E 03.2051E 03.1230E 04.5267E 02 -.5062E 02 -.4852E 03 -.7572E 03 -.8383E 03 -.7207E 03 -.3836E 03.1177E 03.5809E 03.6470E 03.1671E 04 -.2012E 02 -.2814E 03 -.8065E 03 -.11COE 04 -.1203E 04 -.1090E 04 -.7646E 03 -.2263E 03.4236E 03.12281 04.1658E 04.1384E 03 -.3600E 03 -.8856E 03 -.1204E 04 -.1315E 04 -.1204E 04 -.8856E 03 -.3600E 03.1384E 03.16586 04.1228E 04.4236E 03 -.2263E 03 -.7646E 03 -.1090E 04 -.1203E 04 -.11OOE 04 -.8065E 03 -.2814E 03 -.2012E 02.1671E 04.6479E 03.5809E 03.1177E 03 -.3836E 03 -.7207E 03 -.8383E 03 -.7572E 03 -.4852E 03 -.5062E 02.5267E 02.1230E 04.2051E 03.4824E 03.3654E 03.8701E 02 -.1567E 03 -.2760E 03 -.2313E 03.4933E 01.2445E 03.2286E 03.5813E 03.3877E 01.2685E 03.3414E 03.2230E 03.1034E 03.1549E 03.3522E 03.5099E 03.4604E 03.26816 03.10111 03.0000E 00.3877E 01.3474E 03.9340E 03.1435E 04.1571E 04.1299E 04.8048E 03.3591E 03.10111 03.00001 00

UNIVERSIFOF MIo GA -HIGAN2 39015 03695 6145 TAU =.010000, GR - l.00OOOOE 05, PR =.030000 VORTICITY FIELC IS GIVEN BY.OOOOE 00.1852E 03.6246E 03.1383E 04.2223E 04.2699E 04.2509E 04.1706E 04.7176E 03.8502E 02 C0OOOE 00.1852E 03.4145E 03.6729E 03.6768E 03.3475E 03 -.3024E 01 -.8540E 02.1500E 03.4333E 03.4118E 03.8502E 02.9453E 03.3192E 03.3694E 03 -.2154E 02 -.3553E 03 -.4126E 03 -.2443E 03.8987E 02.500GE 03.7056E 03.4464E 03.1958E 04.3423E 01 -.2799E 02 -.7217E 03 -.1104E 04 -.1281E 04 -.1112E 04 -.6255E 03.1219E 03.7815E 03.1167E 04.2662E 04 -.1547E 03 -.3392E 03 -.1187E 04 -.1607E 04 -.1809E 04 -.1602E 04 -.1115E 04 -.3182E 03.4762E 03.2063E 04.2686E 04.4519E 02 -.4592E 03 -.1307E 04 -.1794E 04 -.1995E 04 -.1794E 04 -.1307E 04 -.4592E 03.4519E 02.2686E 04.2063E 04.4762E 03 -.3182E 03 -.1115E 04 -.1602E 04 -.1809E 04 -.1607E 04 -.1187E 04 -.3392E 03 -.1547E 03.2662E 04.1167E 04.7815E 03.1219E 03 -.6255E 03 -.1112E 04 -.1281E 04 -.1104E 04 -.7217E 03 -.2799E 02.3423E 01.1958E 04.4464E 03.7056E 03.5000E 03.8987E 02 -.2443E 03 -.4126E 03 -.3553E 03 -.2154E 02.3694E 03.3192E 03.9453E 03.8502E 02.4118E 03.4333E 03.1500E 03 -.8540E 02 -.3024E 01.3475E 03.6768E 03.6729E 03.4145E 03.1852E 03.OOOOE 00.8502E 02.7176E 03.1706E 04.2509E 04.2699E 04.2223E 04.1383E 04.6246E 03.1852E 03.OOOOE 00 NEW FIELD OF TEMPERATURE IS GIVEN BY -. OOOE 01 -.1OOOE 01 -.1000E 01 -.1OOOE 01 -. OOOE 01 -.1OOOE 01 -.OO00E 01 -.OO 00E C1 -.1OOOE 01 -.1OOOE 01 -.100OF Cl -.8652E 00 -.8641E 00 -.8531E 00 -.8285E 00 -.7896E 00 -.7416E OC -.6955E 00 -.6643E 00 -.6559E 00 -.6647E 00 -.6742E 00 -.7355E 00 -.7353E 00 -.7132E 00 -.6635E 00 -.5871E 00 -.4972E 00 -.4134E 00 -.3559E 00 -.3363E 00 -.3486E 00 -.3674E 00 -.6067E 00 -.6097E 00 -.5805E 00 -.5125E 00 -.4092E 00 -.2893E OC -.1772E 00 -.9606E-01 -.6067E-01 -.6808E-01 -.9282E-01 -.4721E 00 -.4809E 00 -.4526E 00 -.3786E 00 -.2636E 00 -.1283E O0.1643E-02.1019E 00.1552E 00.1605E 00.1370E 00 -.3210E 00 -.3376E 00 -.3184E 00 -.2506E 00 -.1379E 00.6519E-08.1379E 00.2506E 00.3184E 00.3376E 00.3210E 00 -.1370E 00 -.1605E 00 -.1552E 00 -.1019E 00 -.1643E-02.1283E 00.2636E 00.3786E 00.4526E 00.4809E 00.4721E UO.9282E-01.6808E-01.6067E-01.9606E-01.1772E 00.2893E OG.4092E 00.5125E 00.5805E 00.6097E 00.6067t 00.3674E 00.3486E 00.3363E 00.3559E 00.4134E CO.4972E OC.5871E 00.6635E 00.7132E 00.7353E 00.7355E 00.6742E 00.6647E 00.6559E 00.6643E 00.6955E OG.7416E OC.7896E 00.8285E 00.8531E 00.8641E 00.8652E 00.100OE 01.1000E 01.1OOOE 01.1OOOE 01.lCOOE 01.1OOE 01.1OO0E 01.1OOOE 01.1OOOE 01.10OOE 01.10OOE 01 IT = 7 NEW FIELD OF STREAM FUNCTION IS GIVEN BY.OOOOE 00.OOOOE 00.OOOOE 00.OOOOE 00.OOOOE 00 COOOE 0.OOOE 00.OOOOE 00.OOOOE 00.OOOOE 00.0OOOL 00.OOOOE 00 -.9260E 00 -.3123E 01 -.6917E 01 -.1112E 02 -.1350E 02 -.1255E 02 -.8531E 01 -.3588E 01 -.4251E 00.00OOt 00.OOOE 00 -.4727E 01 -.1138E 02 -.2020E 02 -.2753E 02 -.3029E 02 -.2730E 02 -.1949E 02 -.9727E 01 -.2232E 01.000OE 00.OOOOE 00 -.9792E 01 -.2117E 02 -.3475E 02 -.4495E 02 -.4871E 02 -.4445E 02 -.3329E 02 -.1860E 02 -.5833E 01.000 00.OOOOE 00 -.1331E 02 -.2847E 02 -.4545E 02 -.5778E 02 -.6232E 02 -.5738E 02 -.4438E 02 -.2676E 02 -.1031E 02 0000 00.0000E 00 -.1343E 02 -.3055E 02 -.4895E 02 -.6234E 02 -.6732E 02 -.6234E 02 -.4895E 02 -.3055E 02 -.1343E 02.OOOOE 00.00JOE 00 -.1031E 02 -.2676E 02 -.4438E 02 -.5738E 02 -.6232E 02 -.5778E 02 -.4545E 02 -.2847E 02 -.1331E 02.COOOE 00.OOOOE 00 -.5833E 01 -.1860E 02 -.3329E 02 -.4445E 02 -.4871E 02 -.4495E 02 -.3475E 02 -.2117E 02 -.9792E 01.OOOOE 00.OOOE 00 -.2232E 01 -.9727E 01 -.1949E 02 -.2730E 02 -.3029E 02 -.2753F 02 -.2020E 02 -.1138E 02 -.4727E 01.0000E 00.OOOOE 00 -.4251E 00 -.3588E 01 -.8531E 01 -.1255E 02 -.1350E 02 -.1112E 02 -.6917E 01 -.3123E 01 -.9260E 00 -.COOOE 00.OOCOE 00.OOOOE 00 OOOOE 0 0.000 OE 00 O0 OOE 00 OOE 00E.000E 00.0 00 OOOOE 00.OOOOE 00.00OOE 00 THE U COMPONENT IS.OOOOE 00.OOOOE 00.OOOOE 00 OOO0000 0 0.00O OE 0 0.000O E 00.OOOE 00.OOOOE 00.OOOOE 00.OOOOE 00.0000O 00.OOOOE 00 -.2632E 02 -.6290E 02 -.1095E 03 -.1448E 03 -.1542E 03 -.1362E 03 -.9675E 02 -.4833E 02 -.1047E 02.OOO E 00.OOOOE 00 -.4801E 02 -.9656E 02 -.1476E 03 -.1774E 03 -.1828E 03 -.1649E 03 -.1281E 03 -.7777E 02 -.2745E 02.0000L 00.OOOOE 00 -.4680E 02 -.9107E 02 -.1334E 03 -.1590E 03 -.1687E 03 -.1590E 03 -.1323E 03 -.9106E 02 -.4305E 02.000E 00.OOOOE 00 -.1960E 02 -.4974E 02 -.7454E 02 -.9106E 02 -.9740E 02 -.9387E 02 -.8273E 02 -.6406E 02 -.4142E 02.COOOE 00.OOOOE 00.1666E 02.9270E 01.5927E 01.2273E 01.5960E-05 -.2273E 01 -.5927E 01 -.9270E 01 -.1666E 02.COOOE 00.OOOOE 00.4142E 02.6406E 02.8273E 02.9387E 02.974GE 02.9106E 02.7454E 02.4974E 02.1960E 02.COOOE 00.OOOOE 00.4305E 02.9106E 02.1323E 03.1590E 03.1687E 03.1590E 03.1334E 03.9107E 02.4680E 02.OOOOE 00.OOOOE 00.2745E 02.7777E 02.1281E 03.1649E 03.1828E 03.1774E 03.1476E 03.9656E 02.4801E 02.000OE 00.000OOE 00.1047E 02.4833E 02.9675E 02.1362E 03.1542E 03.1448E 03.1095E 03.6290E 02.2632E 02.COOOE 00.OOOOE 00 OOOOE 00.OOOE 00.OOOOE 00.OOOOE 00.0000E 0 0 O.OOOE 00.OOOO E 00.OOOOE 00.OOOOE 00E THE V COMPONENT IS.OOOOE 00. 0000E 0. 0 0 OOO 0 OOOE 00.OOOOE 00 OOOOOE 00 OOOE 00 000 00 OO OOOOE 00.0OOOE 00.COOOE 00.OOOOE 00.1507E 02.3068E 02.4281E 02.3600E 02.8190E 01 -.2682E 02 -.4883E 02 -.4358E 02 -.1953E 02.0000E 00.OOOOE 00.5650E 02.8019E 02.8634E 02.5402E 02 -.8845E 00 -.5717E 02 -.9381E 02 -.9230E 02 -.5362E 02.OOOOE 00.OOOOE 00.1048E 03.1289E 03.1261E 03.7366E 02 -.2120E 01 -.8078E 02 -.1366E 03 -.1460E 03 -.1013E 03.0000E 00.OOOOE 00.1424E 03.1661E 03.1546E 03.8834E 02 -.1797E 01 -.9371E 02 -.1608E 03 -.1793E 03 -.1420E 03.OOOOE 00.OOOOE 00.1568E 03.1848E 03.1670E 03.9597E 02.3974E-06 -.9597E 02 -.1670E 03 -.1848E 03 -.1568E 03.000OE 00.OOOOE 00.1420E 03.1793E 03.1608E 03.9371E 02.1797E 01 -.8834E 02 -.1546E 03 -.1661E 03 -.1424E 03.OOOOE 00.OOOOE 00.1013E 03.1460E 03.1366E 03.8078E 02.2120E 01 -.7366E 02 -.1261E 03 -.1289E 03 -.1048E 03.OOOOE 00.OOOOE 00.5362E 02.9230E 02.9381E 02.5717E 02.8845E 00 -.5402E 02 -.8634E 02 -.8019E 02 -.5650E 02.OOOOE 00.ObOOE 00.1953E 02.4358E 02.4883E 02.2682E 02 -.8190E 01 -.3600E 02 -.4281E 02 -.3068E 02 -.1507E 02.OOOOE 00.OOOOE 00.OOOOE 00.OOOE 00 00E 00.OOOOE 00. OOOOE 00.OOOOE 0.OE 00.OOOO E 00 OOOO E 00.00 00OE 00 NU = 1.251920