THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COTLLEGE OF ENGINEERING NATURAL CONVECTION INSIDE A HORIZONTAL CYLINDER William R. Martini This dissertation was submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan, 1956. Committee in Charge: Associate Professor Stuart W. Churchill, Chairman Assistant Professor Thomas C. Adamson Assistant Professor John W. Carr III Professor Donald L. Katz Professor Harold A. Ohlgren Assistant Professor Alexander Weir August, 1956 IP-176

ACKNOWLEDGEMENT We wish to express our appreciation to the author for permission to give this thesis limited distribution under the Industry Program of the College of Engineering.

ACKNOWLEDGI'IENTS Without the encouragement and active assistance of my wife and of my parents, I would not have entered school for an advanced degree. Without the aid of the Standard Oil Company of California Fellowship in the year 1954-55, and the Horace H. Rackham Predoctoral Fellowship in the year 1955-56, I would not have been able to devote full time to study. I would like to express appreciation for the advice and encouragement given by the members of my doctoral committee. I would especially like to thank Dr. Stuart W. Churchill who, as chairman, devoted considerable time and thought to this problem and offered many suggestions that shaped the character of this thesis. I would like to thank Dr. John W. Carr III whose interest and advice in connection with the numerical solutions of the thesis aided greatly in understanding the problemn and in obtaining meaningful answers. I would also like to acknowledge the influence and encouragement of Dr. G. K. Batchelor who suggested the circular geometry studied in this thesis. Finally, I would like to acknowledge the considerable financial support given me by the Industry Program of the University of Michigan in the reproduction of this thesis. ii

TABLE OF CONTENTS Page ACKNOWLEDGEMENT......................... i AUTHOR' S ACKNOWLEDGEMENTS. ii LIST OF TABLES..................... v LIST OF FIGURES..................... vi NOMENCLATURE.......................... x ABSTRACT............ xiii INTRODUCTION.......................,. 1 PART I. LITERATURE REVIEW. 4 Heat Transfer Correlations............. 5 Mathematical Methods......... 10 Analytical Approach........... 11 Numerical Approach. 15 PART II. EXPERIMENTAL INVESTIGATION........... 28 Description of Apparatus... 29 The Geometry................. 29 Design of the Equipment....... 29 Construction of Experiment..32 Temperature Measurement..... 33 Velocity Measurement............. 38 Data............... 43 Velocity Data................ 43 Temperature Data.......... 43 Analysis of Data........... 49 Velocity Correlations............. 49 Temperature Correlations........... 53 Conclusions.................. 69 PART III. NUMERICAL INVESTIGATION............ 72 Formulation of the Differential Equations...... 74 Basic Equations................ 74 Discussion of Dimensionless Groups.. 77 Simplification of Basic Equations...... 79 Conversion to Difference Equations......... 81 The Grid................... 82 Translation to Difference Equations...... 82 The Iteration Process............. 85 Analysis of Stability............... 88 Stability Analysis of Equations (52) to (54) using the von Neuman Method......... 89 iii

TABLE OF CONTENTS (Cont.) Page Stability Analysis of Equations (52) to (54) using the Matrix Method........... 95 Comparison of Analysis of Stability by the Matrix Method and by the von Neuman Method. 102 Computational Results............... 102 Iteration on T and p Together......... 103 Iteration on p Alone.............. 127 Iteration on T Alone.............. 133 Conclusions.................. 155 CONCLUSION........................... 156 REFERENCES.............. 159 APPENDIX A. FLOW DATA..... 162 APPENDIX B. TEMPERATIRE PROFILE AND RESULTS OF COMPUTATION. 166 APPENDIX C. COMPUTATIONAL PROGRAM DIRECTIONS......... 174 APPENDIX D. CALIBRATIONS.......... 180 iv

LIST OF TABLES Table Page I Calibration of Flash Commutator........... 41 II Local Experimental Velocities............ 50 III Geometric Data.................... 55 IV Analysis of Temperature Data for Experiment 4.... 59 V Summary of Averaged Measures of Heat Transfer Rate.. 64 VI Dimensionless Groups for Experiments 3-15....... 68 VII Boundary Conditions for Some of the Problems Attempted..................... 104 VIII Selected Tangential Velocities p at End of Problem 18....................107 IX Comparison of Solution Progress for Small Time Steps. 114 X Maximum Correction Per Iteration for Small Time Steps....................... 115 XI Estimate of Experimental Velocities for Experiment Number 4...................... 131 XII Temperature Profile Data for Experiments 3, 4, and 6. 167 XIII Temperature Profile Data for Experiments 7, 8, 9 and 10..... 168 XIV Temperature Profile Data for Experiments 11, 12, 13, and 14........................ 169 XV Temperature Profile for Experiment 15........ 170 XVI Results of Computation for Experiments 3, 5, 6, and 7....................... 171 XVII Results of Computation for Experiments 8, 9, 10, and 11.......... 172 XVIII Results for Computation for Experiments 12, 13, 14 and 15........................ 173 XIX Drum Locations for the Program............ 176 v

LIST OF FIGURES Figure Page 1 Geometry of Experiment by Jones and Masson...... 7 2 Rate of Heat Transfer Across a Tall Rectangular Air Space, Robinson and Powlitch........... 9 3 Sample of Grid for Heat Equation Calculation Example. 17 4 Cartesian Co-ordinate Space-Time Nomenclature.... 17 5 Bound Circles for One-Dimensional Heat Flow...... 26 6 Geometry of Experiment................ 30 7 Plan of Test Object.................. 31 8 Diagram of Experimental Setup........ 34 9 Thermocouple Positions on Test Object..... 35 10 Gas Thermocouple Positioning Device......... 37 11 Gas Thermocouple Location for Each Point Number... 39 12 Dust Particle Illumination System.......... 40 13 Sample Dust Particle Picture............. 44 14 Experimental Flow Pattern (Experiment 4)....... 45 15 Horizontal Gas Temperature Profiles (Experiment 4).. 46 16 Vertical Gas Temperature Profiles (Experiment 4). 47 17 Three Dimensional Representation of Temperature Surface Observed in Experiment 4.. 48 18 Correlation of Selected Local Velocities....... 51 19 Velocity and Temperature Profiles for Horizontal Radius at Cold Side (Experiments 8, 9, and 10).... 52 20 Correlation of Boundary Layer Circulation....... 54 21 Slope Number Arrangement............... 56 22 Local Nusselt Number as a Function of Azimuth Angle (Experiment 4).................... 60 vi

LIST OF FIGURES (Cont.) Figure Page 23 Smoothing of Local Nusselt Numbers.......... 62 24 Smoothed Local Nusselt Number as a Function of Azimuth Angle...................... 63 25 Rate of Heat Transfer as a Function of Wall Temperature Difference................ 65 26 Dimensionless Heat Transfer Correlation........ 67 27 Comparison of Experimental Heat Transfer Correlations. 70 28 Polar Coordinate System................ 74 29 The Polar Grid.................... 83 30 Grid Point Numbers.................. 84 31 Polar Coordinate Space-Time Grid Nomenclature..... 89 32 Assumed Round Off Error Frequency Distribution.... 91 33 Plan for Error Matrix-Error Vector Multiplication.. 97 34 Matrix Method Stability Diagrams........... 100 35 Transient Solution of Temperature Equation (Problem 16)...................105 36 Transient Solution of Tangential Velocity Equation (Problem 16).................... 106 37 Comparison of Transient Solutions of the Temperature Equation for Different Time Steps (Point 29)..... 108 38 Maximum Correction in T During Iteration (Problem 18). 109 39 Maximum Correction in p During Iteration (Problem 18). 110 40 Transient Solution of Temperature Equation (Point 29, Problem 20)..................... 112 41 Transient Solution of Temperature Equation at- Very Small Time Steps..................113 42 Progress in Solution of Velocity Equation (Problem20).....................116 vii

LIST OF FIGURES (Cont.) Figure Page 43 Maximum Correction in T During Iteration (Problem 20). 118 44 Maximum Correction in p During Iteration (Problem 20). 119 45 Vector Representation of Velocities at End of Problem 30...................... 120 46 Isotherms at End of Problem 30............ 122 47 Maximum Correction During Iteration (Problem 31)... 123 48 Transient Solution of Velocity Equation with T Constant (Problems 41, 42, and 43)......... 128 49 Maximum Correction in p During Iteration (Problems 41 and 42)................. 130 50 Assumed Velocity Distribution (Problems 26, 30, 39, and 48)....................... 132 51 Transient Solution of Heat Equation (Problem 36)... 134 52 Maximum Correction in T During Iteration (Problem 36). 135 53 Isotherms for Pure Conduction (Problem 36)...... 137 54 Transient Solution of Temperature Equation (Problems 26, 39, and 48)................... 138 55 Maximum Correction in T During Iteration (Problem 26). 139 56 Transient Solution of Temperature Equation (Problem 39).................... 141 57 Calculated Isotherms (Problem 39, t = 125)...... 142 58 Maximum Correction in T During Iteration (Problem 39). 143 59 Transient Solution of Temperature Equation (Problem 48 Initial Solution)............. 144 60 Transient Solution of Temperature Equation (Problem 48 Complete Solution)............ 145 61 Maximum Correction in T During Iteration (Problem 48). 147 viii

LIST OF FIGURES (Cont.) Figure Page 62 Calculated Isotherms (Problem 48) t = 401..... 148 63 Three Dimensional Representation of Temperature Surface (Problem 48, T = 401)............. 149 64 Experimental Dimensionless Isotherms (Experiment 4).. 150 65 Comparison of Experimental and Computed Local Dimensionless Temperature Gradients......... 151 66 Experimental Flow Lines and Velocities (Experiments 3, 6, 7, and 8)........... 163 67 Experimental Flow Lines and Velocities (Experiments 9, 10, 11, and 12)...........164 68 Experimental Flow Lines and Velocities (Experiments 13, 14, and 15)........... 165 69 Calibration of Copper-Constantan Thermocouples.... 181 70 Calibration of Chromel P-Alumel Thermocouple..... 182 ix

NOMENCLATURE (A, B, C, D, E, F, G)ijk = coefficients of equation (67). b = error in C = maximum number of internal grid points in a row or column. C = heat capacity, Btu/#OF c = error in q IV= total differential operator D = diameter of cylindrical space, ft. D = exact solution of differential equation (see page E = distance above prism (see Figure 1), ft. Eu = Euler Number = gc Pref /Dg Pref Gr = Grashof Number based on diameter of cylinder D3gp(AT)/v2. GrL = Grashof Number based on width of gas layer = L3g(AT)/v2. GrH = Grashof Number based on height of gas layer = H3g(ZAT)/v2. G. = acceleration due to gravity in ith direction times p, ft (lb)/sec2 ft.3 g = acceleration due to gravity. = acceleration due to gravity in x direction ft/sec. gc = conversion constant (lb mass) ft/(lb force) sec2. H = height of gas cavity, ft. (H, I, J, K, L, M, N) = coefficients of equation (68) ijk h = heat transfer coefficient = (heat transferred across space)/ (area of hot surface)(AT), Btu/hr ft. OF. hC = heat transfer coefficient at cold wall. hH = heat transfer coefficient at hot wall. he = local heat transfer coefficient. k = thermal conductivity, Btu/hr ft~F. x

kc = equivalent thermal conductivity including conduction and convection. k = mean value of thermal conductivity. m kw = thermal conductivity at wall conditions. L = width of gas cavity, ft. N = numerical solution to the difference equation. Nu = Nusselt Number for circular geometry = D (I ) NuL = Nusselt Number based on thickness of gas layer = hL/k. P = pressure, lb. force/ft.2 PD = pressure difference = P - Po P = hydrostatic pressure. P ijk= radius of eigenvalue bound circle (see equation 69). i jk Pe = Peclet Number = D4 Ig/a. p = tangential velocity, ft/sec. &Qi = radius of eigenvalue bound circle (see equation 70). ijk - q = radial velocity, ft/sec. R = radius of cylindrical space, ft. Re = Reynolds Number = D \fDg/v. r = distance from center of cylindrical space, ft. T = temperature. TH = temperature of hot surface. TC = temperature of cold surface. T = effective temperature of the surroundings where P = Po. t = time, sec. ui = ith directional velocity, ft/sec. V = total velocity vector, ft/sec. xi

xi = ith directional distance, ft. X = number of grid subdivisions in x direction. Y = number of grid subdivisions in y direction. a = thermal diffusivity = k/pCp. = volumetric expansion coefficient expansion coefficient, oC-1. 7 = c./T/(tx) A= exact solution to the difference equation (when used alone). Ass = steady state exact solution to a difference equation. AT = wall temperature difference, TH - TC0, C or OF. 8 = symbol for small error. Sx/By = a divided difference approximation of a derivative (see page = error. 0 = azimuth angle. = eigenvalue.,u = viscosity, lb mass/ft sec. v = kinematic viscosity, ft2/sec. p = density, # mass/ft.3 pO = density at T = To and P= P. PD = density difference = p - p 0 = angle between radius and temperature profile (see Figure 21). Embellishments T', p', etc. mean values a step ahead in time. T, p, etc. means dimensionless quantities, see equation (31). i, j, k directional indicators, see Figure 4 or 31. xii

ABSTRACT The velocity and temperature fields were measured at the center of a horizontal cylinder 4" in diameter and 36" long. The cylinder was filled with air at atmospheric pressure. The right half wall of the cylinder was heated and the left half wall was cooled. Wall temperature differences ranged from 20 to 2000C. The temperature field was measured by thermocouple traverses. The velocity field was measured by multiple exposure pictures of fine dust particles floating in the air space. It was found that all heat transferred by convection was accomplished by a rapidly rotating annular boundary layer. Inside this 1/2" thick boundary layer the central core remained essentially stationary due to a vertical temperature gradient produced by heat transfer from the inside edge of the boundary layer. As the wall temperature difference increased above 20'C, the overall rate of gas circulation decreased. The average Nusselt number based upon the cylinder diameter, the wall temperature difference and the average temperature gradient over either the hot or the cold surface was found to be constant at 7. The local Nusselt number was observed to be a definite function of the wall temperature difference at some points and essentially not a function of wall temperature difference at others. In addition to the above experimental investigation, an attempt was made to compute from the experimental boundary conditions the velocity and temperature fields for comparison with the experimental xiii

data. Five partial differential equations in five unknowns were developed which involve the conservation of heat, mass and momentum as well as the equation of state. These equations were reduced to three in number by assuming a perfect gas and essentially constant pressure. These three equations in three unknowns were converted to explicit difference equations for a cylindrical grid with variable radial mesh spacing (80 internal points in the grid). The calculation was done on the IBM-650 Magnetic Drum Computer. It was found that the above calculation procedure could be made stable only by abridging it so as to make it linear. For instance by specifying the velocity field obtained from experiment as constants in the numerical solution, one set of linear difference equations in temperature was obtained. The solution of this set of equations gave a good approximation of the temperature field obtained from experiment. It was also found that the largest stable time step for the linear equation could be computed from stability criteria for the heat equation since in this instance the fluid velocities have no effect on the allowable time step. The above new findings were explained in part by a matrix method of stability analysis. xiv

NATURAL CONVECTION INSIDE A HORIZONTAL CYLINDER INTRODUCTION In the world the amount of heat transferred by natural convection greatly exceeds that transferred by forced convection. Consider the heat conduction across the dead air space in a building wall, the ordinary steam radiator or other similar space heater, the circulation of air within a room, the winds and the weather. In all these examples natural convection heat transfer and the resultant convection currents play an important role. Nevertheless, natural convection has been a relatively neglected part of the science of heat transfer. As a result, the engineer who is charged with the responsibility of designing equipment to meet certain performance specifications naturally chooses to use forced convection heat transfer even when a simpler and cheaper design based upon natural convection may be possible. On the other hand, the engineer who is asked to predict effects such as room drafts or heat transfer relative to a hot or cold object in a room must rely heavily on his judgment. Sometimes natural convection heat transfer attains critical importance. For instance, nuclear reactors are usually cooled by forced convection but, in the event of a pump failure, the reactor must be shut down and the heat that continues to be produced must be removed by natural convection. There is no way to turn off the heat. Insufficient heat removal would cause melting and rupture of the fuel elements and result in an extremely serious accident. Consequently, the rate of heat transfer by natural convection must be known in advance. - 1 -

-2Natural convection heat transfer is also of critical engineering importance in the blade cooling of high performance gas turbines. It has been proposed that each turbine blade have a dead end channel drilled into it. All holes would open on to a large hole running through the axle of the turbine. Here, because of centrifugal force thousands of times stronger than gravity, natural convection becomes an efficient and simple cooling system. Another interesting application of the principles of natural convection heat and mass transfer which aroused the author's interest in the problem is a device developed by Busey and Hammond3 to dispose of a dangerous mixture of hydrogen, oxygen and gaseous fission products that comes off of a low temperature water-solution type nuclear reactor. In principle this device is a closed space in which natural convection currents are induced by the presence of a hot and a cold surface. The convection currents in the space dilute the gas from the nuclear reactor to non-explosive concentrations. A platinum catalyst applied to the hot surface causes the hydrogen and oxygen to react at any temperature or concentration. The heat of reaction maintains the temperature of the hot surface; the cold surface must be water cooled. Once the reliability and capacity of this device is established, the gases given off by the above water solution nuclear reactor will cease to be a troublesome liability and become a valuable source of neutron free gamma rays. The purpose of the present work is to learn more about the mechanism of convective heat transfer in a closed space. This purpose was carried out by observing the distribution of temperature and velocity distributions inside of a circular air space and by analyzing the experimental results to obtain an experimental correlation of the data. Also, a numerical solution of the partial differential equations that ex

-3press conservation of heat, mass and momentum was carried through to partial success in an attempt to duplicate experimental measurements for velocity and temperature measurements froa the surface temperatures, the diameter of the space and the property values of the air. This work was done using an IBM-650 Magnetic Drum Calculator. The thesis is divided into three parts: (1) Literature review and an explanation of the numerical analysis methods to be employed. (2) Experimental investigation including correlation and interpretation of results. (3) Numerical investigation including development of equations, analysis of stability, results of calculations and comparison with experiment.

PART I LITERATURE REVIEW -4

-5For the purpose of this thesis it seems appropriate to divide any review of previous work into two parts: (1) Experimental results and correlations related to the author's experimental work. (2) Methods of solution of the differential equations expressing conservation of energy, momentum and mass that describe the experimental system. HEAT TRANSFER CORRELATIONS Natural convection heat transfer coefficients for gas in a horizontal cylinder whose wall temperature varies around the circumference have not been found. Consequently this review will be concerned chiefly with vertical air spaces, i.e,s rectangular air spaces whose vertical walls are maintained at given temperatures and whose horizontal walls are insulated. One major objection in the comparison of all these data is the shape of the cavity. The following work was occasioned by the use of tall, vertical air spaces in the insulation of building walls. In addition rectangular shapes are simple both experime tally and theoretically. On the other hand, the author's geometry is circular but can be roughly approximated by a square cavity. Griffiths and Davis13 conducted some early experiments in connection with natural convection heat transfer across enclosed vertical air layers. They investigated cavity thicknesses from 2" to 1/2" and heights from 2' to 4' and temperature differences from 150 to 1000F. They found that equation (1) fitted their data quite well. h = 0.188 (T) (1) (For nomenclature, see page xi.) This equation has the same form as

the ronsrla given in erry's andbook29 for natural convection 6 -heat transfer to a vertical flat plate where the plate is in free air. However, the coefficient is 30% lower. Mull and Reiher20 did some careful work which has been recorrelated by Jakob16'17, employing a little bit of additional data (21 points in all). Jakob's correlation is given in equations (2) and (3). 1 -1/9 -c = 0.18 (GrL) (H/L) when 20,000 < GrL < 200,000 (2) kc = 0.65rL) (r) /(L) /9 when 200,000 < GrL < 10,000,000 (3) The generality implied in the above equations by the use of dimensionless groups is not completely substantiated. All data are for air at near ambient conditions. Also, there is no assurance that the equations will hold below an R/L of 3 because no dataare available below this ratio. Incidental to some other work, Peck, Fagan and Werlein28 measured the heat transfer coefficient between vertical metal plates 14 7/8" square which were enclosed by glass plates. For air at one atmosphere, they found that equations (4) or (5) represented the data well. NuL = 1,38 GrL 0.176 (/H)O0O406 (4) 0.18 h = 0.28 (~T)18 (5) The temperature difference between the plates varied from 150 to 1150F and the distance between the plates varied between 0.75 and 1.5". Equation (5) is based on experimental data only. Equation (4) is a generalization of equation (5) based on the theory of dimensionless groups.

-7Jones and Massonl8 performed a study of heat transfer from a horizontal cylinder and a prism to surrounding walls. The geometry that they used for the prism is shown in Figure 1. FIGURE 1 Geometry of Experiment by Jones and Massonl8 E Outside walls at room fT~ + /temperature Electrically heated prism They used an interferometer to measure the temperature gradients through the air cavity, and from the gradients they computed the quantity of heat transferred across the space by convection, Their study covered the following variables: L from 0.135" to 1.015" H from 2" to 4" (AT) from 500F to 112~F E from 0.125" to 1.866" If E is not small and if L > 0.025 ft., equation (6) is given as a correlation of the data. h = — + kw (0.37 tanh 36 L) (GrH)n (6) where

-8n = 0.25 - 0.0008L -1.43 (7) Robinson and Powlitch31 of the National Bureau of Standards have done some very careful work on the insulating value of air spaces. Their study covered the following variables as far as convection heat transfer is concerned: L from.62 to 1.69" H = 60" Width = 32" AT = not stated but small They found that by plotting L3 aT versus hL (h being evaluated at 500F), all data fell very close to one smooth curve as is shown in Figure 2. A very small correction is given for air space mean temperatures different from 500F. Equations (2) and (3), (5), (6) and the correlation on Figure 2 are compared directly with each other and with the results of this thesis on Figure 27. Finally, Ostroumov26 studied the temperature gradients set up in a horizontal cylindrical cavity filled with glycerine. This cavity was 35 mm in diameter and 35 mm long and was milled in a metal block. The ends of the cavity were enclosed with optically flat glass plates. Employing only one lens, a camera, a point source of light and a grid made up of vertical and horizontal rods, Ostroumov was able to obtain photographs where the shadows of the grid wires were made to represent lines of constant vertical and of constant horizontal gradient of refractive index in the glycerine. The refractive index gradient is, of course, related to the temperature gradient in the glycerine. The precise arrangement of the optical components is not apparent from the

II I III1 I 1 I I I IIII I 1 I il 0 N 0.7 HORIZONTAL.._ HEAT FLOW t 0.5 D 0.3 02. CONDUCTIQN ONLY 0.11 111111111I1111I 1I 1 10 100 1000 AT) L3, (OF)in3 Figure 2 Rate of Heat Transfer Across a Tall Rectangular Air Space, Robinson and Powlitch31

- 10 - description but the method appears to have promise for this type of investigation. No data or correlations were reported. MATHEMATICAL METHODS In order to avoid misunderstanding, the general form of the equations that govern natural convection as well as many other related phenomena are given below: Conservation of energy (a = constant) a =aV2T (8) Conservation of momentum (v = constant) "u-i G _ ge aP + + (9) - -- P xi 37xi Conservation of mass _ V goP) = (10) 6t Equation of State (for a perfect gas) P/p = KT (11) There will probably never be a general solution to the above equations. For many years mathematicians and engineers have been picking at the edges of this problem. Some have used the Analtical approach. That is, they simplify equations (8) to (11) by showing in advance that for the physical system under consideration certain equations or certain terms of the equations have negligible effect upon the solution. Then the exact solution to the simplified set of differential equations is obtained analytically. Others have used the Numerical approach. They do not need to simplify equations (8) to (11). However, approximations

-U - must be made because a grid must be constructed covering the space under consideration. The partial derivatives must be approximated in terms of the values of the variables at the grid points and in terms of the spacing of the grid points. Finally, considerable error may be incurred in solving the n algebraic equations in n unknowns. Neither of these two approaches to the problem is a panecea. However, in the future answers to more and more physical problems will be obtained using the numerical approach. In the following paragraphs the work of those investigators who solved approximations to equation (8) to (11) is reviewed. Analytical Approach Schmidt and Beckman33 presented an analytical solution of heat transfer by natural convection from a vertical flat plate which agreed well with their experimental data taken for air. They assumed that density was constant except as it perturbed the hydrostatic head. The effect of velocity on pressure was assumed to be negligible. They also assumed that the horizontal velocity has little effect and that vertical hesat conduction is negligible. All these assumptions are valid except near the bottom of the plate. Nevertheless, in this highly restricted case of a flat plate, a transformation to a new set of variables is necessary to obtain two ordinary differential equations which must be solved using infinite series involving eigenvalues. These eigenvalues depend upon the boundary conditions and are evaluated by trial and error. The computation and tabulation of the functions that constitute the solution are due to Pohlhausen33. Ostrach22 has redone the above solution in a more rigorous manner. In this same report Albers shows how to solve the two non

- 12 - linear ordinary differential equations of the Pohlhausen solution numerically on the IBM Card Programmed Electronic Calculator. The only parameter remaining in the equations is the Prandtl number. Albers checks the values given by PohLhausen for Pr = 0.72 and gives tabulations of the special functions needed in the solution for Pr = 0.01, 0.733, 1, 2, 10, 100 and 1000 in order to cover all possible fluids. Also, Ostrach23 analyzed the case of natural convection flow and heat transfer in the space between two vertical parallel long plane surfaces whose temperatures are maintained constant but not necessarily equal. The effect of heat sources in the fluid and of frictional and compressive heating were taken into account. Later Ostrach25 analyzed fully the above problem with the added complication that the walls may now have linearly varying temperatures provided both have equal slopes. He further included in his theory the case of mixed free and forced convection and the case of the fully enclosed space. He also treats the case of "convective inversion" where his modified Grashof number changes sign (i.e., where "gravity" changes direction or where the volumetric expansion coefficient changes sign or where the axial temperature gradient changes sign). These last two references involve exceedingly complicated methods of analysis and the results of this work are still only applicable to a relatively small class of practical geometries. But, for the parallel plate geometry and within the restrictions given above, all factors are taken into account. Lietzke19 checked experimentally the equations derived by Ostrach. Leitzke constructed a very thin, vertical annulus open to the air at the top and bottom. The annulus was thin enough so as to approximate parallel flat plates. The inside surface was heated electrically and the outside surface was cooled by water.

13 - Only the rate of heat transfer across the annular air space was observed by measuring the input of electrical energy and applying the proper corrections. Herman14 made some progress in solving equations (8) to (11) for natural convection outside a horizontal cylinder. He simplified equations (8) to (11) by using the concept of a boundary layer. Then by means of some clever transformations, was able to obtain two equations describing a measure of the flow and temperature fields as a function only of the Prandtl number. By doing this Herman was able to come to some useful conclusions about the phenomenon but was unable to solve the partially reduced equations directly. Batchelor1 treated theoretically the analytical approach to equation (8) to (11) for enclosed vertical air spaces like those found in building walls. For the case he covered that would most closely describe the experimental work of this thesis (i.e., H/L specified and Gr - co ), Batchelor reasoned that the convecting boundary layer would become small in comparison to dimensions of the core. He also reasoned that the core of gas inside the boundary layer would have a uniform temperature and therefore would have a flow pattern determined only by the drag of the boundary layer. For a circular geometry this reasoning would indicate solid wheel rotation of the core. Recently there has appeared a series of papers from the Molotov State University (U.S.S.R.)8'34'36 which outlined a method for approximate solution of the natural convection heat transfer problem using a power series of the Grashof number. The solution ca be made as exact as necessary by evaluating the coefficients of higher and higher powers of the Grashof number. These coefficients quickly become very complicated. They contain many terms involving the Prantdl number, the

temperature gradient, geometric ratios and property values. Zhukhovitskii36 has treated the case of a fluid in an infinitely long cylindrical cavity in a material having a horizontal temperature gradient perpendicular to the axis of the cylinder. He evaluated the fluid velocity and temperature and the wall temperature as a function of a power series in Gr. That is 7 = (l)Gr +a(2)Gr2 and T = T(O) + T(l)Gr + T(2)Gr2. The formulas for the coefficients of Gr (i.e., T(O), T(1), T(2), V(1) and (2))were given and contained the variables mentioned above. These formulas appear to be applicable in cases of high fluid viscosity or low temperature difference or both. It is claimed that the equations presented agree well with the experimental work of Ostroumov26. A variation of the above analytical approach is semi-empirical in that the shape of both the velocity and temperature profile is assumed on the basis of experiment. An excellent explanation of this method is given by Eckertll. Also Rutkowski32 has developed an independent method based upon using the error function as an approximation of the temperature profile. The use of the error function is claimed to have some basis in theory. The above two references refer to the solution of heat transfer from a vertical flat plate in air. In the above papers every mathematical trick has been used to obtain the analytical solution for natural convection near such elementary geometries as a vertical flat plate, a horizontal cylinder and two vertical infinite flat plates and for temperature boundary conditions equally idealized. It seems that further progress toward an approximatee nalytical solution for any arbitrary shape and for any arbitrary temperature distribution will be very slow indeed.

- 15 - Numerical Approach The numerical approach to the solution of partial differential equations holds promise of solving a great variety of problems heretoofore not amenable to the analytical approach. For many years it has been known that the basic partial differential equations can be approximated using finite difference equations and that a solution can be computed for any specific case. Using a desk calculator the Job, even for the simpler cases, was completely out of reason. Now high speed digital computers are available to do the massive amounts of computation involved. At the present time papers are appearing on how best to handle the solution of the heat equation by numerical methods. The heat equation is obtained from equation (8) by assuming that the velocity of the fluid is zero everywhere. (See equation (12)). Heat Equation, av2T = rT/at (12) This discussion will be restricted to two dimensional cartesian coordinates so oV = b2T/ x2 + 8r/6y2 (13) So in view of the complete lack of any literature on the solution of equations (8) to (11) by numerical means, the methods of solution of the heat equation will be discussed in some detail so that the reader will be better prepared to understand the attempts that were made to solve equations (8) to (11) numerically. Formulation of Difference Equation Solution Methods Now before a calculation can proceed, the boundary conditions

- 16 - for the space under consideration must be known. This knowledge must include the value of ca for the space and the value of T or its derivative all around the edge of the space. (T must be specified at at least one point.) Then a grid must be imagined to be placed over the space. The computation will find values of T only at the grid points. For example, if the space were rectangular the grid could be drawn as shown in Figure 3 and could represent a rectangular iron bar. Suppose now that when t < O, T for all grid points is zero. Then at t = O,T for the exterior points along the x and y axes becomes 100, while T for the remaining external points of the grid remains at O. The problem now is to compute the values of T for each internal point as the iron bar heats up and comes to equilibrium. Digital machines must approximate derivatives with divided differences. One has a choice how to take the differences. For instance, referring to Figure 4: T T_ Ti, J, k +1 TiJk, forward difference Ft 5 7 at = Ti, J, k + 1 - Ti,j,k-l, central difference It - Ti, j, k - Ti,j,k-1l backward difference At Of course the above definitions hold equally for space derivatives as well as time. It can also be shown from above that:;aT _ Ti, J+l,k + Ti,j.-l,k 2TiJ,k More complex representations of the second derivative using higher dif

-17y Y.. Figure 3 Sample of Grid for Heat Equation Calculation Example ijjq k t ijk i+ijk i, j-ik e~~~~q Frrqn:b.&.~.i~... qw I~ i Figure 4 Cartesian Co-ordinate Space-Time Nomenclature

- 18 - ferences are possible but are seldom used in calculations. Now there are a number of ways that these divided differences can be combined to represent the differential equation (12). They are: (Let Ax = Ar for simplicity) Tk k= +Tjik +T +T -4T (14) At (x Ti,j,k+ Ti,j,k- c 2[Tij+,,k+Ti,j-.,,k+Ti+,,j,k +Ti,,j,k- 4 Ti,j,k] (15) TiJ,k+,-Tijk =.a) [Ti,j+,,k++Ti,j-Ik+I+Ti+,,ik++Ti-,,j,k+r4Ti,jk (16) At (Ax)2 Ti j,k +Tijk = a [T,,] jk+Ti2T,,ir] ('7) At ax In the above equations, i and' can represent any internal point in Figure 2. For instance, if equation (14) is chosen to represent the differential equation (12), T at any grid point at time t + k is a function of a, ax, At, and T at the same grid point and at the four neighboring grid points at time t. The time step At is important because as will be shown abundantly later on,

- 19 - ht < (x)2/4a (19) in order to insure a convergent solution. Equation (15) has been proposed since the central divided difference is a better approximation to the derivative than is the forward divided difference employed in equation (14). However, it can be shown by the methods to be presented that equation (15) is unstable for all values of At. Equations (14) and (15) are said to be explicit because by using equation (14) or (15) quantities in future time (k + 1) can be calcuated by substituting into one equation known quantities in present time (k) and past time (k-l). Equation (16) is said to be implicit in both x and y directions because the Laplacian is written in terms of the temperatures in future time (k + 1). Obviously, each time step is accompanied by the solution of a linear system of equations, one equation for each interior grid point employed in the solution (in Figure 3, 28 equations in 28 unknowns would be required). However, it can be shown that there is no restriction on At due to stability considerations. However, large time steps will not be as accurate as smaller steps when the transient solution is desired. If only the steady state solution is desired, let At - co. That is, let the right hand side of equation (16) be set equal to zero. Equations (17) and (18) are recommended by Douglas and Peaceman7. They suggest that equation (17) and (18) be used alternately with the same time step as part of a double step. Notice that equation (17) is explicit in the y direction and implicit in the x direction. Equation (18) is the reverse. This scheme greatly cuts down the work because the solution of the simultaneous equations involving one row or column of points need be considered at one time and a simple method of calculating the points exists. At the same time the time step is not limited

-20 - by stability considerations. Types of Errors Involved in the Solution of Difference Equations Before proceeding further, some definitions are in order. Let, D. = exact solution to the differential equation A = exact solution to the difference equation obtained using a hypothetical machine carrying along an infinite number of digits N = the actual solution to the difference equation obtained using a real machine The main problem in numerical analysis is to keep (D-N) below scme practical maximum. But, (D-N) (D-A) + (A-N) The quantity (D-A) is called the truncation error and arises because a finite number of mesh points in both space and time are used. To find the condition under which A- D is a problem of convergence, the quantity (A-N) is called numerical or round off error. To find the conditions under which (A-N) is small throughout the entire region of integration is the problem of stability. Practically, of course, the solution N is employed when the solution D is unattainable so that the error, (D-N), cannot usually be computed directly, Now there are four fundamental questions that should be answered before machine calculation begins on the solution of a partial differential equation. They are: (1) Convergence -- Does the solution, A, to the difference equation approach the solution, D, to the differential

equation as the mesh size in time and space approaches zero in somne manner? The classic reference on this subject is by Courant, Fredricks and Lewy6 (2) Stability -- Do round off errors committed at each step of an iteration procedure tend to dampen out? (3) Time -- Which of several methods of solution will lead to the answer in the shortest amount of machine time? (4) Error -- What is the error in the final answer? If a direct method of solution is used, what is the accumulated round off error?. If an iteration method is used, how close has the solution approached to equilibrium? The answer to question 1 will not be treated here in detail because usually if the equation system is stable in answer to question 2, it is also convergent. Nevertheless, if equation (19) is not satisfied continuously as At and Ax approach zero, then the solution, A,, to equation (14) does not approach the solution D to equation (12). Thus compliance with equation (19) insures convergence but it is always obtained as a criteria for stability instead of convergence5'21'30, As the engineering aspects of perfonrming a computation are developed, more thought will be given to question 3. For instance, explicit procedures are relatively easy to program but because of stability considerations they may have an extremely long running time on the computer. The literature on new ways of solving equation (12) with less computing time is constantly increasing7 35. Question 4 is equally as important as question 2. Both questions will be treated more fully in the following paragraphs. Analysis of Stability.

- 22 - Stability may be studied by three procedures: (1) Consider a unit error introduced at an arbitrary mesh point and follow its progress. (2) Make a Fourier expansion of a line of errors and follow the progress of the general term of the expansion. (3) Determine the error matrix for the calculation and determine under what conditions the eigenvalues of the matrix are less than one. (1) Unit Errar The first procedure has been well worked out by Eddy12 Although Eddy's notation is difficult, the principle is not. Introduce a small error e at one point in the grid and develop an equation that will relate this error to all other variables including the time step. If the error magnitude eventually increases, then the calculation procedure is not stable. (2) Fourier Expansion The second procedure was developed and used by von Neuman. An explanation has been presented by 0'Brian, Hyman and Kaplan. To investigate the stability of equation (14), assume T to be the exact solution, A, to the difference equation. Then let T + 8 be the numerical solution, N, to the difference equation, where 8 is due to accumulated round off error. The quantity T + 5 must satisfy equation (14) as well as does T. Substitute T + 8 into equation (14) and subtract out equation (14) and obtain equation (20). A't (Akx)2 [5ii+l,k + 5ij-lk + 5i+lJ,k + 8i-1Jk (-48;,' ] (20)'

23 Because equation (14) is linear, it has the same form as equation (20). Let 5i,j,k = eat eibx eicy (21) which is merely a short hand notation for a Fourier Series. Consequently, i,j, k1 = ea(t*t)eibxeicy = eastijk and ai+l,j,k = eiCri,jk,,8i-.l,,k, e Yz,J,k etc. By substituting the above expressions into equation (20), one obtains, et 1 = e (+eax2 + -ib + eict + e ic44 (22) or e = 1 + (a (2cosbC + 2cos cy - 4) e + 1+(AX)2 In order that an arbitrary error introduced into the equation at t = 0 does not grow, it is merely necessary and sufficient that I eaAt -< 1. Thus for stability, the absolute value of the right hand side must be less than unity no matter what the frequencies and phase relations of b and c are. Therefore: At < (AX)2/4a (19) as was predicted previously. The statements mentioned previously about the stability of calculation methods using equations (14), (15), (16), or (17) and (18) can be verified by this method. An important refinement to the above theory is also given for

- 24 - grids with a small number of grid points. Equation (22) can be rewritten in the forn ea.tl,, ax)2t (cos b Ax + cos cAy - 2) or eaatm =4cs inF b2x + sin2 I] (22a) If only a smal number of mesh points are used only a small number of frequencies in the Fourier expansion exist. Thus if X is the number of subdivisions in the x direction and Y the number of subdivisions in the y direction, then sin2(bA x/2) < sin2 [ -l) /2 ] and sin2 (cAy/2) < sin [-1),/2Y]. For instance in the grid shown in Figure 3, X = 5 andY = 8, and equation (22a) becomes eeht.l = (Ax)2 (.90 +.961) or when leaatl < 1 then At (Ax)2/3.72a This refinement will be important in the explanation of the mathematical results of this thesis. Price and Slack30 also discovered a way of arriving at equation (19). Their method does not lend itself so easily to extension to equations more complicated than equation (12). (3) Matrix Method In the matrix method4 one seeks to construct a matrix such that when it is multiplied by a vector made up of all the errors in the iterated quantities at time t will produce an error vector at time t + At. If all the eigewnalues of said matrix are less than one, then the

- 25 - magnitude of the error vector decreases at each iteration and the calculation process is stable. Gershgorin's Theorem (sometimes referred to as the Hadamard Theorem)27 will be employed to give a bound to the eigenvalues without actually having to find them. The theorem states that the real (or complex) eigenvalues of a matrix A are located within circles whose center is at aii and whose radius is n n i laijl or lai I whichever is least, The symbol aij denotes an element of the matrix A. The error matrix will now be computed. Equation (20) for one space dimension only, is 8i,k+l 7- i+l,k + 78i-.,k + (1-27)bi8k (23) where y = att/(t&x)2. Let i = 0, 1, 2......n, k = 0, 1, 2, 3...... Assume that the boundary conditions have zero error, so 8i,k = 0 when i = 0 or n, of k = 0. This simplification to one space dimension is made in order that the structure of the error matrix to be presented can be written clearly on one piece of paper. For a series of points, equation (23) can be written as follows: 802 = 0 812 = (1-27)811 + 7821 822 = 7 11 + (1-27)821 + 7831 832 = 7 821 + (1-27)831 + 741 I N,,-~.,2- 78-2,. + (1m-27)%.1,.

- 26 - Which is equivalent to (I-2y) y,, T (I-2y) T 8a, ~2 y (I-2~) y 81 32 y (I-2y) x x = 2 I1 42 values for matrix B are less than or equal to one. Gershgorin's Theorem states that the eigenvalue is inside a circle with radius 2,y and center 1-2y. Figure 5 shows the situation. Clearly, a sufficient condition for stability is 2, ~ 1. For a two dimensional case a sufficient condition is 4r L 1 which is the same, again, as equation (19). FIGURE 5 Bound Circles for One-Dimensionl Heat Flow Imaginary Axis Circle Bounding Bigenvalues ircle Bounding Area of Stability,l l' Ei~ I>~ IReal Axis

- 27 - Of course the Fourier expansion or von Neuwn method and the matrix method lead to the same result for equation (12). For more complicated equations the stability criteria obtained by the two methods do not entirely agree. Analysis of Error The machine calculations performed for this thesis attempt to get at the steady state solution of equations (8) to (1l) by following a transient solution from an arbitrary set of initial conditions to equilibrium. Consequently, it would be useful to know within limits the steady state error, ess, where ess | 64S I(24) After extensive mathematical analysis, Collatz5 derived the simple mathematical formula for ess. He employed equally spaced grid points and a time step of one. Logically then =C Tk+l Tk max(25) Fss = (25) 6at where C is the maximum number of internal grid points in a row or a col1umn. If a direct method of solution were employed to solve the n equations in n unknowns, the analysis of error would be most important. A treatice on this type of error analysis is given by Hildebrandl5.

PART II EXPERIMENTAL INVESTIGATION -28

-29DESCRIPTION OF APPARATUS The Geometry In any study of natural convection heat transfer the geometry employed is, of course, the first consideration. It was decided to investigate natural convection inside a horizontal cylinder filled with air. This cylinder was to be long enough so that at the cente; flow in only two dimensions need be considered. In order to furnish a driving force for natural convection heat transfer, the surface of the cylinder was to be divided by a vertical plane into two parts that could be maintained at different temperatures. Figure 6 gives the dimensions of the space that was actually constructed. Desig of the E cuent A four inch standard iron pipe three feet long was cut in two lengthwise and welded into two boxes as is shown in Figure 7. The inside surface of the pipe was machined down to a constant diameter after the boxes were fabricated. Notice that thermocouple wells were introduced from the back side of the pipe and serve to stiffen the structure. Cooling coils were placed in the top of both boxes as well as were pressure taps and as were 1"11 pipe nipples for filling and emptying the boxes of the heat transfer liquids. Heat was supplied to the hot box by two 500 watt strip heaters (Chromolox SE 1850) bolted to the bottom. Notice that the hot box had steel cooling coils so that Dowthemn could be used as a heat transfer medium. Heat was supplied to the system by the electric heater. Energy input was measured by a watt meter and regulated by a variable transformer. The heat fron the electric heaters boils the heat transfer liquid (acetone, water, Dowtherni E, or Dowthenm A) that fills the bottom of the box

-30Tube for dust injection Lucite window waxed, in ploce and polished to fit pipe contour Asbestos Insulation Notes: I) Full scale 2) () Means thermocouple Figure 6 Geometry of Experiment

End View Pressure taps Side View ~ ~~ F 3.6" ~/coils Fill nozzle AL _ ~cooling n coil ~ _ ___,___,,_ _____ Thermo - couple ~ Section (-' well Section U-H''~'~~'"' Draw well, Thermocouple wells Is -- 12" Top View 0 "-'~Ii''...'' r / U0011 II. Section B-lB' |Channl for -lateral rChannel for'lateral illumination Figure 7 Plan of Test Object

-32-. to a depth of 2 to 3 inches. Before heating, the box is evacuated so that the vapor will readily condense on the inside surface of the box thus keeping the entire box at the saturation temperature of the vapors. A 30" mercury manometer aided in operating the experiment and in addition served as a safety valve. Some of the energy supplied to the hot box is transferred to the liquid in the cold box. It was found that for best results the box had to be filled almost full so that the cooling coils were immersed in the liquid. Originally it was planned to allow the liquid in the cold box to boil and to condense the vapors on the cooling coils but the above arrangement proved more satisfactory. Acetone was the original fluid employed in the cold side and proved satisfactory for all experiments. To minimize leaks, the cold side was maintained at several inches of mercury vacuum. The rate of heat removal through the cooling water was detected by measuring the flow rate with a calibrated flow meter and by observing the temperature of the liquid entering and leaving the cooling coil. Several hours were needed for the apparatus to approach steady state. Consequently, short term variations in line voltage and water pressure had no effect upon the experimental concditions and did not need to be tightly controlled. Nevertheless, the watt meter and the flow meter were examined often and long term variations in line voltage and water pressure were corrected for. Construction of Experiment The two boxes were positioned with respect to each other by using a cylindrical template whose OD matched the ID of the machined out pipe. The half inch space between the boxes above and below the pipe was filled with layers of wet asbestos felt insulation and allowed to dry. Then the unit of the two boxes and template was positioned onto a four

-33inch thick pad of Foamglas (Owens-Corning trademark). Four inches of this same insulation was used to surround the two iron boxes on all sides with the exception of the window at the center of the cold side and at the two ends and the space between the pipes at the top. The blocks of Foamglas were cemented together with flashing cement and covered with aluminum foil. On the top and around the cold side window coarse fiberglas insulation was used so it could be removed easily for the necessary repairs. The openings at the ends were closed with 5" x 7" glass Kodak Metallographic plates in order that the velocity and the temperature field could be observed inside the pipe. A diagram of the experimental setup is shown in Figure 8. Temperature Measurement All temperatures were observed with thermocouples employing a portable precision potentiometer*. In all cases a cold junction was used iltmersed in melting ice. For measurement of the surface temperature of the pipe and of other temperatures around the test object, 30 ga copperconstantan duplex thermocouple material manufactured by Leeds and Northrup was used. Before placement, all twenty thermocouples were calibrated against a standard mercury thermometer. The results of this calibration are given in Appendix D. Figure 9 shows the position of these thermocouples about the test object. Notice that thermocouples 1 to 5 and 11 to 15 are inserted into wells in the two boxes and in the inlet and outlet pipes of the two cooling coils. Thermocouples 6, 7, 9, 10, 16, 17, 19 and 20, were soldered to the boxes at the points indicated. Soft solder was used for the cold box and silver solder for the hot box. Thermocouple 8 was destroyed and 18 measures the air temperature outside the insulation. *Leeds and Northrup #8661.

Manometer byepass Rota meters Mercury manometers Thermo Hose con. to Alternate connection uple wells water supp Vent water supply To drain hose Vacuum pu mp Nozzel for filling Microswitches Lucite \ and emptying ight ure 8 Diagram of Experiental etup horn Cam turned_ Insulation by synThree way switch Batteries Flash bulb Electric heaters holder Figure 8 Diagram of Experimental Setup

(9) Cooling water in (5) Cooling water in hot side cold side ~ FiurCooling water inTemcul Cosition onTst cold side

-36The air temperature inside the pipe was measured by an independent thermocouple circuit made of chromel P-alumel. Forty gauge wire was used inside the pipe in order to disturb the flow as little as possible. Twentytwo gauge chromel P-alumel wire was used for the rest of the thermocouple circuit in order to decrease its resistance and thus enhance the sensitivity of the available potentiometer. This circuit as built was calibrated against a standard mercury thermometer. Over the range tested, the potential agreed very well with the published calibration of Leeds and Northrup. Consequently, the published calibration was used (see Appendix D for comparison). In contrast to all other thermojunctions used, this circuit was soft soldered instead of being welded. However, for the last two experiments the solder melted and a welded thermojunction had to be employed. The 40 ga. thermocouple wire with the junction at the center was stretched along the axis of the pipe. A length of welding rod fixed to slide in grooves in the cardboard end plate holders supported the thermocouple wire at each end. A small spring at one end of the stretched wire maintained the tension. Photographs of this arrangement with and without the glass cover plate are shown in Figure 10. Notice that the rod can slide in six grooves, three horizontal and three vertical, cut in the end plate holder. The end plate itself is a glass photographic plate. A piece of twenty squares per inch graph paper was placed over the emulsion when it was exposed so that a fine grid was produced on the glass. The two glass end plates were exposed fraom this one glass negative. Subsequently, as is shown in Figure 10, it was found that paper markers had to be glued to the end plates to point out what points were to be employed in taking data,

I-J(D Q O HO (D CD CD 0 CO) CD Ct CD P tC) CD

-38To position the thermocouple at a given point, the rods at both ends of the pipe had to be placed in the proper slots. Then each rod was moved by hand until the stretched wire pointed toward the gimn grid line. Figure 11 shows as closely as could be ascertained the position of the thermocouple in the central cross section of the pipe as related to the point number assigned for ease of data taking. Each temperature traverse of 17 to 19 points separated at a minimum distance of 1/20" required approximately ten minutes to complete. Thus the temperatures at the intersection of the different traverses were measured at different times and furnish information about the long term stability of the system. Velocity Measurement The velocity of the convection currents at the center of the pipe were observed by recording the movement of laterally illuminated dust particles suspended in the air. The dust employed was titanium dioxide pigment having a particle diameter of 3 to 7 microns. Using Stokes Law, this gives a settling rate of.05 to.2 inches/second. This is sufficiently below most of the convection current velocities measured in this series of experiments so that no routine correction was made for this effect. The particles were illuminated through a plexiglas window fitted flush with the inside surface of the pipe. A light horn shown in Figure 12 fabricated of plexiglas conducted the light from fast peak flash bulbs (Sylvania SM) into the center of the test space. Several sources of illumination were tried. The most satisfactory of Qiich was found to be a 250 watt mercury vapor lamp. However, preliminary experiments showed that the energy of the light itself was sufficient to cause a considerable convection current to be set up. This effect was unfortunate since the

-.39345,~0~ 97Q435.850 W/ _. _3590 -;/ 5 —1 2 -5. 7 _-6 69.50 15 7 9 113040 2 4 6 7 8 9 10 11 12 1416 1" —-- I". 1 4 6 i 9 1 0 9 14 1618,.00 35 8 9 10 11 12 13 15 119 \-277~ l" - I,10 _ 11 -10 r13. ~ 1,2 131517 230' 4 6 7 8 9 10 I 12 14i 250~ -12 13 12 159'50~ 214.50 -187 Notes: I) Fu II scale 2) Azimuth angles given for intersection of profile lines with surface Figure 11 Gas Thermocouple Location for Each Point Number

Commutator for Firing Flash Bulbs in Sequence Light Horn for Flash Bulbs Figure 12 Dust Particle Illumination System

-41120 cycle/second pulse of the light source was observable in the track of some of the more rapidly moving particles. Thus the time of exposure can be determined. For the reasons mentioned above, most of the velocity data were taken using the multiple flash technique. The flash bulbs were fired by microswitches riding on a wheel driven at 37.5 RPM by a syncronous motor. The wheel has a step in its rim. The flash bulb was fired when the lever of its microswitch fell into the step (see Figure 12). In order to calibrate the system, the wheel was painted black and a white strip was applied radially. A quadruple exposure picture of the wheel was taken by the light of the flash bulbs. The angular displacement between the images of the tape on the photographic film is then a measure of the time interval between flashes. The calibration used in the experiments is shown in Table I. One of the most serious drawbacks to this method as it was applied here was that reflection from walls obscured the dust tracks near the walls where, of course, they are of most interest. The walls were painted black and dusted with lamp black but this was not sufficient to completely eliminate the problem. What is required is a small slit on the TABLE I CALIBRATION OF FLASH CCEUTATOR Angular DisplaceFlash ment from First Time Interval frcm Number Flash Image First Flash 1 00 O sec. 2 8.0356 3 19.75.0878 4 32.25.l432

-42order of 1/161 wide rnuing around the entire circumference of the test space. A sheet of light would then be made to pass through this slit without any impingement upon the inside surface. Such a setup would require a complete redesign of the equipment and therefore was not attempted. With the unaided eye the tracks of the dust particles are discernable over the entire field. However, the recording of these tracks upon film required special techniques and then was not entirely successful. The camera employed was a Speed Graphic equiped for 4" x 5" film plate holders. The lens was a 1o. 32 Kodak Anastigmat F-4.5, 6-3/8" focal length with a solenoid tripped shutter. Most of the successful pictures were developed with freshly mixed X-500 developer made by the F R Corporation. The developing was done for the recommended amount of time. In order to take a picture of the flow pattern, the glass cover plates were replaced with opaque plates. One opaque plate had a hole in it centrally located with respect to the pipe. Into this hole was inserted a short black cardboard tube acting as a lens shade. This arrangement effectively cut out room light so that the shutter could be open for a considerable time without fogging the film. And, more important, it sealed the space from drafts that would have destroyed the flow pattern. It was observed very early that if the ends of the test object were left open to the air the dust blown into the illuminated central plain of the object would be dissipated in a few seconds starting at the bottom and working toward the top. If the ends were closed the dust would remain within the field of view for a matter of minutes. A small quantity of titanium dioxide pigment was blown into the test space at the center through a small tube extending vertically

-43from the top of the test object. During the 30 seconds it takes to reestablish the flow pattern checks were made to see that the film was loaded, the shutter was cocked and set on bulb, F-4.5, the flash bulbs were in place and the comutator was running. Then by pressing the button on the battery case, the shutter was opened and the flash bulbs fired. Releasing the button closed the shutter and permitted the camera to be removed. DATA Velocity Data As was explained in the Experimental Setup Section, the primary data for velocity is a multiple exposure picture. (For example, see Figure 13). This picture was analyzed by tracing out the flow lines and by measuring the distance between the first and last image of the same dust particle to get a measure of the magnitude of velocity. The velocities are represented at the horizontal and verticle diameters in Figure 14 by by vectors. The heavy line connects the points of these vectors. Flow tracings for all experiments for which data was collected are given in Appendix A. Temperature Data The primary data for temperature is in the form of velocity traverses that were taken over a period of time. These data appear in Appendix B. The temperature data for each experiment was plotted up as show in Figures 15 and 16. The letter designations for the various profiles in Figures 15 and 16 refer to the nomenclature given in Figure 11. In order to get a better idea of the temperature surface expressed by the profiles in Figures 15 and b16, a three-dimensional model was constructed as shown in Figure 17.

~ 11:'i.!:2,~ ~~~~~~~~~~I~ii WC0009:~~~~~~~~~~~~~~~~~~~~~~: i~!:i~i'!:::'.':........i: X,~!!i L -:4. ~11iu 13 Sp Ds Particle- ~t i~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~i~ ~'<::.:;7.'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ei-;~~~~~~~~~~~~~~~~~~~~~~ i.... -; ~.~,~,.~il ia pl,;2:':"c e ic ur

Flow profile Flow lines Scale for velocity vectors which are drawn from diameter to flow profile L,. I I 0 5 in/sec Figure 14 Experimental Flow Pattern (Exp. 4)

70 60 0 50 -2 -I 0 +1 +2 DISTANCE TO RIGHT OF VERTICAL DIAMETER (IN) Figure 15 Horizontal Gas Temperature Profiles (Exp. 4)

l I" I, I O F 60 0 30 -2 -I 0 2 DISTANCE ABOVE HORIZONTAL DIAMETER (IN) Figure 16 Vertical Gas Temperature Profiles (Exp. 4)

........ =:--;~;;~ —,,_ _,-r;;:::Li;:;. —-—: —; -.-. ~~; i- Z-z ,; _5 - =1 - -"".;;~i; 1 - I; ~ I;; iji $ I iii i s% I IlllglBPB~EBi IC I ~-~r I, i - i~, I,1 I i. — ---~-r iiL Bsaasi -sa ~1 03 I n f I r 4 --— l_.i~-i~i)i~li-lii~li-~~ LBig -PPall-ih —~La — -q " - I bSb Figure 17 Three Dimensional Representation of Temperature Surface Observed in Experiment 4

ANALYSIS OF DATA Velocity C orrelations Only incomplete measurements of velocity were possible with the experimental equipment. (See Appendix A) An attempt was made to organize this Limited information along two lines; (1) a correlation of local velocity with the wall temperature difference, (2) correlation of the rate of circulation through a given radius with the wall temperature difference. Local Velocity Correlation Some local velocities were measured as shown in Table II. The only fairly complete set of localielocities is near the cold side of the horizontal diameter where the illumination is the most intense. One set, however, was obtained along the vertical diameter. These are shown plotted against wall temperature difference in Figure 18. On the horizontal diameter the scatter of data is rather large, however, the trend of the data seems to indicate that all local velocities increase slightly or stay constant until a wall temperature difference of 200C is reached. Then all local velocities decrease as the wall temperature difference increases. This observation is contrary to expectation. One would expect that the thickness of the annular velocity boundary layer would decrease and thus cause the local velocity at points far from the surface to decrease and at points near to the surface to increase. The thickness of the velocity boundary layer does decrease (see Figure 19). However, judgement should be reserved on the decrease in all local velocities with increasing temperature difference until more reliable data are available for velocities close to the wall.

TABLE II LOX)CAL EXPERIMENiTAL VELOCI TIES Exp. On Horizontal Cold Side Radius On Top. Vert. Rad. No. =*4.42 =.44 P.46 | r.48 T = 30T 4 O"/sec -.75* -2.2 -3.4 - -.48 38.90C 6 -1.25 -2.15 -3.00 -3.85 - -.40 20.7 7 -1.30 1 _ _ -.53 5.72 8 -1.20 -1.55 -1.80 -1.90 -1.80 -.51 1.93' 9 -.45 -1.02 -3.1m4 -3.58 -2.32 -.31 51.58 10 -.20 -1.38 -2.87 -2.22 -.88 - 90.5 11 0 -.58 -1.38 -2.59 -.26 75.2 12 0 -.45 -1.95 - - +.36 110.52 13 - - - +.30 111.24 14 0.43 1.20 - - + 50 174.10 * for velocity nomenclature see Figure 26.

-51I I I1 1i1 1 1 1 1 1 111 Vertical radius -top 100 Horizontal radius u cold side I-42 w Z T - =.30 + Vertical radius, top Note: _ Distance from Center Diameter of circle 2 10 100 200 Figure 18 Correlation of Sele C) Figure 18 Correlation of Selected Local Velocities

Exp 8, AT= 1.930C Exp 9, AT=51.58 ~C Exp IO,AT=90.5~C 70 - n ~~~~~~~~~~~60< 4~~~~~4 w 5050 r w 30 40 o 0 2 O0 ~I a- l CL C) -3 LL~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -4 _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ i,!,, I I I I I ~0 0.5 I 0 0.5 I 0 0.5 DISTANCE FROM WALL (IN) Figure 19 Velocity and Temperature Profiles for Horizontal Radius at Cold Side (Expts. 8, 9, and 10)

Rate of Circulation Correlation In order to keep from accumulating mass in one part of the space at the expense of another, the net mass rate of flow through any diameter must be zero. Because of incomplete data, however, only the mass rate of flow in the boundary layer at the center of the cold side were computed for experiments 8, 9 and 10. In only these three experiments was the maximum velocity in the boundary layer definitely measured. These velocity profiles shown in Figure 19 were integrated taking into consideration the temperature profiles also shown on Figure 19 to give a net rate of circulation of the boundary layer for a unit length of cylinder in terms of cubic inches of air at 200~C per second. The results are shown in Figure 20. Notice that circulation rate decreases as the wall temperature difference increases. Temperature Correlations The temperature correlations are based on an analysis of the temperature profile data that is given in Appendix B. For experiment 4 these data are graphed in Figures 15 and lb. Method of Computation The temperature profile plots similar to Figure 15 and 16 were examined for all experiments. It was found that for almost all experiments, three to five experimental points near the surface lay on a reasonably straight line. The slope of a straight line through these points as well as the confidence limit of the slope was computed by the method of least squares9. The points that were chosen are the same for every experiment and are tabulated in Table III according to the nomenclature show in Figure 21. The point nearest the surface at slope numbers 2 and 8 was discarded when computing the slope because its reading

I I I I At temperature of boundary layer LI Normalized to 200C z 0 j Note: Rote of circulation of boundary layer O measured here LLU0 0.5 ELD Tc i )TH IT TH 0 I I I I I I j I I I 0 50 100 AT (~C) Figure 20 Correlation of Boundary Layer Circulation

TABLE III GECiMETRIC DATA Slope Points Slope Correc. Weighting Number Location Used Factor Factor 1 Profile E* Top 1,2,3 * 1.000.07778 2 Profile D Top 2,3,4,5.887.21111 3 P rofllo A Cold Side 1,2,3,4.887.17361 4 Profile B Cold Side 1,2,3,4 1.000.314861 5 Profile C Cold Side 1,2,3.895.17361, 6 Profile D Bottom 13,14,15,.887.21528 16,17 7 Profile E Bottom 17,18,19 1.000.07778 8 Profile F Bottom 13,14,15,.887.21250 16 9 Profile C Hot Side |13,14,15 |.891 | 17361 16,17 10 Profile B Hot Side 16,17,18,19 1.000.15000 11 Profile A Hot Sida 14,15,16,17.887.17222 12 Profile F Top 1,2,3,4,5.887.21389 * See Figure 10

-56Line of zero temp. grad.lope no. I-b%.igr 1Sop ubrAragmn:A Figure 21 Slope Number Arrangement

was consistently off, indicating wall interference. As is seen in Figure 21, the temperature profiles do not always intersect the surface at 900~ The temperature gradient can be expressed as follows: VT * T + _ a T ar r d8 If a constant wall temperature is assumed near the wall, then near the wall I dTI 1 dT I4a>> I- 4I Thus for slope 3 in Figure 21, XVT aaT OT _ aT V T ~ m a r - m cos 4'3. The cos ~ is the slope correctionfacet given in Table III. From the normalized slope measurements the Nusselt number, Nu, is obtained by multiplying by D/AT. The rate at which heat is transferred by convection is computed by multiplying the normnalized slope measurements by the thermal conductivity of the gas at the wall assuming that the wall is at either the hot or cold temperature for the experiment. The heat transfer coefficients were obtained by dividing the rate of heat transfer byLT. All the above quantities refer to specific positions on the surface. In order to compute average quantities for the hot and for the cold surface, the normalized slope measurements were multiplied by their corresponding weighting factor and summed. The weighting factors are proportional to the arcs shown for the slope numbers in Figure 21 and are adjusted so that the weighting factors for slope numbers 1 to 6 which are in the cold side add up to one and the factors for slope numbers 7 to 12 which are in the hot side add up to one also. From the averaged normalize slope, or temperature gradient over the hot side and over the cold side

-58the Nusselt number, the rate of heat transfer and the heat transfer coefficient were computed as outlined above. The confidence limits on all the quantities above are a very important part of the analysis and warrant some special discussion. The 540 confidence limit was computed as outlined by Duncan9. Thus, within the range given by the confidence limit the true slope will lie 504 of the time providing the probable error in all experimental points is equal. The fact that only 3 to 5 experimental points were used has been taken into account in the calculation. iven the temperature gradients are averaged to get the effective temperature gradient for the hot and cold side, a truncation error enters in because only a small number of temperature gradients were observed. This truncation error is not included in the confidence limit computed for the averaged values. The theory given by DuncanlO for the sums of approximate quantities was employed. The results of the above outlined calculations for the case of experiment 4 is given in Table IV. The same table for all the other experiments is given in Appendix B. In Table IV and in some of the following graphs, negative values indicate that the wall is heating the fluid; positive values indicate that the fluid is heating the wall. Local Nusselt Number Correlation For the local data, the Nusselt number is probably the most usable since it tends to bring the data for all wall temperature differences into the same range of values. Figure 22 gives Nusselt numbers for experiment 4 with their confidence limits plotted against the azimuth angle as defined in Figure 11. The information shown on Figure 22 could be given for all twelve experiments in one graph. However, because of scatter in the

TA:B iS IV ANALYSLS OF TEMPERATURE DATA FOR EXPERINENT 4 Slope Normalized 50%Nusselt No. 50% Conf. Rate of 50% Conf. Heat. Trans. Conf Slope Nonnalized ileat. Trans Linit Coefficint Limit No. Slope ~C/in limit OC/in Nu = D t Heat miTrans Limit oeftfi C __ r by Cony. BTU/hr ft2 TU/hrft OFBU/hr ft2OF..T SrI:',/rf 2 1 54.00 + 30.02 5.98 + 3.32 17.61 + 9.79 0.252 + 0.140 2 99.79 + 3.41 11.04 + 0.37 32.55 + 1.11 0.465 + 0.016 3 81.40 + 3.01 9.01 + 0.33 26.55 + 0.98 0.379 +.021h 4 76.40 + 3.49 8.45 + 0.39 24.92 + 1.14 0.356 + 0.016 5 71.51 + 1.29 7.91 + O.14 23.32 + 0.42 0.333 + 0.006 6 34.61 + 0.45 3.83 + 0.05 11.28 + 0.15 0.161 + 0.002 Ave. C.S. 70.62 + 1.71 7.82 + 0.19 23.03 + 0.56 0.328 + 0.008 7 - 48.00 + 2.31 - 5.31 + 0.2o - 17.47 + 0.84 - 0.249 + 0.012 8 - 114.93 + 5.75 - 12.72 + 0.64 - 41.83 + 2.09 - 0.597 + 0.030 9 1- 9.91 + 4.69 - 10.50 + 0.52 - 34.54 + 1.71 - 0.493 + 0.024 10 - 100.4 + 12.69 - 11.11 + 1.40 - 36.54 + 4.61 - 0.522 + 0.066 11 - 70.57 + 4.59 - 7.81 + 0.51 - 25.69 + 1.67 - 0.367 + 0.024 12 - 33.46 + 1.06 - 3.70 + 0.12 - 12.18 + 0.39 - 0.174 + 0.006 ve. H.S. - 79.00 + 1.98 - 8.74 + 0.22 - 28.75 + 0.72 - 0.411 + 0.010 A T = 70.02 OF / 0110668 /C k' 0.32616 BTU 0.36396 C HR FTZoc/in H HR FT2 oC/in

10 i | | \' [- Hot side z w 2 z 0 1210 I 240 Hu 2LAZIMUTH ANGLE,8 (~) i -2 con | - Asbestus insulation /, 0 Z -4 //////////__ o t -6 L L Pipe wall o -6 - 8 tCold side -10 50/ % conf. Figure 22 Local Nusselt Number as a Function of Azimuth Angle (Experiment 4)

-61-. data it is not clear whether the local i4usselt number is a function the wall temperature difference or whether it is not. Figure 23 shows the local Nusselt numbers with their 50% confidence limits plotted against the wall temperature difference AT. Since the cold wall was maintained at substantially the same temperature for all experiments, the temperature difference alone will detenrmine the property values of the gas at a given point. (Note: According to the data given in Appendix B the metal walls are substantially isothermal having an extreme variation at AT = 2000C of 250). Consequently, Pe, Re, Pr, Gr and AT/T which are used to characterize natural convection are in this experimental series only a function of AT. Nevertheless the scatter of data in Figure 23 is quite large. This scatter is not the result of experimental error in measuring the temperature gradients because in most cases the 50% confidence limit is negligably small. Above AT = 100C it is reasonable to assume that the scatter of data is random due to a large number of unknown causes. Below AT = 100C a different mode of heat transfer may begin to function but the two experivnental runs made in this region are not enough to define it. Somewhat arbitrary lines were drawn through the data in Figure 23 to indicate trends. These smoothed results are shown in Figure 24 in the same form as Figure 22. A comnparison of Figure 24 with the original data showed that the main characteristics of the data were preserved. Overall Correlations For experiments 3 to 15 the overall convective heat transfer data is summarized in Table V. These data can be plotted in a number of ways. Figure 25 shows the rate of heat transfer to the cold side as

-62COLD SIDE HOT SIDE A T,(C) A T,(~C) 10 100 I0 100 1000 - I I Irill -- i I liliIl I I 111l I Il lilil I I ill I I I 1lll +10'_____ _ ~Slope 7 +10 | Slope I - Slope 2 Slope 8 +10l -I10 2,+5 n I > I Iv IS J -15 Slope 3 I Slope 9 z +105 r +5 zo w cn w o3 c __ Slope 4 Slope 12 5 10 S I0 o " +5 Slope 5 +10 0' ~I I I~a Slope II Figure 23 Smoothing 6 f Local NusseltSlope 12..0Slope 6 0 EXPERIMENT NUMBER Figure 23 Smoothing of Local Nusselt Numbers

12 i ~~ I oA= 2000C 20 W - Z 0 (I) D SLOPE NUMBERS-J 0 0 0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~0 U)~ ~ ~ ~~~~~~~~~~~~~~2 w 3 005 AT 2000C IJ~ ~ ~~~~ 120 240 AZIMUTH ANGLE.,O (0) Figure 24g Smoothed Loal Nusselt Number as a Function of Azimuth Angle

TABLE V SUMMAIY OF AVERAGED MEASURES OF HEAT TRANSFER RATE Exp. TC TH 4T Average Normnnalized Overall Nusselt Rate of Heat Trans. Heat Trans. Coeff, I No|. o 0T Temp. Grad. ~C/in Number BTU/HR FT2 BTU/HR FT2OF OC OC O C Cold S. Hot Sd S Hot To Cold SS. Hot To Cold S. 3 1 41.5 65.5 24.0 33.60~ 1.42 41.65~ 1.31 6.03+ 0.25 7.47~ 0.24. 11.50+ 0,.49 15.20+0 0.266+ 011 0.362_0.011 4 24.5 63. 4 38.9 70.62+1.71 79.00+1.98 7.82~0.19. 8.74+0.22 23.03+0.56 28.75 0.72 0.328 40008 10413+0.010 5 16.5 42.2 25.7 35.31~ 0.88 410.81 + 0.90 5.92 + 0.15 6.83 + 0.35 11.29 +0.28 14.02 04311 0.2 44 006 0.303.007 6 19.5 40.2 20.7 26.05 0.8l 29.48 0.57 5 42+0.17 6.13,0.12 8.36 ~0.27 10.06 PJ9 0.224 007 I.270-+o05 7 17.12 22.84 5.72 4.74+ 0.42 5.25 +0.28 3.57 +0.32 3.95 +0.21 1.51 +0,13 1.70+0,. 09O.147+.014 0165_0.009 8 16.20 18.13 1.93 2.72+0.38 3.07 +0.10 6. 0+0.86 6.85+0.23 0.86 +0.12 0.98+00310.2490.035 0.2830O.010 9 20.62 72.47 51.85 63.20+ 2.501 85.07 1.50 5.25 +0.21 7.06 0.12 20.33 0.80 31.7Q0_O.561 0.2180.008 10.30.006 10 26.10 116.62 90.50 165.53~5.36370.84+2.57 17.87 ~0.25 8.13 +0.12 54.17 ~1.76 73.25~1.10 0.332_0.010 0.9440.007 11 25.00 100.20 75.2 125.78~4.97 6.97 ~4.02 7.20+0.28 8.41+0.23 41.02+1.62 58.73~1.60 0.303+0.012 0.34+0.012 12 36.80 147.32 10.52 148.76+6.1180.79+3.70 5.79.+0.24 |7.04 +0.14 50.29+ 2.08 87.28~1.78 0.253+0.010 0.4390. 13 35.48 146.72,24 155.74~5+.6 9.48 ~3.58 6.03+ 0.22 7.33~+0.14 52.20~1.90 91.2 61,72 0.26310. 100456.009 14 46.135 220. 25 f.10 239.72~ 8.622 3.78+ 7. 64 5.93 ~0.21 7.26t+0.19 83.36 ~3.00 P79.59-+4.67 0.2660.C10 0 573+0.015 15 50.45 254.10.6 1 288.3110.18131.30~ 8.451 6.09 ~0.22 | 7.32~ 0.18 100.42+ 3.54 233.30_5.67 0.274+0.0100. 637-015

RATE OF HEAT TRANSFER (BTU/HR FT2) 0' i I' I I I R) ~dD CD (DD H' ~-3 F- \ c+ ('D" iCD CD ct Fj- I-3 ob a ~(D Fl ~ - o 0 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~= —= (D CD00 __ I> o -- r~ D ~~~~~~~~~~~~~~~~~~~C 0.......... I I II FO ~ ~ ~ ~ I I, III!

-66a function of the wall temperature difference. All data were taken for air at one atmnosphere inside a 4,S3" ID pipe, nevertheless the results are to be plotted as Niu - f(Gr x Pr) (See Figure 26) as is the custom. If all twelve experiments are considerecdNu $ f(PrGr), although there is some scatter of data. However, if only experiments 3 to 7 are considered, a tight correlation is obtained. On consulting the records it was found that a leak occurred after experiment 7. This leak was in the cooling water lines at the top of the test object and many days were required to dry the apparatus out. At this time also the cylinder was repainted and dusted with carbon black. There is a possibility that this accident changed the characteristics of the apparatus. It is shown in the Discussion of Dimensionless Groupsin Part III that Gr(Pr). Pe(Re) LT/T. Thus the three dimensionless groups that govern heat transfer are usually all treated alike. If data were available in which Pe, Re and AT/T are varied independently over a good range then the effect of each group could be determined. For experiments 3 to 15 the dimensionless groups have the values listed in Table VI as well as Table V. Comparison with Previous Work As is explained in Part I, all the previous experimental correlations of the coefficients of convective heat transfer in a closed space have been for air near ambient conditions. The geometry most comnparable with the present work is the tall rectangular air space with heated and cooled vertical walls and insulated horizontal walls. The previous investigators usually report their correlation of the heat transfer coefficient in terms of AT, H/L, and the property values of the gas. For H/L = 1 and for ambient air, the correlation of Peck, Fagan and ierlein28

10 9 4/ 8 —1 I I AIA I 12 exp /3/ = / 3-7 c Hot side.01.1 I 10 (Pr)(Gr) Figure 26 Dimensionless Heat Transfer Correlation A Exp~ /, oCodse 3-7 incl A Hot side'Figure 26 Dimensionless Heat Transf er Correlation

T~.A.BLIE V I DDIMENSIONLESS GROUPS FOR EXPERIMENTS 3-15 ExP a x 104 ~ v x 104 Pa x lob Re x lob Pr Gr x lO8 AT N. ft2/sec ft2/sec T T C OK HK Cod S Hot S. Cold S. Hot S. Cold S. Hot S. Cold S. Hot S. Cold S. Hot S. Cold S., Hot S. 3 2.70 3.10 1.90 2.16 2.22 2.64 1.56 1.77.70~.697.186.221.07.64 0 4 2.44 3.06 1.72 2.lb4 2.00 2.51 1.41 1.76.70o5.700.260).356 4,1309.15 5 2.32 2.71 1.65 1.91 1.90 2.22 1.35 1.57.710.705.161.200.0887.01 6 2.36 2.67 1.67 1.88 1.94 2.19 1.*37 1.54.707.704.132.156.0707 7 2.33 2.42 1.65 1.70 1.91 1.98 1.35 1.40.708.702.0357.0378.0197.13 8 2.32 2.35 1.65 1.67 1.90 1.93 1.35 1.37.710,710.0121.O124,oo067.06 9 2.40 3.23 1.68 2.25 1.97 2.65 1.38 1.85.700.697.336.512.177.5 10 2.47 4.36 1.74 3.00 2.03 3.58 1.43 2.46.705.688.615 1.350.302.3 11 2.45 3.72 1.73 2.57 2..01 3.06 1.42 2.11.706.69o,508 Oe8895.252 12 2.62 5.55 1.085 3.77 2.15 4.56 1.52 3.09.706 0680.820 2.50.356.6 13 2.60 5.52 1.84 3.76 2.1i4 4.53 1.51 3.08.707.682.820 2,50.360.6 14 2.75 8.30 1.95 6.00 2.26 6.81.1.60 4.93.709.723 1,390 8.57.545.5 15 12.83 19.64 1.9 7.10 2.3 7.90 1.6 5.. 703, 737 1,1635 113.10.618.8 Note: D = 4.s305"1 g 32.16 ft./see.2

-69isP 0.18 h 0=. 28 ( AT) (5) The Griffiths and Davis13 correlation is, h ~ 0.188 ( A T)'25 (1) The J. kob 17 correlation is, h.0945 ( AT)1/3 (3a) The correlation derived by Robinson and Powlitchi is not formulated but is given by Figure 2. The correlations mentioned above are plotted on Figure 27 along with the experimental values for h. and hH taken from Table V. A line was drawn by eye through the experimental points. Because of the scatter of the experimental points this correlation must be regarded as preliminary. Nevertheless, it has the following fomnnulas: hH = 0.115 ( AT)0'25 for 25 < AT < 4000F hH = 0.27 for 3 < AT < 25OF (26) hc a 0.27 for 3 < AT< 400oF Notice that the present work gives heat transfer coefficients of the same magnitude as have been measured previously for somewhat similar geometries. However, for the cold side at least it does not show any real dependence on the wall temperature difference. It is difficult to say why the heat transfer coefficient observed in this work is practically constant. Possibly the circular shape of the space surpresses the fonnation of eddies in the boundary layer. Conclusions The experimental data and the analysis of it presented in this section is exploratory in nature. Even so a few conclusions can be drawn from the analysis.

h For ambient air at H/L=I [0 0~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~0 CMNo 0.2 Equa. 5 IZ ~~~~~~~~~~~~~~~~~~~~~~~This wcr 010 10060 M 0 ao n~~~~~~~~~~~~~ 0.2 Fg Equa. 3a 0.1 2 ~~~~~~10 100 0 AT (0C) Figure 27 Comparison of Experimental Heat Transfer Correlations

-71(1) The local fluid velocities and the rate of circulation in cubic inches per second decrease as the wall temperature difference increases above 20~C. This conclusion is not absolutely certain and is subject to future verification. (2) For this series of experiments the local Nusselt number is sometimes dependent upon the wall temperature difference and sometimes is not. A smoothed representation of this rather complex relationship is shown best in Figure 24. (3) The overall convective heat transfer coefficient for air at 1 atm pressure and near ambient temperature in a 4.3" ID horizontal cylinder that is divided by a vertical plane through the axis into a heated and cooled semi-circle is 0.27 BTU/hr ft2OF based upon the wall temperature difference and the cooling surface area, if the wall tenmperature difference is less than 25~0F. For 25 < AT < 400OF, hH 0.115 ( AT) 25 and he remains 0,27.

PART III NUMERICAL INVESTIGATION -72

-73Just as Part II includes sections on how the work was accomplished, what work was done, and what significance the work has, so also will Part III contain the same sections. Part III is an account of only some of the more successful and relevant calculations that were made in an attempt to solve the differential equations that describe natural convection (see equations (8) to (1l)). The goal of these calculations was to compute the temperature and velocity fields corresponding to the boundary conditions of an actual experiment and compare these fields with their experimentally measured counterparts. Several steps are necessary in the solution of a differential equation system by numerical analysis. They are: (1) Formulate the differential equations and the boundary conditions that describe the physical system under consideration. (2) Convert the differential equations to difference equations and formulate a method of solution. (3) Analyze the method of solution for accumulated round off error if a direct method is employed such as Gauss elimination. Analyze for stability if an iteration method is employed. (4) Pick a grid and, if necessary, a time step in accordance with the above analysis. Try to estimate the truncation error that will be incurred by using this approximation for a continuously varying function. (5) Code the problem for the available digital computer. Make an estimate of the machine time required. (6) "Debug" the problem. This task is the most time consuming and vexing part of entire process.

-74(7) Take the data. Naturally a large amount of data would have to be required to make the above process economically sound. Beutler2 has treated the economics of computation. Eighty hours of machine time on an IBM-650 Magnetic Drum Calculator were used to produce the work described herein. FORMULATION OF THE DIFFERENrIAL EQUATIONS Basic Equations Since the object of these calculations is to reproduce experimental data, it is extremely advantageous to work with polar coordinates. The nomenclature that will be used is shown in Figure 28. FIGURE 28 POLAR COORDINATE SYSTEM T TH Down In polar coordinates, equation (8) becomes,

-75a + - + q = a (v2T) (27) equation (9) becomes, + r + q r= - g sin @ + (28) v I pr T o(28) [v p +1 ~ )] and also, + + 6 + q r= gcos C _ + (29) v V2q + 6 (V)] and equation (10) becomes,,(PO) + 6r (=. t- (30) Theoretically at least, it is possible to solve any heat transfer or fluid flow problem with the above equations. Only in the most elementary cases have enough approximations been allowable so that the equations could be solved analytically (see Ref. 8,14,22,23,24,25,33,34,36). However, with the availability of high speed computing machines the general set of partial differential equations will eventually be solved by numerical methods. Since a machine solution may be expensive to arrive arrive t, it is advantageous to consolidate the variables so that one calculation will be applicable to a large class of real situations. The variables in

-76equations (27) to (30) will be consolidated by non-dimensionalizing them. Let, T = T/Tef P/Pref P = P/Pref q r - q =I p =, D, t t ) (31) If the reference temperature, Tref, and the reference pressure, Pref, are specified, the reference density, Pref, is fixed by the equation of state. The reference quantities are entirely arbitrary. Also equaticns(32) and (33) will be simplified by assuming only mild variations in density so that, V2q >>l3 r ] (32) and V2P > >3r a (33) Therefore, neglect the second term in parenthesis in equations (28) and (29) and substitute the transformations (31) into equations (27) and (30). The results are: 65 6P_ drapl Eu rP + P_ on+q S_ sin 9 (34) +Re.. 62j + 6- 6- E6 Eu a E + + ~ + 4 q =- Cos _- + (35) Re - + 71

-77Fr +i g + q a+I = 1, pl; + _; 1 r a; (36) rS + + FE = - (37) where Eu = g Pref/Dg Pref' Re = D 405/v, Pe = D,g/ia. Thus only Re, Pe, and Eu need be specified in order to describe the system. Notice that equations (34) to (37) constitute a system of four equations in five unknowns, p, q, T, p, and P. The equation of state relates T, P and p and gives a soluable system of five equations in five unknowns. A program based directly upon these equations was written for the IBM-650 and "debugged". It employed an explicit iteration method to arrive at the steady state solution. However, trial runs of the program as well as analysis of stability showed that the time step was much too small to be useful if one in fact does exist. Discussion of Dimensionless Groups Notice that if the Reynolds number, Re, the Peclet number, Pe, and the Euler number, Eu, are specified along with the boundary conditions for T, p and q, then the system is fully described. It has become the custom in studies of natural convection heat transfer to correlate the Nusselt number, Nu, against the product of the Prandtl number, Pr, and the Grashof number, Or. The Grashof number is

-78also considered a must in any work involving natural convection. However, these more customary dimensionless groups can be obtained by combining Eu, Re and Pe and the boundary conditions. For instance: GrL = 3 (38) V2 Let L, being a distance be converted to D. Also for a perfect gas B = 1/T, K1, so Gr= (Dg) T = Re2 AT (39) 2 T T Once the shape of the temperature boundary condition mentioned above is specified, the dimensionless group LT/Tref would be a measure of the magnitude of the boundary condition. The construction of the experimental apparatus described in Part I fixes the shape of the temperature boundary condition, thus tT/Tref has meaning in the case under consideration. Also the Prandtl number can be obtained as follows: Pr e = a = v = k.L L (40) Re MJ~g/v = P P k k Notice now that the usual correlating parameter, (Pr)(Gr), equals (Pe) (Re)(Tr/T). This is the simplest combination of the basic dimensionless groups in equations (34) to (37) that describe natural convection heat transfer. The Nusselt number is a direct result of the temperature distribution inside the sis tern. The local Nusselt number at angle Q is, DNu~ - =O O; _(aD kw = ( DT (141)

-79where the subscript w refers to the surface. Thus the common formula for correlating natural convection heat transfer data, Nu = f [ (Pr)(Gr) (42) may be in agreement with equations (34) to (37) because it relates a dimensionless group derived from solution of these equations to a simple combination of the chief dimensionless groups that fix the solution to these equations. Notice that the velocity boundary conditions and the Euler number, Eu, are not used to construct any of the canmmon parameters used to describe natural convection heat transfer. The reason for this fact is right in the words "heat transfer". Natural convection currents could probably be correlated employing Eu and some dimensionless measure of velocity along with other dimensionless groups. Simplification of Basic Equations In most instances of natural convection the pressure is very nearly the same as if the fluid were at rest. Also in many instances of natural convection only mild temperature differences are encountered so that the fluid density is practically a constant when compared with the fluid velocity. The following simplification of equations (34) to (37) is based upon the above two approximations. Suppose in Equation (28), p and q are everywhere zero, so: O = -g sin 9'gc bP (43) pr au r Physically, p and q would be zero only if the temperature were the same

-80. throughout or the temperature gradient pointed straight up throughout. In either case the temperature in a horizontal plane will be constant at To, and from equation (43) it is always true that -g sin GPo -.C PO o (44) pr 6( where the subscript o refers to the hydrostatic condition. Now let PD = P - PO and pD = P - Po and for u A O, -g( sini)p - = a D+ p)i (45) It is now assumed that Po aPD because 6PD/2O exists only because of gas movement which is of very low order. Neglect 6PD/89 and substitute 6Po/69 from equation (44) into equation (45)...g(sinQ) r - = p (PD + PO)g sins + gsing po = - PD g sin 6 So -g sin - c = po - g sin P (46) pr ae p and for essentially constant pressure, p is proportional to ~K'l, so.-g sin 8 j f Q =c [T O]g sin Q (47) Consequently, when this additional simplification is made in equations

(28) and (29) and equations (27), (28), and (29) are non-dimensionalized as shown before, the results are: + & g + =~ g sin L M-T1 + Re 2V (48) CT'os= g c T' v + Re'lvm ~ (49) Ad + + xPe" a =PT T (50) Finally, assume that density variations are very much less than velocity variations. So equation (30) non-dimensionalizes and simplifies to: i + = 0 (51) Now the equation system is reduced to four equations in three unknowns, T, p and q. Equations(48) and (49) express conservation of momentum, equation (50) expresses conservation of energy and equation (51) expresses conservation of mass. Since momentum, energy, and mass must all be conserved and since p is usually more significant than q, equations(48), (50) and (51) were employed in the program. CONVERSION TO DIFFERENCE EQUATIONS Before equations (48), (50) and (51) can be converted to difference equations, a grid must be set up. Also the manner in which the difference equations are used is fully as important as the difference equations themselves in describing the computation method. The above points will be described in the following section. The reader is advised to review the section entitled, "Formulation of Difference Equa

-82tion Solution Methods" in Part I (page 15) before proceeding. The Grid A plot of the grid used is shown in Figure 29. A cylindrical grid was used in order to fit the experimental geometry. Notice that the radial spacing is not uniform but it is crowded toward the outer edges where the rapidly rotating boundary layer was observed. From the standpoint of numerical analysis, the circular grid has two peculiarities not found in square grids. One is that velocity of gas at the central point cannot be expressed so it is fixed at zero. Physically, this restriction is equivalent to having a co-axial rod in place the size of the small circle at the center of Figure 29. The second is that the row of points at a given radius has no natural end. Consequently, two more columns (rays) of points must be added so that each circle shown in Figure 29 will have ends. Calculation starts in ray 2 with values in ray 1 fron the previous iteration as external points. Calculation ends in ray 1 with values in ray 2 from the previous iteration as external points. At the end of each iteration the values in the extra rays are brought up to date. This then gives 120 values for each variable that must be stored. It was found convenient in the planning of the program to write out the grid point nunbers as shown in Figure 30. Translation to Difference Equations In order to obtain a set of equations that can be solved explicitly, the forward divided difference in time and the central divided differences in space were used. Therefore, equations (48), (50) and (51) are written below in difference forn for point nunber 12 in Figures 29 and 30. Similar equations can be written for any other interior

rl-~~~~~~~~~~~~~~~~~~~~~~~~~~~~r 0 c0 0 I \ O!U!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ u~l!~~~~~~~~~~' U'!~~~~~~V!e ire..

Figure 30 GRID POIT NUMBERS Ray No. 1 2 3 4 5 6 7 8 9 10 1 2 1 1 11 21 31 41 51 61 71 81 91 101 111 2212 22. 2 42 62 72 82 92 102 112 ~~~7289..2 109 3 3 13 2 33 33 53 73 3 93 103 113 4 4 14 24 14 4 4 64 74 84 94 104 l l4 c 5 5 15 25 35 55 6 75 95 105 115 0 o 6 6 16 26 36 46 56 66 76 86 96 106 116 7 7 17 37 47 57 67 77 87 97 107 117 8 8 18 28 8 48 8 68 78 88 98 108 118 9 9 19 29 39 49 59 69 79 89 99 109 1 10 10 20 30 40 50 60 70 80 90 100 110 120 I —-~~~~~~~~~~~~~~~~~~0 I1 I120-II

-85point, T12- + 112 F T1 1 T11 P12 ___ - _ 13'T12 12 - _1 _T13 _Tl Pe-f1,...+V~L~3 22] + (52) _r3 r_ 2 r_2 +2 () P12-P12. -+ - +2 +q1~~2~2 1 -!P23' P12 P12 - Pl sin 2L + Rer r3 2 r2 r L To....... + (53) 1 4 1_ _ P-22 +P-2 1' 124 it4 - v 2 a e 3 o (54) The primed values above refer to qntities at a step ahead in time at. There are several ways that the iteration process can pro

-86ceed using the above equations. They are: (1) The total step, or Jacobian, process. (2) The single step, or Gauss-Siedel, process. (3) The half grid process. The total step process consists of computing T12', T22', T32, etc., and storing each in a separate register until their counterparts, the quantities T12, T22, T32, etc., are no longer needed in the computation before the old values are replaced by the new. The single step process consists of computing T12' and putting it in place of T12 and then computing T13' and putting it in place of T13, and so on until all internal points have been replaced by new values. The half grid process is a variation of the single step process in that only half the values assigned to the grid points are changed during each iteration cycle. For instance in Figure 30 the interior points that are underlined would be changed during the odd numbered iteration cycles and the other interior points would be changed during even numbered iteration cycles. The half grid process was chosen because no storage external to the grid storage was required and yet most divided differences employed were made up from values of the same generation. Consequently, it would seem that this method would have less error than the single step method where the divided differences are computed from varying proportions of old and new values depending upon what point is being calculated. The iteration process proceeds as follows: (1) Calculate, using equation (52), (T12) and all other (T')'s whose grid points are underlined in Figure 30.

-87(2) Find and print out |T'I T max (3) Revise the boundary conditions on T grid by: a. Take the average of points 12, 22, 32. 102, and put into points 11, 21, 31...101. b. Make points 1 to 10 the same as points 101 to 110 (see Figure 30). c. Make points 111 to 120 the same as points 11 to 20. (4) Calculate, using an equation similar to equation (52), (13) and all other (Ti')'s whose points are not underlined in Figure 30. (5) Find and print out IT' TI(x. (6) Repeat step three. (7) Calculate, using equation (53), (P'12) and all other (p')'s whose grid points are underlined in Figure 28. (8) Find and print out IP' - Pmax (9) Revise the boundary conditions on the p grid by: a. Points 11, 21, 31... 101 are always zero on both the p and q grid by reason of the grid used. b. Make points 1 to 10 the same as points 101 to 110. c. Make points 111 to 120 the same as points 11 to 20. (10) Calculate, using equation (54), (q')'s for the underlined grid points. Note that in rar 2 points 18, 16, 14 and 12 are calculated successively and in ray 3 points 23, 25, 27 and 29 are calculated successively and so on. (11) Repeat step 7 for non-underlined grid points. (12) Find and print out Ip' - Plmax (13) Repeat step 9. (14) Repeat step 10 for non-underlined grid points. This time

-88in ray 2 points 13, 15, 17 and 19 are calculated successively and in ray 3, points 28, 26, 24 and 22 are calculated successively. (15) Punch out T's and p's for grid points 22 to 29. (16) Repeat steps 1-15. The above iteration process coded out to around 950 commands on the IBM-650. The program, as optimized using S.O.A.P., required two and one half minutes of machine time for each iteration. During this time the computer generated a better estimate of the steady state value for each of 270 points in the grids. Many iterations must be performed before the calculation system begins to approach steady state. Directions for using the program are given in Appendix C. ANALYSIS OF STABILITY Since an explicit iteration process was employed in this mathematical investigation in an attempt to arrive at the steady state solution of the differential equations as the limit of transient solution, it is of utmost importance to perform a thorough analysis of stability before proceeding with the work of programming and computing using a machine. Much time and effort can be saved thereby not to mention many wasted hours of machine time. Stability analysis for calculation systems based on partial differential equations more complex than the heat equation (see equation (12) have not been presented in the literature. However, the reader is strongly advised to study the section entitled, "Analysis of Stability" in Part I before proceeding on.

-89Stability Analysis of Equation (52) to (54) using the von Neuman Method. - - Consider equation (53) which has been converted to the case of constant grid spacing (Ak= c1 and a- = c2), and to the point nomenclature of Figure 31. _p______-_J__- +. i+ljk Pi-lj,__ At r2j...2a. = sin i [TLi,-io] + Re'l Pi,+l,k+Pi,j-l,k-2Pi,j,k (55) 1 rPi+l, j-,k+Pi-1, j, k-2Pi, j, rj FIGURE 31 POLAR COORDINATE SPACE-TIME GRID NOMENCLATURE r. At i-ijk,fX'

-90Now, according to the method, let the numerical solution be represented by T + a, p + b, and q + c where a, b, and c are round off errors. Assume that the round off errors are independent of each other and have a distribution as shown in Figure 32. The numerical solution as well as the exact solution to the difference equation satisfies equation (55). The error equation then is: b -i'j'jk + li bij+,kb] i, k-l bi+l'j'k-bi'l'j,1 + C,jk Pi,j+l,k'Pi,jl,k 1 + bik i+l,j,kPi-1,, = 2hr j, 2 J (56) sin Qij,, + Rel bi,j+lsk+bi,j-l,k-2bi,j,k + bi,j+2,kk bij..lk L9 -J-1 k+1 J'j ( rj I I, I! 2. + 2....... 2ak 1 bi+l,,k+bi-lj,k-2bijk } Notice that in deriving equation (56), To is treated as a constant. However it actually is not a constant but is an average of all temperatures in a horizontal plane. This average would have much less error than any one of the grid points used to find it. Moreover, in a cylindrical grid, the grid points are not positioned in horizontal layers. So in the machine calculation values of To at a given elevation are determined by linear interpolation of the values of T on the vertical diameter of the grid. There are many types of errors involved in the use of this method of calculating T,. These errors do not prove to be critical in the analysis and will be ignored henceforth.

-91FIGURE 32 ASSUMED ROUND OFF ERROR FREQUENCY DISTRIBUTION Frequency'Magnitude of Round Off Error Let bijk = e e i e (57) Then ea(t+At) eipr i7 = bi,j,k eat Consequently equation (56) becomes: L t 3+ie' - Pi j Ieij L-e' + [ + ciJk I1P 1 F81 ]= sin i a.k1 bi,j,k L6r j 6 Toj, r-i~~ 1 ~~~(58) e- +e-(Z'P)2 + I e -e 1 eiy Y+ +e-) -2 x, ~7 + )

-92This reduces to e 1i+ ij k + PiJk i sin A n +.4i~j ~ k 8+ iJ,k ( i,j,k To ijk Re'l c ~ 2 + a (Rel ()2 2 + __ The quotients c/b and a/b in the above equation will now be evaluated using the frequency distribution assumed in Figure 32. Since b can be as close to zero as desired and since a and c can be at almost maximum magnitude, the quotients can be made as large as desired. However, the chances of such an event happening is vanishingly small. Assume that a bound for 90% of the quotients is satisfactory. If.le<lbl<e and la or cl<e then 90% of the quotients are considered and I c/bLnax. = a/b max. = 10. Now substitute the maximum magnitude for sin, cos, and the above mentioned quotients into equation (59) and apply the triangular inequality. I eOt-l < - 10 sin-Qi + Rei( + +1 2(60) + -1 + 10 + Pi~k + |ijk 1 r-j $ i,j,k | i,j,k i'g- ]

-93Now according to the von Neuman method, it is necessary and sufficient that Ieatl1 < 1. Thus e(,At-Il< 2. Consequently a bound for stability of equation (53) is 10 sin i 1 ( + 1 _4 \ To Re \(j)2 j(Lk) (j)2 j 5 rj ijk + 10 l8Fli k + + i (61) for all i, j, and k. Similarly, the stability criterion for equation (52) would be: Pe1i( 4 + 1 +_ 4 + lo 4f!(Z)2 rj(-r) (rjA)2 rj 59i,J,k (62) + 10 | + Pi+,k,j,k < 2 or i,j,k -rA r for all i, j, and k. Equations(52), (53), and (54) are the difference equations employed in the program. Equation (54), of course, is not time dependent and consequently requires no stability criteria. The stability criteria for equations (52) and (53) are given in equations(61) and (62) respectively for the case when T, p and q all change every iteration according to the iteration process which has been outlined. If it were desired to specify the temperature, T, at every point in the grid and compute velocity only, then only equations (53) and (54) would be employed and the error, a, would always be zero. The convergence criterion would then become,

-944 + 1'4 + (63) 8k Piijk 2 | Fr ij,k IrjL; A r i -k C If in addition it were desired to ignore q and make q = 0 at every point in the grid, then the convergence criterion would further simplify to: ~I~fI - & t r AO Pe-/ 4 +1 + 4 ) 1 I p e \(2 r (s)2 (r&)2/,8 i,J,k (64) + i, < 2 _ rAt On the other hand if it were desired to specify velocity, p and q, at every point in the grid and compute temperature only, then only equation (52) would be employed and the errors, b and c, would always be zero. The convergence criterion would then become: (4e' s 44) Pi rjIk I fi k 2 \(r) 2 ~(t~i) (~Ag)2/? If p = q = 0 and if the grid has X subdivisions in the radial direction and Y subdivisions in a tangential direction, then the criterion becomes 14 sin2 [(X-1)l/2X] + 1 4 sin2[(Y-)n/2Y1 2 (66) \ (JY)G2 rj(&s) (YrjA) / 2 At Notice that in no matter what form the above criteria are in, they say that if At is made small enough stability is assured. Of course,

-95sufficient digits must be carried in the machine so that at least a three or four digit correction is made at each point during the iteration cycle. Stability Analysis of Equations (52) to (54) using the Matrix Method4 By substituting T + a, p + b and q + c into equations(52) and (53) and by carrying through the same process illustrated under the von Neuman analysis, the following error equations are derived: Error equation for equation (52). aijk + 1 = aijk (1-Aijk~) + ai,j-l,kBi,j,k) + ai,j+l,k Ci,j,kAt + ai+l,J,k Di, j,kAt (67) + ai++,j,k Ei,Jkat + bi j,k Fi,jkt + Ci, j,k Gi,j, kt where Ai,j,k = 2Pe'1 (/(j+l r-j) (+ -j.1) / Qj k [ - +/()22 Bi, Jk =''_ + pe'l [(2/(rj rJ-j.1)(rj+l r ) rj+l - rjl.1 - r/rj (J+l r j-l)] CiJjk,,k + Pe ((2/(+l rj)(PJ+1 - + 1(i +1 - rj.1)2 rJ (rj+l rj-1) ]e Pe-1 -1 Di; -i,j,k -=Pe''~j~k'2' - rj2(AG)2'J' 2rj= rr(tQ)2

Fi, k =_- i+15Jlk' i-ljk ik) Gijk = - (, 9 ) Error equation for equation (53). bij k+l = bi j k (l-Hi,,kA) + T bkk Ii,k + b i,k + bi,k Ki,,k + bi+-l,j,k Kijk (68) where + Ci,jlk Mi,j,k + ai,j,k Ni,,k) Hijk:= 2Re'l + 1 1)(Fj+l- j- ) (Fj+l F)(j+(Pi+ljk Pi-ljk) +'j (A() 2rjAQ i J, k - 2 - _. -k... + ReO l - 2...-.?j....) Ji = _qiJk + Re-1 2 + I - r+-rjl &(j+l-rj) (Fj+ j-.1). ((rj+l. j 1) / Pi j k Re- PJk ReKiljlk -~P +....2 Li,j,k 2jtQ Aj2(Q)2 27rjLQ Yj (A)G Mji k ij+,Pi,j-l,k, ijk in rj+j- r-J- 1 In the same way as was done for the heat equation (equation (12) ) in Part I, the coefficients of the error equations (67) and (68) can be arranged in matrix form as shown in Figure 33.

i= 1,23 j 1,2,3 T Factors p Factors q Factors II 12 13 21 22 23 31:3233 11 12 13 21 22 23 31 32 33 11 12 13 21 22 23 31:32 33 2 C~~~~ 0 ~ ~~2 2 Cr D2.B R E F22 022 w 2 2 22 22 a 023 D2IJX 032,~- I'521 ~~1 1 x a6 I IL121 b12 0C II 0~~~~~~~~~~~~~~~~ I io w - _ &~zz424 L~a~~M~I C x b21123 I I w 2229 22 _22 2 IT23 b 0 232 k=I ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~l.. k~~~l "~~C' o0 C22E. k=l c22 k=2 Figure 33 Plan for Error Matrix-Error Vector Multiplication

-98But the general matrix is too complex to write out in detail in a reasonable sized figure. Consequently the notation had to be compressed and specialized for a 3 x 3 matrix in space. When T, p and q factors are considered as well as the T and p equations, the matrix grows to a flat 27 x 18. Even with this many elements, only two rows and two columns are as complete as they would be inside a larger grid. These are shown in Figure 33. The normal notation is transformed to the compressed notation as follows: Si-elS11 12 S13:Sij.l Si,J Sij+1 s21 S22 S23 Si+1,j S31 S32 S33 and 1-Ai,j,kZ = R22, 1-Hi,j,kFt = S22, Bijkt - B22, etc. Notice that if the matrix multiplication is performed for the row containing R22 with the error vector at k = 1, the error a22 at k = 2 is obtained. This process is the same as indicated in equation (67). Similarly, b22 can be obtained by multiplying the row containing S22 by the error vector at k = 1. (Compare equation (68).) The error c22 is evaluated from the algorithm employed to determine q's from p's. In this way all components of the error vector at k = 2 are computed. Then in turn the error vector at k = 2 is multiplied by the error matrix to get the error vector at k = 3 and so on. If the error matrix has all eigenvalues less than one, all components of the error vector should decrease. However, if one eigenvalue is greater than one, then the absolute value of the error vector will constantly increase and make the process unstable.

-99Up to this point the stability analysis by the matrix method is on a solid foundation. This method however can be improved by obtaining a closer bound to the magnitude of the largest eigenvalue of the matrix. Of course, actual computation of the eigenvalues of the matrix is possible but the work involved is much greater than the solution of the differential equation. One easily computed bound on the eigenvalues of a matrix was given by Gershgorin (sometimes attributed to Hadamard 27). This theorem states that the eigenvalues are in circles of center 1-Aijk Af and radius [IDI + BI + c + ElB + IFI + IGO]i i,k At or LDi+l,jkl + IBi,J+lkl + ICi,j-1.,k + Ei-l.,j,kl + liN,j,k Jt whichever is smaller. Or the eigenvalues are in a circle of center 1-E j,7k At and radius [III + I JI + KI + ILI + IMI + N]iJ,k t or [IIi,+ll + Iji,j-lkl + IKi+l,J,kI + I-,i-,,kI + IFi,j,k i whichever is smaller. For convenience, let Pi. jk At be the computed radius for the eigenvalue bound circle whose center is l-Aijk At and let Qi jk At be the radius for the circle with center 1-Hi jk At. Consequently, i'AiJk t' kijkI < Pijk iAt (69) and Il-Hi,j,k At - xi,J,k I Qi,J,k At (70) where X is the eigenvalue. ijk

-100For equation (69) the situation is seen clearly in Figure 34. In contrast to Figure 5, the circle bounding the eigenvalues and the circle bounding the area of stability are not necessarily tangent. In fact P > A. Figure 34 Matrix Method Stability Diagrams Imaginary Axis Circle Bounding Eigenvalues Real Axis _/'Circle Bounding Area of Stability Large At Small At Clearly if P>A, then no matter what the time step is, the circle bounding the eigenvalues of the error matrix lies partly outside the area of stability. It should be emphasized that the eigenvalues may still all lie within the area of stability. Thus if the matrix method indicates that the calculation process is stable then it is stable. If the matrix method indicates that the process may be unstable, then it still of course may be stable. Suppose for example that: Pi,j,k = [B + + I DI + I EI + I FI + IGI i l,J,k (71) If the radius happens to be smaller so much the better. Notice now that if I I 2Pe-1 ij, Pe-1 + (72)

-101and Pe-1 Pijk ~~~(ji~~~~~~.~~~1~~2L)2 bl F(73) Then, A= jBJ + ICI + IDI + IEI However, the coefficient F and G still make P>A and make the matrix analysis inconclusive. If it were desired to specify the velocity, p and q at every point in the grid and compute temperature only, then the error equation for equation (52) would not contain F and G because the errors b and c would be zero. Now since the velocities are fixed, it would be possible to arrange the spacing of points in the grid so that inequalities (72) and (73) would be satisfied. Then P = A and At < 1/A. This criterion again is the same as for heat conduction. That is, if inequalities (72) and (73) are satisfied and if p and q are fixed, the upper limit on the time step that will insure stability is the same as for no velocities at all. On the other hand, suppose for example that Q,jk 5[III + IJI + IK + IL+ IMI + INIijk (74) An examination of the definitions of H. I, J, K, L, M, and N shows that in order for Q<H and thus give a bound circle entirely within the area of stability, the first requirement is to have ap8/0 negative. Having a8/80 negative will increase H and allow for a possiblity of having a Q_5 H. If the temperature is fixed (i.e., a = 0) and if inequalities similar to (73) and (74) are observed, then for certain stability it must be that

-10211 _ > _P al a < 0 ( rj | bg | and < O (75) Of course, ap/a8 cannot be negative all the way around a closed circle. Also, from experimental observation Ig __a (76) Thus the matrix method states that the overall stability of equation (53) is in doubt for all real conditions. Comparison of Analysis of Stability by the Matrix Method and by the von Neuman Method There is a definite conflict between the two methods of analysis. The von Neuman method, as extended by the author, states that the solution of both equations (52) and (53) can be made stable if At is taken small enough. The matrix method says that the solution of equations (61) and (62) is not certain to converge for reasons related to the nonlinearity of the equations. On the positive side, the matrix method states that equation (61) will be stable if it is made linear by fixing p and q and if inequalities (72) and (73) are observed and if a small enough time step is employed. COMPUTATIONAL RESULTS After writing the program according to the outline given on pages 85 to 88 and after the program was "debugged" and thoroughly checked by hand, then it was ready to solve problems. The two purposes for doing the calculations were to:

-103(1) Peproduce the gas velocity and temperature distributions measured for experiment 4 by substituting the proper boundary condition into the program. (2) Investigate the stability of the calculation method employed. In performing these two purposes, the problems listed in Table VII were treated in three main categories: (1) Iteration on T and p together. (2) Iteration on p alone, (3) Iteration on T alone. Table VII also gives the boundary conditions employed for each problem so that comparison and reference may be made easier. Iteration on T and p Together This mode of solution uses the unmodified program outlined on pages 85 to 88. Data Figures 33 and 36 show the transient solutions for some grid points during the first run to produce useful results. The locations of the points referred to are given in Figure 29. These plots are possible because the temperature T and the tangential velocity p for point numbers 21 to 29 in the grid were punched out at the end of each iteration. Notice that here instability causes the value assigned to a point to increase or decrease exponentially.

TABLE VII. BOUNDARY CONDITIONS FOR SOME OF THE PROBLEMS ATTEMPTED Variable Internal Points Computed Initial Conditions Problem -l- 1 T _ _ No. T p q Pe Re- T T T p 16 x x x.0002293697.0001610922 1.061346790 0.938653209 0.5+ 1 0 0 17 x xx I it 0.2 1 0 18 xlxl " " i t 0.5, 0.2 1 0 0 20 x x x.0016109220 1 ".1, 0.2, 0.5 1 0 0 26 x I DNA* " 0.1, 0.2 1 Table XI 0 30 x x x ".0001610922 1.003333333 0.996666667 0.5 1 0" 31 x | xl | I "i I, " 1i 0.2, 0.5 Prob. 30 t = 28.5 0 H 34 x x x ".0016109220 1.061346790 0.938653209 0.01, 0.995 Prob. 20 t = 15.5 36 x.00100 DNA 0.9 0.5, 0.2 1 0 0 39 x.0002293697 DNA i 0.5 1 Table XI 40 x x x ".0016109220 I " 0.01, 0.001, 0.0001 Prob. 20 t = 15.5 41 x x DNA.0001610922 " " 0.25 Prob. 39 t = 53.5 0 42 x DNA "t " 0.1 0 43 x DNA It " 0.05 0 48 x.0002293697 DNA i It 0.90 1 Table XI * DNA, does not apply. + First iteration at At =.1.

-105I I. 1.00 Pt. 25 ~P t. 262 Pt. 27 Pt. 28.9 98 L I I9I7I l Pt. 29 TE~TIME t Figure 35 Transient Solution of Temperature Equation (Prob!.16)

-1060,15 0.10 Q05 I =. Pt 29 Pt. 27 0 5 10 TIME, T Figure 36 Transient Solution of Tangential Velocity Equation (Problem 16 )

-107This calculation has the same boundary conditions as does experiment #4. Problem 17 and 18 were performed at various time steps to check problem 16. All three transient solutions for the temperature of point 29 are compared on Figure 37. Notice that there is essentially no difference between them. It seemed reasonable that if a solution were converging, the magnitude of the maximum correction in the mesh point values should be decreasing. Figures 38 and 39 show how these maximum corrections behaved for the temperature equation and for the velocity equation respectively in problem 18. It was also observed that the values of velocities, p, on the vertical diameter was always an order of magnitude smaller than they should normally be. (See Table VIII.) TABLE VIII. SELECTED TANGENTIAL VELOCITIES p AT END OF PROBLEM 18 Distance from Wall.02 r.04 r.06 r Ray #1 -0.000465 -0.000595 -0.000536 Ray #2 -0.0389 -0.0347 -0.0300 Ray #3 -0.207 -0.182 -0.144 Ray #4 -0.0839 -0.152 -0.169 Ray #5 -0.0337 -0.0511 -0.0389 Ray /6 -0.000425 -0.000599 -0.000631 Ray #7 -0.0380 -0.0377 -0.0366 Ray #8 -0.215 -0.178 -0.137 Ray 19 -0.111 -0.199 -0.209 Ray #10 -0.0731 -0.0437 -0.0360

-108i.0c 1.99 A t 0.2(Prob. 17) Itat=0o.5 i- (Prob. 16).98 AT=0.5 To t=9.5 AT =0.2 after (Prob. 18).97 0 5 10 TIME t Figure 37 Comparison of Transient Solutions of the Temperature Equation for Different Time Steps (Pt. 29)

-1090.03 _E 0.02. I. I I 0O 5 10 TIME, t Figure 38 Maximum Correction in T During Iteration (Prob. 18)

-1100.06 I4- 0.04 Ia. 0.02 0I I,. I I I., I'0 5 10 TIME, t Figure 39 Maximum Correction in p During Iteration (Prob.18)

-111At this juncture it was reasoned that possibly the boundary conditions called for velocities higher than could be handled by the existing program. Thus if a fluid of higher viscosity were specified, the steady state gas velocities for the system would be slowed down and at the same time the velocity for the points on the vertical diameter would pick up faster since the only force acting to change the values of these points is viscous drag. Since Re-1 = v/-g D, Re-1 was increased by a factor of 10 in problem 20. A plot of the transient solution of the temperature equation for point number 29 only is shown in Figure 40. One would expect that this point would reach equilibrium early and afterwards would change very little. Strangely, the temperature gradually rises until eventually point number 29, which is right next to a cold wall, is hotter than average. Finally, the temperature curve breaks sharply and obviously becomes unstable. As shown in Figure 40, this break point was run through at At =.2 and At =.1 with very little difference in the results. Now the von Neuman type analysis extended by the author states that stability should exist at some time step provided sufficient digits are carried in the machine. Accordingly, several more attempts were made to reverse the trend toward instability starting at t = 15.6 in problem 20. As can be seen from Figure 41, there is no plotable difference between the solutions computed atAt =.1 to.005. The calculation was also carried out starting from the same initial conditions forAt =.001 and.0001. The computed value of the temperature at point 29 obtained using varying time steps is shown in Table IX.

(0o'qoxj6'6z ) uo*! b fenb l axn;lagadmaI jo uoT;nTos cua su'3u l, Oft aanOTa I I l I l I I I I II I lI I I I 1 1 1 i 1.00 - - 1.08 1.07 1.06.99 1.05 0J I-i — t:.A — =.. 1.04 c At i A2 ZT=.I A.-At=.2 —------ [ I I I I IiI I l I. I I 1 1 1 1 11= 1.02 1.0 I.97 1.00 0 5 10 15 20 TIM E, t

1.0019 n 1.0015 — Prob 20, + t= 0 r- _ Prob 40, o A t =0.01 o. Prob 34, A A-T =0 01 X ~os I 1.0010 15.5 15.6 15.7 15.8 TIME, t Figure 41 Transient Solution of Temperature Equation at Very Small Time Steps

-114TABLE IX. COMPARISON OF SOLUTION PROGRESS FOR SMALL TIME STEPS......______ Temperature of pt. 29 l t =.1 At =. 01 1 t=.001.0001 Prob. 20 Prob. 34 Prob. 40 Prob. 40 15.5000 1.000953591 1.000953591 1.000953591 1.000953591 15.5001 1.000953997 15.5002 1.000954403 15.5005 1.000955621 15.5010 1.000957659 1.000957651 15.5020 1.000961727 1.000961711 15.5050 1.000973934 15.5100 1.000994276 1.000994283 15.5200 1.001034990 1.001034995 15.5500 1.001157248 15.5600 1.001360450 1.001361209 Except for At =.1, these values are identical to nine significant figures. However, the maximum corrections during an iteration are not all decreasing even at the very smallest time step tried. For instance Table X shows that the maximum increase for one pass is increasing and the other is decreasing for both the temperature and velocity equations. Usually in such a situation eventually the corrections for one of the passes reverses its trend (see Figures 44 and 55) but the calculation was moving too slow in computor time to run this one out. Returning to a description of problem 20, Figure 42 is a plot of the buildup of a ray of tangential velocities. These curves begin to look peculiar when some of the velocities start decreasing. The break in the temperature and in the velocity curves occur at the same time. As in problems 16 to 18, the maximum corrections during an iteration are instructive. The maximum corrections in T is plotted in

-115TABLE X. MAXIMUM CORRECTION PER ITERATION FOR SMALL TIME STEPS For LA =.001 (Problem 40) lT - T Pmax/L|t - P /mxt/ L l1st Pass 2nd Pass 1st Pass 2nd Pass 15.501.00oo4o68.004324.005256.003816 15.502.00oo4068.00434.005246.003819 15.503.00oo4069.004324.005236.003822 15.504.004069.004324.005227.003824 15.505.004069.004323.005217.003827 15.506.oo4069.004323.005208.003830 15.507.004070.004323.005199.003832 t!!! t I t t I 15.520.004072.004321.0050o86.003862 For Ct =.0001 (Problem 40) llIT- I,/ | j|p Pimax/A l t 1st Pass 2nd Pass 1st Pass 2nd Pass 15.5001.00406.00432.00525.00380 15.01.00406.00432.00524.00381 15.02.00406.00432.00523.00381 15.029.004o6.00432.00522.00381

.015 A~~~~V-t==.l la. AI=.I AI=.2 / >.. -- - I.010 I.. - / -d r //A w Ozer~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ q~~~~~~~~ 0 - 0 5 10 15 20 TIME, t Figure 42 Progress in Solution of Velocity Equation (Prob. 20)

-117Figure 43. Notice that most of the time the maximum error is constantly decreasing. There seems to be no reason for the sudden jump at t = 4.2. Notice also how rapidly the maximum temperature correction increases when At =.5 after an induction period of two iterations. Also after t = 12.5, the maximum correction gradually increases until the break. Even after the break there appears to be an intermittent tendency for the maximum temperature correction to decrease. Figure 44 is the counterpart to Figure 43 for the tangential velocity p. Notice here that no evidence of a trend toward decreasing velocity correction is evident. The first jump at t = 4.2 is present as before. The second jump could be blamed on the change in time step even though this change did not effect the maximum temperature correction. The gradual rise sets in sooner, t = 10.5. Notice that the maximum correction during the second half of the iteration lags behind the first half until t = 15.9 at which time the corrections begin to increase exponentially and drive the temperature equation into unquestioned instability at t = 16.9. Problem 20 was an attempt to decrease steady state velocity by increasing viscosity arbitrarily. Another attempt was made in problem 30 to decrease the computed velocity by changing the temperature boundary conditions to comply with experiment 8 which has a wall temperature difference of only 20C instead of 40~C as in experiment 4. The progress of solution as well as the maximum corrections for both T and p seem to behave similar to problem 20 except that no sharp break and exponential runaway occurs. Nevertheless, the problem stopped because of scaling difficulties. A plot of the velocity vectors for all grid points shown in Figure 45 shows immediately that some velocities in ray 4 have become

.05 D04 K.02.01 Ct=-0.5 0 0 5 10 15 20 TIME, t Figure 43 Maximum Correction in T During Iteration (Prob. 20)

0.008 AT =.5 Ist Iteration Ist Hblf 8 0.006 Iteration 8 S~econd half Iteration 0.004 10. 0.002 80)0000 000 00 I0 I 0 I IIAt =.I I II=.2 0 5 10 15 20 TIME, t Figure 44 Maximum Correction in p During Iteration (Prob.20)

of caalqoad Jo pusa as sa9GI.TOOGA Jo uoT.;'eruasadias ao;DGA 5t alnsls. 01 5~~~~~~~~~~~ 6 8 I F 9 -0oz71

-121extremely large and positive; all other velocities are fairly reasonable. This same thing eventually happened to the velocities of points 27, 28 and 29 of ray 3 in problem 20 (see Figure 42). Figure 46 shows that an area of abnormally high temperature also exists in the region of instability. Notice in Figure 45 that the radial velocity component, q, causes the velocity vectors to appear more erratic than they actually are if the tangential velocities alone are considered. This is true because of the algorithm that was adopted for calculating q. On any one ray, alternate q's were calculated starting from the wall. The remaining q's were calculated starting from the center. If a sufficiently large number of points were used, the errors manifest in this algorithm would not be so apparent. Because of the already circular character of the flow, the effect of the radial velocity, q, should be second order. Accordingly, problem 31 was begun at the conditions existing at t = 28.5 in problem 30 and was carried forward so that q = 0 at all times and places. Nothing untoward happened to any point in the third row. Nevertheless, the solution stopped at t = 27.7 because of scaling difficulties. The final print out of the values at all the grid points shows that points 35, 36 and 37 in the p grid had values 10 times normal and were positive. Points 45, 46 and 47 had also grown to the same magnitude but were negative as is normal for the p's. Here again the consequences of instability are quite localized at least at first. A glance at the maximum correction during an iteration plotted in Figure 47 shows that the maxi-.mum correction in both p and T decreased for a time and then increased rapidly. An attempt to cause the maximum correction to decrease by decreasing the time step was in vain.

-1221.000 0996 66-'01 /~~~. ~ 0/ /.01 r/ 8/ v/1.00 0.9666 1.O333 ~~~~~~~~ Figure,1 46 I0r9o5 0~~~~~~~~~~103 0.96r 10-1 1.0333 ~1. 000

0 o. 0o oo.08 -0 t=. I and.01 for A, ~o r 0001 / 8 __ [ 8 8 1-_~~~.k /8 ~'~~IT -T Imx/,, /A 0.0 0 0 I I I I I i I I I 0 10 20 TIME, t Figure 47 Maximum Correction During Iteration (Prob. 31)

-124Thus, as is seen from the above account, the data taken with the program which was set to iterate on both p and T all eventually became unstable. Higher viscosity, a lower AT and smaller At's were tried but to no avail. Discussion Since no problem in the above account came near settling down to a steady state solution, the chief job of this discussion is to explain this difficulty. The von Neuman type stability analysis made by the author gives values of At well in line with those actually tried. For the initial conditions of problem 16, equation (61) was evaluated at point 26 as follows: 2 9.51 + 0.000161 (10,000 + 119 + 230) < -At At <.179 Point 26 was used because it has the largest sin Q, 0.951, the smallest radial grid spacing, 0.02, and the smallest radius for the small grid spacing, 0.42. Also for the initial condition of problem 16, equation (62) states the At < 0.845. Notice that using equation (61) the increased viscosity in problem 20 leads to a smaller time step, At = 0.0765. Starting at t = 15.6 in problem 20, several attempts were described in the data section to induce stability by decreasing the time step in accordance with the von Neuman type analysis. Time steps as small as 0.0001 were used with inconclusive results. At the above mentioned starting point the following maximum magnitudes were observed: ST = -1.445 at pt. 29 5?

-125+0.0992 at pt. 59 p = -0.05025 at pt. 95 q = -0.0181 at pt. 97 — P = -1.08 at pt. 94 8r B = +0.0732 at pt. 5 Strictly, equations (61) and (62) should be evaluated at every point in the grid and the smallest time step obtained should be used. However, if the maximum value for each variable were put into one equation, a conservative result would be obtained. If this substitution is made, equation (61) gives At <.0855 and equation (62) gives At < 0.101. Clearly then the von Neuman analysis as extended by the author is useless in this application. The Matrix method, as already stated, gives inconclusive results for all cases in this type of calculation. Nevertheless the concept of an error matrix is useful in giving a possible explanation for the behavior of the preceding problems. Evidently certain eigenvalues of the error matrix governing the error growth in one part of the grid are always greater than unity no matter what the time step. They can be made as close to unity as desired by decreasing the time step but this tactic only forestalls the inevitable. Notice in Figure 37 that the use of At =.2 instead of.5 increased the amount of real time that the transient solution remained reasonably stable by only 10 percent. For this slight advantage machine time was increased from 50 to 125 minutes. Notice also that when the real runaway occurs the smaller time steps produced a faster runaway in real time.

-126As is explained in the logical description of the problem on pages 85 to 88, the maximum correction for T and for p were accumulated for each pass through the grid. Two passes were required to change all the T's and two more for the p's. One particular pass always changed the same set of points. The maximum correction for both passes were plotted and connected with a vertical line to show that they belong together. The maximum correction for a given pass during successive iterations does not necessarily refer to the same grid point. Consequently, Figure 39 is very erratic while Figures 38, 43, 44 and 47 have some smooth parts. The real cause for the sudden jumps in the maximum corrections was not determined. A cathode ray read out device connected to a computer would almost be a necessity to hunt for the cause of such an effect. The maximum correction feature was originally incorporated into the program to furnish a criteria for stopping the program after a satisfactory approach has been made to the steady state solution as determined by equation (25). Notice that dividing the maximum correction by the time step employed satisfactorily normalizes the results. Convergence to a steady state solution is possible only when the maximum correction constantly decreases. In the problems so far described, the maximum correction usually decreased or remained constant for some time and then increased rapidly. (See particularly Figure 45.) Now if the eigenvalues of the error matrix are close to but sometimes greater than unity, the true transient solution to the problem will be followed for a time while the error is building up in the non-significant portion of the ten digit number carried in the machine. Once the error

-127becomes significant then the transient solution appears to blow up suddenly when really it has been blowing up from the very first. Iteration on p Alone A study of both the von Neuman and the matrix method of stability analysis shows that anything that can be done to make the equations independent will help the stability of the calculation. If the computations of p and T are stable independently, then a solution to the general problem can be worked out. Of course, if the first iteration for a given set of constant temperatures or velocities is convergent, then alternation of the iterations as is done in the original program would be a convergent or stable process. Ample evidence has already been presented to conclude that the original program is not stable under round off error. The hope is in performing the following calculations that the iterations on velocity will become and remain stable after a number of iterations performed holding the temperatures constant. Problem 41 was done with the set of T's from the end of problem 39. The T's remained fixed and iterations were made on p with q being computed. The time step was 0.25. Problems 42 and 43 were duplicates of problem 41 except that the time step was.1 and.05 respectively and the q's were kept at zero. All three of the above problems ran into difficulty because of scaling and stopped. It was found by inspecting the final print out of all values, which was taken after the program stopped, that in all cases point 29 was the cause of stoppage. Point 28 was next largest (0.156). Figure 48 shows how point 29 progressed toward overflow. Notice the lack of any exponential increase.

-128Overf low Prob 41 0.20 0.20 A.25 / Overflow ~/ ~ AProb 42 AT=.5 g hi Prob43, At=.05 / / same as prob 42 within accuracy of graph a0.15 / I TIME, t Figre 48 Transient Solution of Velocity Equation with T Constant (Probs. 41, 42 and 43)

-129The plot of the maximum correction of p during an iteration cycle shown in Figure 49 exhibits the same type of behavior as was experienced many times before. (See Figures 38, 43, 44 and 47. ) Discussion The von Neuman analysis for stability analysis as extended by the author gives equation (63) as stability criteria for problem 41 and (64) for problem 42. Since equation (64) will give a smaller time step, it will be evaluated at the initial condition for problem 42 or 43 to see what time step would be indicated. As is shown in Table VII, the velocities in problem 39 remained constant. They are tabulated in Table XI and are graphed in Figure 50. The maximum magnitude of ap/6Q is at point 49 with a value of 0.0882. The maximum magnitude of p is 0.090. If the same conservative treatment is adopted as was explained in the previous discussion, equation (64) evaluates as follows: Re-1 [4 + 1 + 4- 41 p + I I <2 -1 4 r.(/hr) (r.A~ r,~ j - i,j,k rALQ- At (64) 0.0882 o.090 2 0.000161 [ 10,000 + 119 + 230] + 0 0882 +.914)< 2 <t 0.78 This time step is of course much larger than those that did not work (At = 0.1 and 0.05). Notice that equation (64) does not contain the factor 10 that was derived by considering only 90 percent of the possible error quotients. Equation (64) considers all the errors and can be derived by a straightforward extension of the example given in Part I.

Legend o Prob 41, A =.25 11 Prob 42,At =O10 o Prob 42,At =.05 0.06 <3. 0.05 3 E 0.04'n A 0[ 0.03 0.02 0.01 - 0 I 2 3 TIME, t Figure 49 Maximum Correction in p During Iteration (Problems 41, 42)

-131TABLE XI. ESTIMATE OF EXPERIMENTAL VELOCITIES FOR EXPERIMENT NUMBER 4 q = 0 for all internal points. p = 0 for all internal points except: Location Location Point No. on Drum Value Point No. on Drum Value 5 1125 -0.055 65 1185 -0.0275 6 1126 -0.065 66 1186 -0.050 7 1127 -0.065 67 1187 -0.0575 8 1128 -0o.o055 68 1188 -0o.o055 9 1129 -0.040 69 1189 -0.0475 16 1136 -0.040 76 1196 -0.045 17 1137 -0.065 77 1197 -o.o65 18 1138 -0.0725 78 1198 -0.080 19 1139 -0.055 79 1199 -0.080 26 1146 -0.025 86 1206 -0.045 27 1147 -o.o65 87 1207 -0.065 28 1148 -0.090 88 1208 -0.080 29 1149 -0.070 89 1209 -0.080 36 1156 -0.025 95 1215 -0.0275 37 1157 -0.065 96 1216 -0.055 38 1158 -0.090 97 1217 -0.065 39 1159 -0.070 98 1218 -0.0675 46 1166 -0.040 99 1219 -o.o60 47 1167 -0.0725 105 1225 -0o.055 48 1168 -0.050 106 1226 -0.065 49 1169 -0.0075 107 1227 -o.o65 55 1175 -0o.55 108 1228 -0.055 56 1176 -0.055 0og9 1229 -0.040 57 1177 -0.050 116 1236 -0.040 58 1178 -0.030 117 1237 -. 065 59 1179 -0.015 118 1238 -0.0725 119 1239 -0.055

-132L,.1.,,,1,I SCALE FOR P 0.05.10 f COLD SIDE HOT SIDE Figure 50 Assumed Velocity Distribution (Probs. 26, 30, 39, and 48)

-133Evidently the key to the trouble is the non-linearity of the velocity equation (see equation (48)). The term (p/r) (8p/84) is a necessary part of the velocity equation. If two sets of p's were carried in the machine and if one set remained constant and were used to evaluate p/r above and if the other set were being revised by iteration and were used to evaluate a/69 and the other partial derivatives in their divided difference approximation, then the equation would be linear. This idea was tested out in principle in the next section. Iteration on T Alone If all the velocities of the system are specified and are held constant then the temperature equation becomes linear. Practically speaking, this case is of the least value because few, if any, situations exist where velocity data can be measured and temperature data cannot. Nevertheless, some encouraging results were obtained with this variation of the basic program. Data The simplest problem done with this program variation was problem 36 in which all velocities were fixed at zero. This problem is nothing more than the solution of the heat equation. (See equation (12)). It was found that wild fluctuations in temperature occurred when At =.9 and At =.5. However, the iteration procedure settled down quickly when At =.2 (see Figure 51). Strangely, however, the maximum temperature correction did not appear to be so wild. According to Figure 52, this correction varied only from.05 to.1. But physically the range for T is.94 < T < 1.06 so.it was relatively very large. ifter the time step was changed to 0.2, the maximum temperature correction

ATL.9 =-2 Change in scales 1.00~~~~~~~~~~- 1.00 1.00 -Iw a:.959 Ui Ld I-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 90- I.F 510 40 60.. TIMEt Transient Solution of Heat Equation (Prob. 6) F90C~~~~~ \I +~ — "' ——'+-"'-t~9 0 ~~~~5 I 0 14 20 40 60 75 TIME, t Figure 51 Transient Solution of Heat Equation (Prob. 36)

-13510o-~~21 ~ ~~~ LV~: \ Maximum ess 8 89 10 I CT 0 20 40 60 TIME t Figure 52 Maximum Correction in T During Iteration (Prob. 36)

-136dropped extremely rapidly and eventually became a case of exponential decay. A plot of the isotherms is given in Figure 53 in order to show what kind of solution was obtained at t = 75 (over four hours' machine time at 45 seconds per iteration). According to equation (25), the error bound on the steady state value may be 81 times the maximum temperature correction. Consequently, according to Figure 52, the values for the grid points may possibly be taken to be the steady state values with an error no greater than + 0.002. That is, the isotherms may be no farther off than one fifth the distance between the isotherms. Thus the biggest error in this calculation is now truncation error since too few radii were employed to give a good approximation of the temperature boundary condition (see Figure 21). Next a set of velocities were assumed based upon the measurements taken in experiment 4. (See Table XI and Figure 50.) This set was only a first approximation of the true values in that only the tangential velocities in the annular boundary layer are included and no effort was made to make the set of velocities conform to equation (51) which is an approximation of the continuity equation. Starting with the same initial conditions, problems, 26, 39 and 48 were carried forward with time steps of 0.2 and less, 0.5, and 0.9. The beginning of all three solutions were similar as is shown in Figure 54. However, as the time step increased above At =.2, the transient solution became less exact. However the larger the time step the sooner the approach to steady state would be accomplished. Problem 26 was only carried out to t = 25.3 and required about 2-1/4 hours of machine time. Figure 54 shows the transient solution for the points along ray 3. Figure 55 shows the maximum correction in T

-137T.99 1.001.01.05.93 9 06 1.061.939 1.061.93=J I.o 61.939, 1.061 Figure 53 Isotherms for Pure Conduction (Prob. 36)

-T- ~ Pt. 23 Pt.2 Pt. 25 Pt. 26 L~ rl I Pt. 27 nLIJ ~A T =.9 (Prob 48) lPt. 28 IEf I —- iA t =.5 (Prob 39) 3 Pt. 29..29 ( Pr ob.26) t. - -2 zT:.. 2 -2 AT.98 0 10 20 TIME) t Figure 54 Transient Solution of Temperature Equation (Probs. 26, 39,and 48).

0.06 QOI _0 N ~~0 0 0 E 8 I 0 0 H 0.00 10 0 10 20 TIME) t Figure 55 Maximum Correction in T During Iteration (Prob. 26)

-140during an iteration considering all points. Notice that the two jumps in Figure 55 correspond in time to the points of inflection of the transient solution of point 29. (See Figure 54.) Problem 39 was carried out to t = 125 and required approximately 3-1/4 hours of machine time. Figure 56 shows the transient solution for the points along ray 3. Notice that after t = 30 or 40 the closely spaced points near the edge of the grid (i.e., points 25, 26, 27, 28, and 29) have practically reached equilibrium and are only affected slightly by the subsequent large change in the values of the widely spaced points in the core (i.e., points 22, 23, and 24). In fact at t = 125 the temperature at point 21 shows no signs of coming to equilibrium. Nevertheless the temperatures at t = 125 were plotted and from these plots the isotherms shown in Figure 57 were constructed. The maximum correction made during an iteration in problem 39 is plotted in Figure 58. After some initial jumping around, the maximum correction drops rapidly and eventually settles down to simple exponential decay. The grid used in this problem (see Figure 29) is 18 grid spacings in diameter. Thus by equation (25) the error bound on the steady state value may be 81 times the maximum temperature correction as shown in Figure 58. The temperatures in problem 39 varies from T = 1.06 to 0.94 approximately. As can be seen from Figure 58, the maximum error involved in taking the transient solution at t = 125 as the steady state solution may be 0.01, still quite a large value. Problem 48 was carried out to t = 401 and required approximately 6 hours of machine time. Figure 59 shows the transient solution for the points along ray 3 for the initial part of the solution. Figure 60 shows the transient solution for the same points plotted to a much larger

II-~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 1.00 Pt. 22 Pt. 23 ~ 99 ~ Pt. 24 w I~-.98 Pt. 25 <~.96 Pt.29.95.94 50 I00 TIME t Figure 56 Transient Solution of Temperature Equation (Problem 39)

-142T.99.$ LOOI'02 3 04.19 J.05,.01.9 VF.~Z~' 06.9.06.06 /igure 5 Calcul~ted Istherms (rob, jg t.1if l~~~~~~~~~~~~~~~

Maximum E s 0o0 I 0 8 0~I 0 C I~' ITT I max/AT 0.001 O I 00001 50 100 TIME, t Figure 58 Maximum Correction in T During Iteration (Prob. 39)

-144Pt 23 1.00 Pt 24P Pt 25.99 Pt 26 9.'.98 Pt 27 w.97 Pt 28.96 Pt 29.95.94 I 0 10 20 TIME t Figure 59 Transient Solution of Temperature Equation (Prob. 48 Initial Solution )

1.00 Pi. 22.99Pt % Pt. 25~~~~~~~~~t.2 Pt. 24~ \~.98r Pt. 26 Pt. 27 Pt. 2.8 W-.97 Ui Pt. 29.95 2-00- 300 0 too -Tit-AE 0 Frb 4B Covml-ete 01 Soluion"f Tmerature Equation ('b Figure 60 Transient ltonfTea

-146time scale. Notice that even at t = 401, points 22, 23 and 24 do not appear to have come to equilibrium. The maximum correction made during an iteration is plotted as a function of time in Figure 61. As in Figure 58, the maximum correction goes into exponential decay. However, the rate of decay decreases slightly at t = 200. According to equation (25), the maximum error involved in taking the transient solution at t = 401 as the steady state value may be 0.0013, a reasonable error. Consequently the values of all grid points at the end of problem 48 were cross plotted to obtain the isotherms shown in Figure 62. In order to visualize the temperature surface that was computed, a three dimensional model of it was made. (See Figure 63.) Discussion Finally, a calculation was carried to the point where comparison is possible with the experimental measurements. So in Figure 64 the experimental dimensionless isotherms for experiment 4 were constructed. There is a close resemblance between this plot and the computed isotherms of Figure 62. The chief difference apparent from these two figures is that the temperature gradient at the center of the circle is inclined 200 from the vertical toward the hot side in the computed results while the temperature gradient is exactly vertical in the measured results. The temperature gradient at the wall is not apparent from a look at Figures 62, 63 and 64. Nevertheless, it determines the amount of heat the space will transmit by convection. Figure 65 shows a comparison of the experimental data from experiment 4 with the end result of problem 48. If the thermal conductivity of the gas is not considered a function of temperature then the integral of the curves in Figure 65

10-2 Lo _~~~~~on~'-~ ~~MAXIMUM Ess 10-3 IO-5L 00 200 300 400 TIME, t Figure 61 Maximum Correction in T During Iteration (Prob. 48)

TO~ - 4 (9~'qoad) smaaq;osI p@~InnoTSO E,9 Garu7, 90'1 901' 90, 90'1,~l' SOII~~~~~~~~~ 90 t?6'j~gi~66 9 k~~~~~ / 6 \ ~~Io —-/

. -.-. -. - -ZZ - - —... — 7- - - -z- —- --- - -- H Figure 63 Three Dimensional Representation of Temperature Surface (Problem 48 t 401)

-150are _ = TH =1.06.99.94 Not e: Figure 64 Experimental Dimensionless Isotherms (Exp. 4)

-151Legend o Exp. 4:a a Prob. 48 20 cn 0 + or) -3 0 120 240 360 AZIMUTH ANGLE (0) L~ocal Dimensionless Temperature Gradients

-152from 0 to 360~ should be zero. Physically this means that no heat is lost. That is, all heat leaving the hot surface goes to the cold surface. The experimental graph very nearly integrates to zero. However, the calculated results show much more heat input than removal. The probable reasons why the calculated results agreed no better with the theoretical measurements may be as follows: (1) The calculation was not carried out far enough. Equation 25 is not true. For instance in problem 48 at T = 80 the temperature T at point 22 is 0.9960 and the maximum steady state error according to equation (25) is 0.0135. However, at t = 401 the temperature at point 22 was O.9812, a change of 0.0148. Furthermore this temperature is still changing with time at a considerable rate. (See Figure 60.) (2) The grid contains too few rays. Twenty or thirty rays instead of ten would probably be sufficient to reduce the truncation error by allowing the temperature boundary conditions to be specified more like the experimental system. The central divided difference approximation of the first derivative has the advantage of symmetry over either the forward or backward divided difference. However, when too few points are used to describe a function, the central divided difference causes the results to be extremely irregular, that is, alternate points appear to go together. (See Figure 65.) (3) The specified velocity field was not good enough. From the existing information more pains could be taken to get a set of velocities that would satisfy equation (51). However the real solution to this difficulty is to find a solution to the velocity equation so that the velocity field can be computed along with the temperature field. Only when the two can be computed will the calculation method have any practical value.

-153Comparison of time steps predicted from the theory presented herein with the time steps actually used in problems 26, 39 and 48 has led to some unexpected conclusions. Equation (65) is the relationship that should govern this case of specified p and q and variable T if the von Neuman theory is valid. For point 28 which has the highest velocity, equation (65) evaluates as follows: -1 _4 1 4 2 Pe + + < - (65) Ar)2 r (r.AG) rXQy jt 3J J 3.0002294 (.46 x.314)2).46(.314) <2 At < 0.669 Nevertheless a time step of 0.9 was used without experiencing instability. The matrix method of analysis states that this variation of the basic program is stable if inequalities (72) and (73) are satisfied and if the time step is as small as it would have to be for pure conduction alone. For the same point 28 above, inequality (73) evaluates as follows: Pe -1 1 Pijk > k (73) (r jA) 2r AG 0.0002294 +0.090.46(.314) ~ 2.00159 j o.045 The inequality is not satisfied so the matrix method still gives an inconclusive answer.

-154However on the basis of the limited experience described in this thesis the matrix method can be modified slightly so that it will more accurately separate the convergent from the non-convergent cases. The modification is simply to remove the absolute value signs from the expressions that are used to find the radii of the eigenvalue bound circles. (See page 99.) This modification would be in harmony with the two main conclusions of this mathematical investigation. Exceptions to these conclusions were not found but probably do exist. These conclusions are: (1) Non-linear differential equations are always unstable when solved numerically by an explicit iteration process. (2) If differential equations of the form of equations (48), (49) and (50) are made linear, they can be solved numerically by an explicit iteration process, and the largest stable time step is found by the same criterion as for the simple heat equation (12). It should be emphasized that in accordance with conclusion 2, the magnitude of p and q have no effect upon stability. The von Neuman analysis for coarse grids outlined in Part I gives a At that is slightly larger than any of the experimentally successful time steps. In the grid shown in Figure 29 there are only five small intervals in the radial direction and ten equal intervals in the tangential direction. pe-1 4 sin2[(Y-l)/2Y] 1 4 sin2 (X-1) t/2Xf\ 2 Pe + < (66) rj() (r )2( (E~) i~(AY)(VA) -- At.0002294 4 sin 72~ + 1 + 4sin 81~ (02)..42(.02) (.42x. 31416)2 At

-155-.0002294 (9045.2 + 119.0 + 224.1) > 2 At At < 0.9287 For an infinite grid with the same grid spacing, Az < 0.8424. The matrix method does not as yet have any provision for small grids. The criteria is: At 1 = 1/2 Pe 1 + 1 n -A k(n) (K)2) (77) or if substituted in as above, At = 0.8519 Thus it is seen that no matter what criteria of stability for the heat equation is used, the time step is about the same. Problem 36 in this section was computed using a Pe-1 of 0.001 instead of 0.0002294. Equation (66) gives a At of 0.213 which is in harmony with the experimentally successful At of 0.20. Conclusions The results of the computations performed with a program devised to solve equations (48), (49) and (51) by an explicit numerical procedure on the IBM-650 lead to the following important conclusions: (1) Only linear equations can be solved by this method. (2) The maximum stable time step is the same as for the heat equation (equation (12)). (3) Equation (25.) is not correct. More investigation is necessary before the true relationship is apparent.

CO11NCLUSION Experimen ta For the case of natural convection heat transfer inside of a horizontal cylinder of essentially infinite length and with the right half wall at a constant temperature higher than the constant temperature of the left half wall, it was found that an almost stagnant central core of gas exists. A boundary layer travels up the hot side and down the cold side and accounts for the convection heat transfer. Small changes in the average velocity of the boundary layer cause flow into or out of the boundary layer and cause a rather comPlex flow pattern in the central core. Since the boundary layer as a whole is hot on top and cold on bottom, the stagnant gas core acquires a vertical temperature gradient and thus is able to resist the drag of the boundary layer which otherwise would cause the core to rotate as a solid wheel. Quantitative analysis of the phenomenon described above revealed that: (1) The rapidly rotating boundary layer rotated less and less rapidly as the wall temperature difference increased above 200C. Although the scatter of data is large this decrease appears to effect the local velocities at fixed pointsin space as well as the overall rate of gas circulation. (2) The average Nusselt number based upon the cylindrical diameter, the wall temperature difference and the average temperature gradient over either the hot or the cold surface was found to be constant at 7. (More precise measurements will probably show a mild variation of iu with AT) -156

-157(3) The local Nusselt number was observed to be a definite function of the wall temperature difference at some points and essentially not a function of wall temperature difference at others. Computational Five partial differential equations in five unknowns were developed which involve the conservation of heat, mass and momentum as well as the equation of state. These equations were reduced to three in number by assuming a perfect gas and essentially constant pressure. These three equations in three unknowns were then converted to explicit difference equations for a cylindrical grid. This set of difference equations was solved by starting with the desired boundary conditions and with an artaitrary set of initial conditions. Then the approach to the steady state solution was accomplished by iterating through a transient solution to the equation. It was found by test that, (1) Only linear equations can be solved by the above method. (2) The maximum stable time step is the same as for the heat equation (Equation(12)). (3) Equation(25)is not correct. More investigation is necessary before the true relationship is apparent. The matrix method of error analysis is able to explain the above limitations in the solution. Given the conditions at any time, an error matrix can be computed which when multiplied by a vector made up of the errors at all the grid points at time t will give the error vector at time t v At. When the velocity and the velocity gradient become too great, somne of the eigenvalues of the error matrix. become

-158greater than one no matter what time step is used. This condition causes the error to increase constantly in some parts of the grid. This local error will eventually overwhelm the rest of the solution.

REFERENCES 1. Batchelor, G. K., "Heat Transfer by Free Convection Across a Closed Cavity Between Vertical Boundaries at Different Temperatures," Quarterly of Applied Mathematics, Vol. 12, No. 3, pp. 209-233. 2. Beutler, J. A., "Programming of Kinetic Calculations for Automatic Computation," Chem. Eng. Prog., 50, #11, (November, 1954). 3. Busey, H. M., Hammond, R. P., "Recombination of Radiolytic Gas for Aqueous Nuclear Reactors," Chem. Eng., Prog. Symposium Series No. 13, Nuclear Engineering Part III, p. 27, American Institute of Chemical Engineers, 25 West 45th Street, New York 36, New York (1954). 4. Carr, J. W., III, Class Notes on Numerical Analysis, University of Michigan (Sp. 1956). 5. Collatz, L., "Numerechebehandlung von Differential Gleichungen," Springer-Verlag/Berlin, Gottingen, Heidelberg (1951). 6. Courant, R., Friedrichs, K., and Lewy, H., "Uber die Partiellen Differengengleichungen der Mathematischen Physik," Math. Annalen, 100, pp. 32-74 (1928). 7. Douglas, J., and Peaceman, D. W., "Numerical Solution of Two Dimensional Heat Flow Problems," AIChE Journal, 1, p. 505 (December, 1955). 8. Drakhlin, E., "Thermal Convection in a Spherical Cavity," Zhur. Tekh. Fiz., 22, pp. 829-831 (1952). 9. Duncan, A. J., "Quality Control and Industrial Statistics," Richard D. Irwin, Inc., Homewood, Illinois, p. 492 (1953). 10. Reference 9, p. 72 (sum of approximate quantities). 11. Eckert, E. R. G., "Introduction to the Transfer of Heat and Mass," McGraw-Hill, New York, pp. 158-171 (1950). 12. Eddy, R. P., "Stability in the Numerical Solution of Initial Value Problems in Partial Differential Equations," Naval Ordnance Laboratory Memorandum 10232 (1949). 13. Griffiths, E., and Davis, A. H., "The Transmission of Heat by Radiation and Convection," Food Investigation Board Special Report No. 9, His Majesty's Stationary Office, London 1922). 14. Herman, R., "Heat Transfer by Free Convection from. Horizontal Cylinders in Diatomic Gases," NACA-TM-1366. Translated by German, "Warmeubergang bei freier Stromunglam wagrechten Zylinder in zweiatomigen Gasen," V. D. I. Forschungsheft, No. 379 (1936). -159

15. Hildebrand, F. B., "Introduction to Numerical Analysis," McGraw-Hill, New York (1956). 16. Jakob, M., "Free Convection Through Enclosed Plane Gas Layers," Trans. A. S. M. E., Vol. 68, pp. 189-194. 17. Jakob, M., "Heat Transfer," John Wiley and Sons, New York, p. 538 (1949). 18. Jones, C. D., and Masson, D. J., "An Interferometric Study of Free Convection Heat Transfer from Enclosed Isothermal Surfaces," Trans. A. S. M. E., 77, pp. 1275-1281 (November, 1955). 19. Lietzke, A. F., "Theoretical and Experimental Investigation of Heat Transfer by Laminar Natural Convection Between Parallel Plates," NACA-TN, 3328. 20. Mull, W., and Reiher, H., "Der Warmeschutz von Luftschichten," Beihefte zum Gesundheits-Ingenieur Reihe 1, Heft 28, Munich and Berlin, Germany, (1930). 21. O'Brian, G. G., Hyman, M. A., Kaplan, S., "A Study of the Numerical Solution of Partial Differential Equations," Jour. of Mathematics and Physics, 29, p. 223 (1951). 22. Ostrach, S., "An Analysis of Laminar Free-Convection Flow and Heat Transfer from a Flat Plate Parallel to the Direction of the Generating Body Force," NACA Rept. 1llll, (1953). 23. Ostrach, S., "A Boundary Layer Problem in the Theory of Free Convection," Ph.D Dissertation, Brown University (August, 1950). 24. Ostrach, S., "Laminar Natural Convection Flow and Heat Transfer of Fluids with and without Heat Sources in Channels with Constant Wall Temperatures," NACA-TN 2863 (December, 1952) 25. Ostrach, S., "Combined Natural- and Forced-Convection Laminar Flow and Heat Transfer of Fluids with and without Heat Sources in Channels with Linearly Varying Wall Temperatures," NACA-TN 3141. 26. Ostroumov, G. A., "A Quantitative Optical Method of Observing Heat and Diffusion Phenomenon for Two Dimensional Problems and Changes in Thickness of Plane Surfaces (The Grid Method)," Doklady Akad. Nauk, 71, pp. 887-889 (1950). 27. Parodi, M. M., "Sur quelques proprietes des valeurs caracteristiques des matrices carrees," Memorial des Sciences Mathematiques, Fascicule CX VIII, Gauthier-Villars, Libraire du Bureau des Longitudes, de L'Ecole Polytechnique, Quai des Grandes-Augustins, 55, Paris, France (1952).

-16128. Peck, R. E., Fagan, W. S., and Werlein, P. P., "Heat Transfer Through Gases at Low Pressures," Trans. A. S. M. E., 73, p. 281 (April, 1951). 29. Perry, J. H. (Editor), "Chemical Engineers Handbook," McGrawHill, New York, (1950). 30. Price, P. H., and Slack, M. R., "Stability and Accuracy of Numerical Solutions of the Heat Flow Equation," Brit. Journal of App. Physics, 3, pp. 379-384 (1952). 31. Robinson, H. E., Powlitch, F. J., "The Thermal Insulating Value of Air Spaces," Housing Research Paper No. 32, U. S. Government Printing Office, Washington 25, D. C. (April, 1954). 32. Rutkowski, J., "Free Convection from Heated Surfaces —Laminar Boundary Layers," Wayne Engineering Research Institute, Detroit, Michigan, (September, 1955). 33. Schmidt, E., and Beckman, W., "Das Temperature —und Geswindigkeitsfeld vor einer Warme abgebenden senkrechter Platte bei naturlicher Konvektion," Tech. Mech. u. Theramodynamik, Bd. 1, Nr. 10, (Okt. 1930). S. 341-349; cont. Bd. 1, Nr. 11 (Nov. 1930), S. 391-406. 34. Shaposhnikov, I. G., "On the Theory of Weak Convection," Zhur, Tekh. Fiz., 22, pp. 826-8 (1952). 35. Young, D., Lerch, F., "The Numerical Solution of Laplaces Equation on ORDVAC," Memorandum Report No. 708, Ballistic Research Laboratories, Aberdeen Proving Ground, Maryland. 36. Zhukhovitskii, E. M., "Free Stationary Convection in an Infinite Horizontal Tube," Zhur, Tekh. Fiz., 22, pp. 832-835 (1952).

APPEN DIX A FL( DATA The flow tracings for experiments 3 and 5 to 15 are presented on the following pages. The light lines represent flow lines. The two heavy lines are the locus of the tips of the gas velocity vectors drawn from the vertical or horizontal diameter. These tracings are derived from a careful reading of the priamary data which are photographic prints like Figure 13. -162

-163Exp.3 Exp.6 Exp. 7 Exp. 8 Scale for Velocity Vectors l,,,, I in/sec 0 5 Figure 66 Experimental Flow Lines and Velocities (Experiments 3, 6, 7, and 8)

Exp. 9 Exp. 10 Exp. II -Exp. 12 Scale for Velocity Vectors I I I,, I in/sec 0 5 Figure 67 Experimental Flow Lines and Velocities (Experiments 9, 10, 11, and 12)

-165Exp. 13 Exp. 14 Exp. 15 Scale for Velocity Vectors I, I,, I in/sec 0 5 Figure 68 Experimental Flow Lines and Velocities (Experiments 13, 14, and 15)

APPENDIX B TEMPERATURE PROFILE DATA AND RESULTS OF COMPUTATION This appendix contains all the original temperature profile data for all experiments. The position of a temperature measurement in this data can be found from the traverse letter and the point number by consulting Figure 11. The results of computation are given for all experiments except Experiment 4. The results of computation for Experiment 4 are given on Table IV. -166

-167TABLE XII. TEMPERATURE PROFILE DATA FOR EXPERIMENTS 3, 4, 5, AND 6 Exp. No. 3 Temperature Survey, *C Exp. No. 4 Temperature Survey,'C Ave. Thermo- Point Traverse Ave. Thermo- Point Traverse couple Readings No. A B C D E F couple Readings No. A B C D E F T. C. Temp. T.C. Temp. Number 1 43.8 ---- 43.5 43.5 52.7 63.9 Number 1 28.2 28.4 9.4 27.9 49.2 62. 2 1 2 ---- 45.8 44.9 43.0 53.8 62.6 1 2 32.5 31.9 32.7 27.9 49.3 60.6 2 41.5 3 46.0 47.5 46.2 48.0 56.5 ---- 2 24.5 3 35.6 35.3 35.8 31.7 54.6 58.8 3 4 47.8 49.5 48.0 49.5 ---- 60.8 3 4 39.2. 40.0 37.0 35.8 ---- 57.4 4 9.0 5 52.0 52.7 49.6 54.3 ---- ---- 4 13.25 5 45.2 ---- ---- 45.5 56.2 54.8 5 10.7 6 55.8 58.5 50.2 57.0 58.8 56.5 5 17.25 6 49.5 44.5 39.9 48.5 54.6 50.7 6 7 57.9 58.0 49. 7 57.4 57.0 54.0 6 7 50.3 43.7 38.6 48.3 51.7 49.0 7 8 57.4 57.0 49. 6 54.8 55.8 52.7 7 8 50.4 43.5 38. 1 46.6 49.7 46.6 9 9 57.2 ---- 49.3 52.0 53.5 51.2 9 9 50.4 43.7 37.7 43.0 46.5 43.2 10 10 56.5 57.2 49. 3 50.0 51.5 48.8 10 10 50.0 43.7 37.3 39.9 43.0 40.0 11 11 56.1 57.2 49.2 49.8 49.5 46.7 11 11 49.5 44.0 36.7 38.7 39.2 35.3 12 65.5 12 56.1 57.7 50.5 48.5 48.0 47.5 12 63.4 12 48.7 44.0 37.5 36.7 37.3 36.0 13 13 58.1 57.5 54.1 46.8 47.3 52.3 13 13 ---- 44.0 41.9 33.9 36.0 40.6 14 14 60.4 57.2 56.0 45.6 46.2 58.2 14 14 52.5 43.4 48.3 30.7 32.8 49.2 15 15 61.7 55.6 62.0 45.0 46.2 60.3 15 15 54.6 45.5 53.9 29.4 31.7 55.8 16 16 63.7 58.8 64.0 44.4 ---- 63.7 16 16 58.6 50.0 59.1 27.7 ---- 60.8 17 17 64.4 60.9 66.4 43.5 ---- 63.7 17 17 61.6 54.3 62.2 26.2 33.7 62.4 19 18 63.9 49.2 19 18 57.4 36.3 20 19 64.7 49.8 20 19 65.7 38.5 Exp. No. 5 Temperature Survey,'C Exp. No. 6 Temperature Survey,'C Ave. Thermo- Point Traverse Ave. Thermo- Point Traverse couple Readings No. A B C D E F couple Readings No. A B C D E F T. C. TempT. C. Temp. Number 1 19.8 20.2 20.6 20.7 34.6 40. 0 Number 1 22.3 22.8 23.8 22.2 31.4 36.8 1 2 21.0 22.0 21.8 21.2 35. 3 40.5 1 2 23.0 23.8 24.5 22.5 32.2 36.8 2 16.5 3 22.9 23.4 22.5 23.5 36.5 38.7 2 19.5 3 23.8 26.0 25.2 24.0 33.6 36.3 3 4 25.9 26.0 23.8 26.1 37. 0 37.9 3 4 26. 1 26.7 26. 7 26.4 34. 1 36.2 4 8.0 5 29.7 28.9 25.4 31.5 37.0 36.5 4 13.9 5 28.7 28.9 28.0 29.7 34.9 35.5 5 10.8 6 32.2 30.2 27.5 35.0 36.3 35.0 5 16.8 6 31.2 30.6 28.8 32.6 34.6 34.6 6 7 33. 1 29.9 26.2 35. 1 35. 3 33.0 6 7 32.6 30.8 28.9 33.1 33.6 33.5 7 8 33.2 29.9 26.2 32. 7 34. 0 32.2 7 8 32.6 30.2 28.6 31.9 32.9 32. 7 9 9 33.2 29. 1 26.2 30.9 32.2 30.4 9 9 32.6 30.3 28.4 30.9 32.0 31.6 10 10 32.7 29. 1 26.0 29.0 30.6 28.5 10 10 32.4 30.3 28.3 29.6 31.2 29.6 11 11 32.5 27.7 25.7 28.2 28.9 26.3 11 11 32.4 30.3 27.8 28.9 29.5 27.9 12 42.2 12 32.5 27.7 26.7 29.4 27.1 26.0 12 40.2 12 32.4 30.3 28.9 27.7 28.1 28.3 13 13 34.1 27.8 29.4 24.5 25.7 30.0 13 13 33.6 30.3 31.5 25.5 27.5 30.7 14 14 36.0 29.7 32.4 22.7 24. 0 35. 1 14 14 34.6 30.9 34. 3 24. 1 25. 7 35.8 15 15 37.0 31.2 35.8 21.5 23.1 39.3 15 15 36.1 31.8 36.3 22.9 24.5 38.1 16 16 38.5 33.4 38.4 20.7 23.7 41.9 16 16 36.7 33.9 38.2 22.2 24.5 39.7 17 17 39.9 35.8 39.9 20.2 25.0 42.7 17 17 38.0 35.6 39.3 21.5 25.5 39.8 19 18 38.1 26.3 19 18 37.5 26.4 20 19 40.2 27.4 20 19 39. 0 27.5

-168TABLE XIII. TEMPERATURE PROFILE DATA FOR EXPERIMENTS 7, 8, 9, AND 10 Exp. No. 7 Temperature Survey,'C Exp. No. 8 Temperature Survey, eC Ave. Thermo- Point Traverse Ave. Thermo- Point Traverse couple Readings No. A B C D E F couple Readings No, A B C D E F T. C. T. C. Number Temp. 1 16.75 18.6 18.0 17.9 21.25 22.9 Number 1 16.28 15.98 16.40 ----- 17.78 17.70 1 2 17.2 18.6 18.3 17.9 20.3 22.9 1 2 16.28 16.08 16.80 15.40 17,70 17.55 2 17.12 3 17.2 18.8 18.7 18.3 20.5 22.85 2 16.20 3 16.42 16.20 16.68 15.40 17.70 17.30 3 4 17.2 19.2 19.3 18.7 21.0 22.7 3 4 16.49 16.25 16.78 15.82 17.74 17.20 4 12.3 5 17.2 19.7 19.3 19.7 21.0 22.6 4 9.50 5 16.82 16.56 16.95 16.12 17.88 17.07 5 13.25 6 18.95 20.5 19.8 20.6 21.5 22.1 5 11.30 6 17.20 16.95 17.25 16.60 18.15 16.88 6 7 20.6 20.9 20.3 21.5 21.25 22.1 6 7 17.92 17.25 17.62 17.25 18.15 16.45 7 8 20.75 21.0 20.3 21.3 21.25 21.85 7 8 18.20 17.38 17.55 17.18 17.88 16.30 9 9 21.25 21.0 20.2 21.0 21.1 21.25 9 9 18.25 17.38 17.40 17.00 17.70 16.10 10 10 22.0 21.0 20.0 20.3 21.1 20.85 10 10 18.22 17.41 17.30 16.88 17.42 15.71 11 11 22.0 21.0 19.7 20.0 20.85 20.6 11 11 18.18 17.41 17.18 16.75 17.18 15.55 12 22.84 12 22.0 21.25 20.3 18.8 20.6 20.6 12 18.13 12 18.22 17.38 17.40 16.30 17.05 15.86 13 13 22.5 21.25 21.1 18.4 20.3 21.85 13 13 18.42 17.30 17.70 16.12 16.90 16.42 14 14 22.7 21.25 21.6 17.9 19.7 22.85 14 14 18.56 17.38 18.00 15.92 16.35 17.05 15 15 22.85 21.5 22.2 17.7 19.7 23.3 15 15 18.56 17.69 18.15 15.80 16.40 17.30 16 16 22.85 22.0 22.5 17.7 19.8 23.7 16 16 18.65 17.95 18.20 15.52 16.53 17.55 17 17 22.85 22.7 22.5 17.7 19.8 23.7 17 17 18.70 18.05 18.42 15.40 16.53 17.56 19 18 22.6 20.0 19 18 18.20 16.53 20 19 22.8 20.6 20 19 18.18 16.67 Exp. No. 9 Temperature Survey,'C Exp. No. 10 Temperature Survey, ~C Ave. Thermo- Point Traverse Ave. Thermo- Point Traverse couple Readings No. A B C D E F couple Readings No. A B C D E F T. Temp. T.C. Temp. Number 1 25.6 26.8 26.3 41.0 50.9 65.0 Number 1 35.3 37.3 36. 1 33.8 87.4 112. 2 1 20.60 2 25.7 29.7 29.3 42.2 52.7 64.5 1 2 37.1 44.3 43.8 33.8 91.2 107.5 2 20.62 3 29.3 34.8 32.6 45.5 56.3 62.9 2 26.10 3 46.5 54.9 52.4 45.6 98.8 103.3 3 21.1 4 35.0 38.7 35.5 49.4 57.0 61.0 3 4 59.3 61.5 52.5 60.3 98.9 98.3 4 7.60 5 44.2 44.2 39.0 53.0 56.4 58.0 4 7.50 5 73.2 68.7 58.7 77.0 98.3 90.3 5 11.60 6 51.1 45.6 40.6 55.4 53.3 54.2 5 14.10 6 81.4 68.7 59.9 89.0 92.2 84.8 6 21.0 7 52.9 44.0 39.3 53.9 52.9 51.9 6 30,3 7 82.7 b6.5 57.3 85.7 86.9 82.1 7 21.75 8 53.,0 44.4 38.9 48.9 54.6 50.9 7 29.85 8 83.2 67.0 57.3 76.0 82.8 77.4 9 37.6 9 52.6 44.4 38.9 45.3 48.0 46.5 9 9 83.0 67. 1 56.5 68.9 76.0 69.9 10 22.0 10 52.3 44.9 39.0 42.0 45.8 43.5 10 10 82.3 67.7 56.0 62.4 69.2 64.6 11 71.3 11 51.8 45.3 39.0 39.8 42.5 40.0 11 11 81.8 66.5 56.2 58.4 63.4 56.9 12 72.47 12 51.9 45.6 40.0 38.5 39.7 37.9 12 116.62 12 80.0 66.9 56.4 55.6 58.0 51.2 13 71.4 13 54.7 46.1 44.2 35.1 38.6 44.0 13 13 84.9 67.6 63.0 49.8 50.0 64.6 14 14 59.6 45.6 53.3 31.5 37.0 55.9 14 14 93.5 67.0 77.9 44.7 48.1 84.8 15 15 61.8 47.4 58.6 29.0 34.8 64.6 15 15 99.5 69.2 88.1 40.6 44.3 99.0 16 *65.5 16 64.8 54.4 65.1 27.2 33.3 69.6 16 118.4 16 103.9 81.0 101.7 35.9 44.3 111.4 17 74.0 17 67.7 58.8 70.6 25.7 32.6 72,3 17 118.0 17 110.0 90.1 110.8 33.1 45.2 114.2 19 65.25 18 65.4 31.8 19 18 101.2 48.7 20 72.9 19 69.6 32.1 20 19 110.0 51.6

-169TABLE XIV. TEMPERATURE PROFILE DATA FOR EXPERIMENTS 11, 12, 13, AND 14 Exp. No. 11 Temperature Survey, ~C Exp. No. 12 Temperature Survey, ~C Ave. Thermo- Point Traverse Ave. Thermo- Poitut Traverse couple Readings No. A B C D E F couple Readings No. A B C D E F T. C, T.C. T.C. Temp T.C. Temp. 1 Number 1 30.8 30.8 33.0 30.5 74. 1 95.5 Number Tmp. 5.6 44.8 48.6 44.9 79.4 111.1 1 2 31.4 37.0 38.2 31.7 78.3 91.7 1 2 45.8 52.7 54.6 44.9 82.2 110.7 2 25.0 3 38.3 44.6 43. 1 47.4 83.7 88.9 2 36.80 3 55.2 61.9 59.9 46.0 92.3 109.8 3 4 47.0 52.6 47.2 51.0 84.8 81.5 3 4 64.8 70,1 66.2 60.8 99.7 108.4 4 8.0 5 62.8 59.7 51.6 67.1 84.7 78.7 4 8.0 5 80.9 81.6 73.6 80.3 102.6 103. 9 5 12.0 6 71.9 61.0 52.9 77.2 76.4 72.2 5 12.45 6 91.2 85.8 75.6 92.9 100.6 97.3 6 26.2 7 74.0 58.7 50.2 74.2 75. 1 69.3 6 41.25 7 94.9 82.0 72.9 94.8 96.0 92. 1 7 25.75 8 73.6 58.8 49.9 65.6 71.5 65.5 7 39.70 8 93.9 86.3 73.0 88.2 93.0 88.8 9 9 73.5 59.1 49. 2 59.8 65. 1 59.2 9 9 92.3 83.9 73. 0 81.3 88.4 83.9 10 10 70.9 59.3 49.7 54.2 59.2 54.1 10 10 90.8 84.5 73.0 76.5 82.8 77.8 11 11 69.9 59.7 49.9 50.8 53.8 48.6 11 11 92.,0 78.9 73.6 72.4 76.8 70.5 12 100.2 12 68.0 59.6 50.3 48.8 49.2 45.4 12 12 92.2 87.4 74.4 69.9 71.8 67.4 13 13 72,1 59,8 58.2 44.7 45.6 54.3 13 13 98.9 87.1 83.5 63.7 68.7 79.3 14 14 81.0 58.6 64.5 39.3 41.0 71.5 14 147.32 14 108.0 86.0 102.0 56.2 62.9 107.5 15 15 85.7 60.7 72.1 35.1 38.1 85.6 15 164.3 15 114.0 89;5 116.3 51.7 56.4 127.,1 16 100.0 16 91.2 68.6 85.5 31.1 38.8 94.7 16 110.0 16 118.3 101.6 129.2 47.2 54.8 146.5 17 99.7 17 97.0 77.3 93.4 29.2 41.1 97.5 17 17 122.4 109.4 137.0 43.0 55.2 147.5 19 18 86.5 46.4 19 18 119.6 58.2 20 19 95.0 47.5 20 19 130.0 59.4 Exp. No. 13 Temperature Survey, OC Exp. No. 14 Temperature Survey,'C Ave. Thermo- Point TraAve. Thermo- Point Traverse couple Readings No. A B C D E F couple Readings No. A B C D E F T. C. Tm.T.C. Temp. Tomp. Number 1 48.6 45.2 48.9 45.8 80.8 122.8 Number 1 71.8 66.2 71.8 66.0 147.0 184. 1 2 50.9 51.4 54.4 45.8 84.8 120.1 1 2 72.2 81.5 80.1 66.0 152.0 182.6 2 46.15 3 83.7 95.8 86.2 75.3 160.7 179.1 2 35.48 3 56.3 61.9 61.4 49.9 98.3 118.0 2 46.15 3 83.2 75.3 160 3 4 66.9 71.7 68.3 64.7 105.4 114.7 3 4 102.3 110.1 96.2 99.4 167.0 176.0 4 9.85 5 81.9 84,0 73.8 82.5 107.6 108.4 4 10.65 5 129.7 124.4 104.8 136.4 168.5 166.6 5 15.20 6 85.0 86.9 76.2 98.2 105.7 101.4 5 17.6 6 147.0 125.3 106.5 152.3 163.6 154.1 6 37.60 7 100. 1 83.3 74.4 98. 7 99. 7 96.8 6 49.95 7 152.2 120.6 101.8 151.2 152.7 126.8 7 37.65 8 98.4 83.9 74.2 90.6 96.8 92.8 7 49.1 8 150.0 120,.9 101.8 134.2 147.7 136.9 9 9 97.6 84.9 73.0 83.7 91.3 86.9 9 9 147.5 121.7 101.8 120.9 136.4 125.1 10 10 97.0 84.9 73.2.77.9 85.0 80. 6 10 10 145.6 122.6 100.9 109.2 124.0 112.9 11 11 96.9 85.5 72.6 73.8 79.9 73,9 1t 11 143.1 125.5 101.2 100.9 112.0 101.4 12 146. 72 12 96.2 85.9 72.9 7b.6 74. 6 69. 3 12 220. 25 12 140.2 124.4 102.9 97.8 102.6 96.0 13 13 103.6 85.2 8-2.8 63.8 71.1 82.6 13 13 148.2 124.1 115.4 89.2 96.7 116.5 14 14 113.6 83.3 99.0 55.1 62.9 110.2 14 14 163.2 121.4 147.0 81.1 90.2 153.1 15 15 118.9 87.8 113,9 50,5 57.4 129,0 15 15 171.7 128.7 164.3 73.2 82.4 186.5 16 119.6 16 125.1 102.0 125.6 45,8 56,9 141.7 16 187.6 16 181.5 158.7 184.1 64.9 81.0 213.7 17 160.0 17 131. 1 109.4 137.1 41.6 60.5 144.4 17 223.9 17 189.0 168. 7 195.3 58.7 82.2 212, 6 19 18 122.3 68.2 19 18 189.0 90.1 20 19 132.4 68.8 20 19 205.7 99.1

-170TABLE XV. TEMPERATURE PROFILE DATA FOR EXPERIMENT 15 Exp. No. 15 Temperature Survey, ~C Ave. Thermo- Point Traverse couple Readings No. A B C D E F T. C. Number 1 81.1 82.5 83.9 82.6 181.0 220.6 1 2 81.6 98.5 94. 1 83.2 182.0 212.4 2 50.45 3 105.7 117.5 103.6 98.9 191.2 205.3 3 4 124.5 130.2 111.5 128.5 196.0 197.8 4 9.00 5 151.6 145.8 120.6 162.3 193.0 185.6 5 19.55 6 173.7 145.8 122.0 180.5 188.7 171.6 6 56.1 7 178.6 139.0 117.5 179.5 177.8 166.2 7 54.05 8 176.7 138.8 117.5 158.3 171.7 155.9 9 9 174.5 140.6 116.9 142.0 155.6 141. 0 10 10 169.8 141.8 116.8 127.8 140.8 127.8 11 11 167.1 142.9 117.0 118.4 127.0 114.9 12 254.10 12 162.2 143.0 118.0 113.7 114.1 106.6 13 13 172.3 144.2 131.9 102.9 106.3 139.7 14 14 189.8 139.7 159.6 90.2 95.4 183.1 15 15 202.7 146.1 185.3 81.6 88.0 215.2 16 225.6 16 214.0 161.5 205.5 73.8 91.3 239.5 17 251.6 17 227.5 173.3 223.9 64.8 102.7 240.5 19 18 185.0 111.7 20 19 200.5 112.1

ANALYSIS OF TEMPERATURE DATA FOR EXPERIMENT 3 ANALYSIS OF TEMPERATURE DATA FOR EXPERIMENT 5 Nusselt No. 50% Rate of 50% Conf. Heat Trans. 50% Conf. 50% Rate of 50% Conf. Heat Trans. 50% Conf. Slope Normalized 50% Conf. Limit Slope Normalized 50% Conf. Nusselt No. Heat Trans. Limit Coefficient Limit N D 8T Conf. Heat Trans. Coefficient Limit Lcmit by C.... Limit N D 8T Conf. by C.... BTU/hr ft2BTU/hr ft2 ~F No. Slope'C/in,'C/in. u = ~ -~r Limit BTU/hr ft2 BTU/hr ft2 BTU/hr ftZ'F BTU/hr ft2'F Slope ~C/in ~ in u = ~-~ 8r Limit BTU/hr ft2 ~F No. BTU/hr ft2 1 38.00 +9.24 6.82 + —1.66 13.00 +_3.16 0.301 +0.073 I 19.00 +-Z. 88 3.18 +0.48 6.07 +-0.92 0.131 +0.020 2 60.04 +-7.54 10.77 + 1.35 20.56 +2.58 0.476 +0.060 2 58.43 +1.13 9.78 +0.19 18.68 +0.36 0.404 +0.008 -- w _ -- -- -- -- [_.[ 3 29.31 +-3.91 5.26 +-0.70 10.04 +-1.34 0.232 +-0.031 3 45.55 +_5.26 7.63 4_0.88 14.56 +_1.68 0.314 +-0.036 ~ 4 37.00 +1.73 6.63 +-0.31 12.67 +-0.59 0.293 _~0.014 4 37.60 4.2.77 6.26 4.0.46 12.02 +-0.88 0.260 +_0.019 5 30.17 +-0.65 5.41 +-0.12 10.33 +-0.22 0.239 +-0.005 5 21.23 +-3.23 3.55 +-0.54 6.79 +-1.03 O. 147 +-0.022 6 9.97 +3.30 1.79 +0.59 3.41 4-1.13 0.079 + 0.026 6 20.08 4-1.00 3.36 +1.68 6.42 +-0.32 0.139 +0.007 Ave. C~ Ave. 33.60 + 1.42 6.03 +0.25 11.50 +0.49 0.266 + 0.011 35.31 +0.88 5.92 +0.15 11.29 +0.28 0.244 + 0,006 C.S. -- _ _ C.S. -- -- 7 -18.46 +-1.60 -3 31 +-0.29 -6.74 +0.58 -0.156 +-0.014 7 -24.00 +_1.15 -4.02 +-0. I9 -8.24 +0.40 -0.178 +-0.009 8 -63.08 +-2.37 -11.40 +0.43 -23.00 +0.86 *0.533 +0.0Z0 8 -68.35 +3.88 -11.45 4.0.65 -23.47 +-1.33 -0.507 +0.029 9 -59.15 +7.17 -10.61 +1.29 -21.59 +2.62 -0.500 +0.061 9 -49.72 +3.22 -8.33 +0.54 -17.07 +1.11 -0. 369 +0.024 O 10 -41.40 +-4.90 -7.42 +-0.88 -15.11 +-1.79 -0. 350 +_.0.041 10 -45.40 +-0.78 -7.60 +-0.13. -15.59 +-0.27 -0. 337 4-0.006 11 -31.57 +3.13 -5.66 +0.56 -11.52 +1.14 -0. 267 +0.026 11 -29.76 +1.40 -4.99 +0.23 -10.22 +0.48 -0.221 +0.010.... --.... O 12 -22.86 +-2.23 -4. 10 4-0.40 -8.35 +-0.82 -0.194 +-0. 019 12 -18.01 +_2.55 -3.02 +-0.43 -6.18 +-0.87 -0.134 4-0.019 Ave. Ave. 41.65 40.81 +0.90 -6.83 +0.15 -14.02 +0.31 =0. 303 +0.007 H.S. +-1.31 -7.47 +-0.24 -15.20 +-0.48 -0.362 +-0.011 H.S. -- -- -- -- BTU BTU k~t =0.36504hrft2.C/in AT 46.26*F D/AT = 0.16751 in/*C k~ = 0.31968 hr BTU k~4 = 0.34344 BTU AT = 43.2~F D/AT = 0.179375 in/'C k~ = 0. 34236 hr ftZ'C/in = I I I I I I ftz * C/in'hr ft2 * C/in -" H O ANALYSIS OF TEMPERATURE DATA FOR EXPERIMENT 6! ANALYSIS OF TEMPERATURE DATA FOR EXPERIMENT 7 ~ [.-a — 4 Slope 50% Conf. Nusselt No. 50% Rate of 50% Conf. Heat Trans. 50% Conf. ~ ~.j Normalized Heat Trans. Limit Coefficient Limit Slope 50% Conf. Nusselt No. 50% Rate of 50% Conf. Heat Trans. 50% Conf. Limit Normalized D 8T Conf. by Cony. No. BTU/hr ft2 BTU/hr ftZ*FBTU/hr ft2 ~F No. Slope "C/in. LcT. it D 8T Conf. Heat Trans [ Slope ~C/in. by Conv. Limit Coefficient Limit ~J Nu - AT 8r *C/in. Limit Nu = ~-~ mr Limit BTU/hr ft2 BTU/hr ft2 xn. BTU/hr ft2 BTU/hr ft2~ BTU/hr ft2 *F I 22.00 +-3.46 4.57 +-0.72 7.06 4-1.11 0.189 +_0.030 I 7'.50 +_6.64 5.64 +-5. O0 2.39 +_2.11 0.232 +_.0.205 2 41.36 +-1.81 8.60 +_0.38 13.27 +_0.58 0.356 +_0.016 2 10.18 +0.36 7.66 +0.27 3.24 +-0.12 0.315 _ +-0. 011 ~J 3 27.51 +5.00 5.72 +-1.04 8.82 +-1.60 0.237 4.0.043 3 3.04 +_.1.43 2.29 +-1.08 0.97 +-0.46 0.094 +-0.044 4 27.80 +-3.21 5.78 +_0.67 8.92 +__1.03 0.239 +-0. 028 4 4.00 +-1.03 3.01 +__0.78 1.27 +-0.32 0.124 4.0. 032 5 15.64 4-0.00 3.25 4-0.00 5.02 +_.0.00 0. 135 +_0.000 5 7.82 +-0.65 5.89 4-0.49 2.49 +-0.21 0.242 +_0. 020 i~ 6 18.53 +_0.74 3.85 +-0. 15 5.94 +-0.24 0.159 +_0.006 6 3.23 +_0.66 2.43 +-0.49 1.03 4.0.21 0.100 +__0.020 C~ Ave. 26.05 + 0.84 5.42 + 0.17 8.36 + 0.27 0. 224 + 0. 007 Ave. C.S. -- -- -- -- C.S. 4.74 +_0.42 3.57 +-0.32 1.51 +-0.13 0.147 +-0.014 7 -20.00 +__1.15 -4.16 +-0.24 -6.83 +-0.39 -0.183 +_.0.011 7 -8.00 +-2.31 -6.02 +-1.74 -2.59 +_0.75 -0.252 +_0.073 8 -51.60 4-2.71 -10.73 +-0.56 -17.61 +-0.92 -0.473 +-0.025 8 -10.50 +-0.30 -7.90 +_0.23 -3.40 +-0.98 -0.330 +-0.010,,o 9 -36.28 +-1.33 -7.54 +_0.28 -12.38 4-0.46 -0.332 +-0.012 9 -6.95 4-0.78 -5.23 +-0.58 -2.25 +_0.25 -0.219 +-0.024 C~ 10 -34.40 +-0.86 -7. 15 +0. 18 -11.74 +-0.29 -0.315 4-0.008 10 -4.60 4-1.81 -3.46 4-1.36 -1.49 4.0.59 -0.145 4-0.057'~' 11 -24.35 +_.2.16 -5.06 +_0.45 -8.31 +-0.74 -0.223 +-0.020 11 -1.01 +0.48 -0.76 +0.36 -0.33 +0.15 -0.032 +0.015 lZ -6.12 +-0.61 -1.27 +_0. 13 -2.09 4-0.21 -0.056 +-0.006 12 -1.50 +__0.20 -1.14 4-0.15 -0.49 4-0.06 -0.047 +__0.006 Ave. -29.48 +0.57 -6.13 +0.12 -10.06 +0.19 -0.270 4-0.005 Ave. H.S. -- -- -- -- H.S. -5.25 +-0.28 -3.95 +__0.21 -1.70 + —0.09 -0.165 +__0.009 I I I! 32076 BTU kI.I 0. 34128 BTU AT = 10. 296*F D/AT =0.752622 in/*C k~ 0. 3186 ~ BTU BTU AT 37.26~F D/AT = 0.207971 in/*C ~C = 0. = = ktI.I = 0.324 hr ftZ *C/in l~r ft2 ~C/in hr ft2 *C/in hr ft2'C/in I I

ANALYSIS OF TEMPERATURE DATA FOR EXPERIMENT 8 ANALYSIS OF TEMPERATURE DATA FOR EXPERIMENT 9 ~-~ Slope 50~, Conf. N.... It No. Normalized 50~ Rate of 50~, Conf. Heat T.... 50~ Conf. 50~ Conf. N.... It No. 50~, Rate Tf 50~ Conf. Heat T..... 50~ Conf. Heat Trans. Heat..... Limit Coefficient Limit Limit Coefficient Limit Slope Normalized ~f Conf. by Cony. Limit Nu - Limit, D 8T Conf. by Cony. BTU/hr It2 BTU/hr ftZ'FBTU/hr ft2 ~F No. Slope'C/in. ~T'C/in. BTU/hr ft2 No. Slope'C/in, ~C/in. 1,4u: ~-~- ~rr Limit BTU/hr ft2 8r BTU/hr fta"FBTU/hr ft2 ~F Limit BTU/hr ftZ I 0.80 + —0.46 1.78 +1.03 0.25 +_0.15 0.073 +_0.042 I 54.00 +_10.39 4.48 +-0.86 17.38 +-3.34 0.186 +-0.036 4.45 F-t 4-0.76 9.91 +-1.69 1.41 +-0.24 0.406 +_.0.069 2 60.88 +- 6.46 5.05 +-0.54 19.59 +-2.08 0.210 +-0.022 [.__[ 3 1.74 4-0.34 3.87 4-0.76 0.55 +-0.11 0.159 4-0.031 ~ 71.70 +_16.39 5.93 +_1.36 23.08 +-5.27 0.247 +-0.057 * 4 1.86 +-0.17 4.15 +-0.37 0.59 4-0.05 0.170 +_.0.015 4 81.60 +- 4.70 6.78 +-0.39 26.26 +-1.51 0.281 +-0.016 5 3.13 +3.35 6.98 4-7.48 0.99 4-1.07 0.286 4-0.307 5 70.39 +- 1.94 5.84 +__0.16 22.65 4-0.62 0.243 +-0.006 6 3.36 4-0.31 7.50 +-0.69 1.07 4-0.98 0.306 +-0.028 6 43.42 +_ 1.32 3.61 +-0.11 13.97 +0.42 0.150 4-0.005 ~/~ Ave. C.S. 2.72 4.0.38 6.07 +0.86 0.86 70.12 0.249 Ave. 70.035 63.20 7 2.50 5.25 +0.21 20.33 +0.80 0.218 +0.008 -- -- -- C.S. -- -- -- -- 7 -1.40 4.0.81 -3.12 +__1.80 -0.45 +__0.26 -0.129 +-0.074 7 +5.00 4. 6.35 +0.42 +-0.53 +1.86 +-2.36 +0.020 +0.025 8 - -6.39 4.0.24 -14.25 +-0.53 -2.04 +-0.76 -0.588 +-0.022 8 -147.59 7 6.98 -12.25 +__0.58 -55.00 72.60 -0.589 4.0.028 -- -- 0 9 -3.08 4.0.19 -6.88 4.0.41 -0.99 4.0.06 -0.284 +-0.017 9 -119.45 +- 3.89 -9.92 +-0.32 -44.51 +-1.45 -0.477 4-0.016 10 -1.68 4-0.40 -3.75 4-0.89 -0.54 4-0.13 -0.155 +-0.037 10 -104.40 +_ 5.33 -8.67 +-0.44 -38.90 +-1.99 -0.417 4-0.021 11 -1.15 4-0.22 -2.57 +-0.50 -0.37 +-0.07 -0. 106 +-0.021 11 -61.56 i- 2.35 -5.11 +-0.20 -22.94 +-0.88 -0.246 +-0.009 O 12 -2.90 4-0.37 -6.47 +_0.82 -0.93 +-0.12 -0.267 4-0.034 12 -33.18 +- 2.03 -2.76 +__0.17 -12.36 +_0.76 -0.132 +__0.008 Ave. -3.07 70. l0 -6.85 70.23 -0.98 70.03 -0.283 +0.010 Ave. H.S..... -85.074. 1.50 -7.06 4.0.12-31.70 +-0.56 -0.340 +_0.006 F-~ H.S. hr ft2'C/in I I [ I AT = 93.33 ~F D/AT = 0.083028 in/~C k~ = 0.32184 ~ k'i.I = 0.37260 BTU hr ft2'C/i —-— ~~ hr ft2 ~C/in hr ft2 ~C/in o ANALYSIS OF TEMPERATURE DATA FOR EXPERIMENT 10 ANALYSIS OF TEMPERATURE DATA FOR EXPERIMENT 11 ~j o FO 50~ Conf. Nusselt No. Rate of 50~, Conf. Heat Trans. 50~ Conf. 50~ Rate of 50~ Conf. 50~ Heat Trans. 50~, Conf. ~[~! Heat Trans. Slope 50~ Conf. Nusselt No. Heat Trans. Slope Normalized Normalized Limit Coefficient Limit. D 8T Conf. by Conv. No. N D 8T Conf. by Cony. Limit BTU/hr ft2 BTU/hr ft2~l Slope'C/in ~C/in. 5TU/hr ftZ ~F ~J ~-{ ~rr Limit BTU/hr ft2 BTU/hr it~ BTU/hr ft2~iBTU/hr ft2 ~F No. Slope ~C/in. Limit Limit Coefficient Limit au = ~-~ 8r Limit BTU/hr ft2 ~<~'C/in..... I 114.00 +-21.94 5.42 + —1.04 37.31 4- 7.18 0.229 +0.044 I 96.00 4- 6.93 5.50 +-0.40 31.31 4-2.26 0.231 +-0.017 2 224.74 +_17.75 11.64 +__0.84 80.09 +- 5.81 0.492 +-0.036 2 187.34 4-23.51 10.72 +-1.35 61.10 +_.7.73 0.451 +_.0.057 3 183.5 +-32.46 8.73 +-1.54 60.06 4-10.62 0.369 +-0.065 3 125.14 4.24.28 7.16 +-1.38 40.81 4.7.91 0.302 4-0.059 166.40 7 8.83 7.92 70.42 54.45 + 2.89 0.334 +0.018 4 146.00 + 4.79 8.36 70.27 47.62 +1.56 0.352 +0.016 _.... _ _ _ _ 5 182.12 4. 5.81 8.66 +-0.28 59.60 4- 1.90 0.366 +_0.012 5 112.85 4. 1.94 6.49 +-0.11 36.81 +_0.63 0.272 +-0.005 [-~ 6 77.97 +_ 4.14 3.71 +_0.~O 25.52 4. 1.36 0.157 +-0.008 6 73.16 +_ 3.44 4.18 +-0. Z0 23.86 +-1.12 0.176 +-0.008 Ave. Ave. C~ 165.53 7 5.36 7.87 +0.25 54.17 + 1.76 0.332 70.010 125.78 + 4.97 7.20 +0.28 41.02 +1.62 0.303 +0.012 C.S. -- -- -- -- C.S. -- -- -- -- x, 7 -64.00 4. 3.46 -3.04 +-0. 17 -27.44 +- 1.49 -0. 168 +-0.009 7 -64.00 +-24.25 -3:66 +__1.39 -25.57 4.9.69 -0.1889 +-0.072 k~) 8 -264.52 +-11.57 -12.58 +-0.55 -113.42 +_. 4.96 -0.696 +-0.030 8 -231.60 +-13.00 -13.26 4-0.74 -92.55 +-5.20 -0.684 +-0.038',,, 9 -219.46 +-11.34 -10.44 +_.0.54 -94.10 + 4.86 -0.578 +__0.030 9 -163.71 +-20.31 -9.37 +_.1.16 -65.42 4-8.12 -0. 483 +_0.060 ~_3 10 -196.20 +_. 5.02 -9.33 +-0.24 -84.12 4- 2.15 -0.516 +-0.013 10 -176.80 +_,. 1.48 -10.12 4-0.08 -70.65 4-0.59 -0.522 4-0.004 O 11 -121.53 4- 4.30 -5.78 +-0.10 -52.11 +- 1.85 -0.320 +0.011 11 -110.63 4- 3.27 -6.91 4-0.19 -48.20 +-1.30 -0.356 +-0.010 12 -99.06 4. 1.64 -4.71 +_.0.08 -42.47 4. 0.70 -0.261 +-0.004 12 -79.77 + 8.78 -4.57 +0.50 -3{.88 4-3.51 -0.235 4-0.026 Ave. Ave. -170.84 -146.97 + 4.02 -8.41 +0.23 -58.73 +1.60 -0.434 70.012 t 2.57 -8.13 4-0.12 -73.25 H.S. 4- 1.10 -0.449 4-0.007 H.S_ -- -- -- -- BTU AT = 162.9'F D/AT = 0.047569 in/'C k~ = 0. 32724 hr ft~'C/in k~ = 0.42876 AT = 135.36~F D/AT = 0.057247in/~C = 0. k~ = 0.39960 hi- ft2 ~C/in hr it2'C/in hr ft2 ~C/in I I BTU { BTU I!kc 32616 BTU I

ANALYSIS OF TEMPERATURE DATA FOR EXPERIMENT 13 ANALYSIS OF TEMPERATURE DATA FOR EXPERIMENT 12 Slope 50% Conf. Nusselt No. 50% Rate of 50% Conf. Heat Trans. 50~o Conf. Slope Normalized 50~o Conf. Nusselt No. 50~e Rate of 50% Conf. Heat Trans. 50% Conf. ~-~ Heat Trans.! Normalized Limit _ D 8T Conf. by Conv.' Limit Coefficient Limit Heat Trans. Limit Coefficient Limit LGi No. Slope *C/in. Nu = AT' ~r No. *C/in. Limit Slope *G/in it N D 8T Conf. by Gonv. in. u = A-'T 8 r Limit BTU/hr ft2 BTU/hr ft2*~BTU/hr ftZ *F BTU/hr ft2 BTU/hr ft2 BTU/hr ftZ*EBTU/hr ft2 *F BTU/hr ft2 1 129.00 +42. 15 5.02 +1.64 43.61 +-14.25 0.219 +0.072 1 175.00 +54.45 6. 77 +2.'12 58.78 +18.42 0.294 +0.092 2 213.37 +-26.83 8.31 +1.05 72.13 +- 9.07 0.363 2 217.04 +-19.25 8.40 +_.0.046 +__0.75 72.90 +- 6.47 0.364 +0.032 3 151.07 +-29.75 5.88 +1.16 51.06 +10.06 0.257 +0.051 3 136.10 +24.45 5.27 +0.95 45.70 + 8.20 0.228 + —0.041 H 4 170.20 +- 2.76 6.63 +-0.11 57.53 +__ 0.93 0.289 +-0.005 4 180.00 +-10.94 6.97 +__0.42 60.45 +- 3.67 0.302 +-0.018 H 5 126.25 +- 4.52 4.92 +-0.18 42.68 +- 1.52 0.215 +-0.008 5 139.66 + 9.67 5.41 +-0.37 46.91 + 3.25 0.234 +-0.016? 6 94.06 + 1.75 3.66 +-0.07 31.79 +_ 0.59 0.159 +0.003 6 100.76 + 0.72 3.90 +-0.03 33.84 + 0.24 0.169 +0.001..... Ave. Ave. 148.76 + 6.14 5.79 +0.24 50.29 + 2.08 0.253 +0.010 C.S. 155.74 + 5.66 6.03 +0.22 52.20 + 1.90 0.261 +0.010 C.S. -- -- -- -- -- 7 -42.00 +-10.39 -1.64 +-0.40 -20.28 + 5.02 0.102 +0.025 7 -83.00 +-40.99 -3.21 +-1.58 -39.98 +19.74 -0. 200 +-0.099 ~.~ 8 -377.26 +__19.61 -14.69 +-0.76 -182.13 +- 9.47 -0.915 +-0.048 8 -338.41 +__12.15 -13.10 +-0.47 -163.00 +- 5.85 -0.814 +-0.029 9 -249.19 +- 9.92 -9.70 +-0.39 -120.30 +- 4.79 -0.605 9 -248.85 +-11.71 -9.63 +-0.024 +-0.45 -119.87 +- 5.64 -0.599 +-0.028 10 -190.80 +- 7.17 ~ -7.43 +-0.28 -92.11 +- 3.46 -0.463 10 -208.20 +-11.84 -8.06 +-0.017 +-0.46 -100.29 + 5.71 -0.501 +-0.028 ~'~ 11 -107.10 +- 5.86 -4.17 +0.23 -51.70 + 2.83 -0.260 11 -132.36 + 2.49 -5.12 +-0.014 +_0.10 -63.75 +- 1.20 -0.318 + —0.006 ~ 12 -32.88 +- 4.18 -1.28 +-0.16 -15.87 +_ 2.01 -0.080 +0.010 12 -64.90 +_ 2.14 -2.51 +-0.08 -31.26 _ +- 1.02 -0.156 +0.005 ~j Ave. F.J -Ave. -180.79 + 3.70 -7.04 +0.14 -87.28 + 1.78 -0.439 +0.009 H.S. -189.48 + 3.58 -7.33 +0.14 -91.26 + 1.72 -0.456 +0.009 H.S. - - - - - - - - -~'C] BTU BTU BTU BTU aT = 198.936~F D/AT = 0.038952 in/~C k'C = 0.33804 k'H: 0.48276 - AT = 200.232~F D/AT = 0 03870 in/~C k'C = 0 33588 hr ft2 ~C/in I [ j hr ft2 ~C/in hr ft2?C/in'. k~t = 0-48168 hr ftz 0C/in ~ U~ ANALYSIS OF TEMPERATURE DATA FOR EXPERIMENT 14 ANALYSIS OF TEMPERATURE DATA FOR EXPERIMENT 15 F-J kJq b-q I [ o 50% Conf. Nusselt No' 50% Rate of 50~0 Conf. Heat Trans. 50~o Conf. 50% Conf. Nusselt No. 50~0 Rate of 50~o Conf. Heat Trans. 50~o Conf. ~ lope Normalized Heat Trans Limit Coefficient Limit Slope Normalized Heat Trans. Limit Coefficient Limit Limit. D 8 T Conf.' Limit 8T Conf. by Conv. by Cony. BTU/hr ft2 BTU/hr ftZOFBTU/hr ft2 ~F Slope'C/in. ~C/in. No. Slope ~C/in. ~C/in. l~u = ~- ~r Limit u = ~-~- ~rr Limit BTU/hr ft2 BTU/hr ft2 BTU/hr ft2~FgTU/hr ft2 ~F F~ I BTU/hr ftZ o I 137.00 +_21.36 3.39 +0.53 47.64 +_ 7.43 0.152 + —0.024 1 102.00 +-47.34 2.16 +- 1.00 35.52 +-16.49 0.097 +_0.045 2 411.72 +_28.78 10. 18 +-0.71 143. 18 +-10.01 0.457 +_0.032 2 457.46 +-30.04 9.67 +-0.63 159.34 +10.46 0.435 +0.029 L-~ 3 232.24 +-53.20 5.74 +_1.31 80.77 +- 1.85 0.258 + —0.059 3 347.91 +_65.17 7.35 +_i.38 121.18 +-22.70 0.331 +-0.062 4 292.00 +- 2.83 7.22 +0.07 101.55 +- 0.98 0.324 +-0.003 4 324.:20 +13.70 6.85 +-0.29 112.92 +- 4.77 0.308 +-0.013 5 160.89 + —14.19 3.98 +-0.35 55.95 +- 4.94 0.179 +-0.015 5 220.~1 +_ 4.52 4.65 +_0.95 76.66 + 1.57 0.209 +-0.004 ~J 6 141.66 +- 9.39 3.50 +-0.23 49.26 + 3.27 0.157 +0.010 6 171.91 +- 5.66 3.63 +_0.12 59.88 + 1.97 0.163 +-0.005 A~ve. 239.72 + 8.62 5.93 +0.21 83.36 + 3.00 0.266 Ave. +0.010 C.S. -- -- -- -- 288.31 +-10.18 6.09 + —0.22 100.42 + 3.54 0.274 + 0.010 C.S. ~3 7 -169.00 +- 6.35 -4.17 +_0.16 -103.30 +- 3.88 -0.330 +0.012 7 -93.90 +-49.60 -1.98 +-1.05 -63.00 +-33.30 -0.172 +-0.091 ~/] 8 -550.04 +-43.27 -13.60 + —1.07 -336.23 +-26.46 -1.073 +-0.084 8 -694.20 +-46. 14 -13.95 +-0.98 -466.22 +-30.99 -1. 270 +-0.085 9 -368.12 +_10.15 -9.10 +-0.25 -225.03 +__ 6.20 -0.718 +-0. 020 9 -419.27 +-20.98 -8.86 +-0.44 -281.65 +-14.10 -0.768 +-0.038 [~0 10 -322.60 +__23.58 -7.98 +-0.58 -197.20 +-14.41 -0.629 +-0.046 10 -251.40 +-12.58 -5.31 +__0.27 -168.88 +- 8.45 -0.461 +-0.023 x~ 11 -196.62 +- 5.51 -4.86 +-0.13 -120.19 +- 3.37 -0.384 +-0.011 11 -280.50 +- 5.24 -5.93 +-0.11 -188.43 +- 3.52 -0.514 +-0.010 12 -82.27 +- 5.17 -2.23 +-0.13 -50.29 +- 3.16 -0.160 +_0.010 12 -157.44 +_ 4.22 -3.33 +_0.09 -105.76 +_ 2.84 -0.289 +_0.008 (..k.) Ave. -293.78 + 7.64 -7.26 +0.19 -179.59 + 4.67 -0.573 Ave. +0.015 H.S. -- -- -- -- H.S. 347.30 +- 8.45 -7.32 + —0.18 -233.30 +_ 5.67 -0.637 +-0.015 BTU I [ [ BTU hr ft2'C/in hr ft2 ~C/in AT = 366.57~F D/AT = 0.021139 in/~C k~ = 0.3483 BTU k~ = 0.67176 hr f'~-~ ~C/in hr ft2' C/in AT:313.38-F! D/AT=0.024727in/~C [ k~:0.34776 BTU I k'H=0'61128

APPENDIX C

-175Whr oet rga A copy of the complete program can'be obtained from the Statistical Research Laboratory, Rackham Building, University of Michigaan, Ann Arbor. It was submitted as?roblem No. 1203 on July 23., 19~0o. it is made up on seven word. load cards ready f or Use on an IMi-6_~0 with 2,000 words of magnetic drim, storage. Location of Programn and Constants on the Dru~m The Program along with the stored parameters occupies storage locations 0000 to 1622. Table C-l gives the location and present contents for the variables Tp, q, for the constants that could be changed without involving reprogramaming. For purposes of scaling in a fixed -point machine., the variables have the following range:.937 < 1.063 given 0~ <p or q.2 expected 0 < <.5:e'.000229369 ReO-l ro.00016109 A e: 62918 Also for the -purpose of keeping track of`[ the decimal point on the IE1i11650, the quantities were scaled as shown in the following table: Ouant Location of Decimal Point

-176TABLE XIX. DRUM LOCATIONS FOR TBE PROGRAM Quantity Location Initial Value s 1020 0938653209+ 1030 0938653209+ 1040 0938b53209+ 1050 0938653209+ 1120 0938653209+ 1070 1061346790+ 1080 1061346790+ 1090 10b1346790+ 1100 1061346790+ All other from 1101 to 1120 1000000000+ 1121 to 1240 all 0000000000-. ~~' s ~~~~1241 to 1360 all 0000000000-'Elevations of Grid Points Above Center of Grid 1361 to 1470 between +.5 and -.5 f for T Equation 1527'j1000000000 + f for P Equation 1528 I Mora'tor 1000000000 + 1529 Punch 0000I,- Timax 159j Band00000I'- 5Imax 1530 0000000000'" LTALY*11 1588 0000000000Temporary Storage Locations (no Carry Over from 1 Iteration to Next) 1531 to 1534 0000000000Moni~ter Punch Band 2527 to 1534 Punch Band for Variables 17to 1584 Grid Radii 1611 0000000000 + 1612 1000000000 + 1613 2000000000 + 1614 3000000000) + 161-5 4000000000 + 1616 4200000000 + 1617 4400000000 + 1618 460000000 + 1619 4800000000 + 1620 5000000000 +

-177How to Use the Program The switches on the console should be set as follows: (a) Programnred switch to RUN. (b) Half cycle switch to RUN. (c) Control switch to RUN. (d) Display switch to PROGRAM REGISTER. (e) Overflow switch to SENSE. (f) Error switc.h to STOP. (i.) Place the Program deck-in the hopper of the 533. (2.) Place standard control panel in the 533 Read-Punch Unit. (In this program all input is by load card which is signaled by a twelve punch in colurnn one. Output is 8-10 digit words per card. This output foxrmat is signaled at this laboratory by a high order eight in the tenth register of the punch band.) (3.) Set storage entry switches to 70 1951 3000. (4~.) Press computer-reset key. ($)Press program start key. (b.) Press start key on the 533. (7.) When read hooper empties, press end of file key. (8.) Put one word load cards into hopper which may be needed to change the program or its parameters (put one punched nonload card at end of deck). (9.) Set Storage entry switches to 70 1951 0001 (10.) Press computer reset key.

- iy8(13.) When read hopper empties, press end of' file key. Computation should now begin. The program now iterates on the T equation in two passes as is outlined'in the thesis. At the end of each rass the monitor punch band as indicated in Table C-i is punched out. (Ignore locations 1531 to 1534 as they contain intermediate results.) At the end of both passes all internal points in the T grid have been changed once,, In the same way the program iterates on the P's., After each pass the corresponding q's are computed. After iteration is complete on both T's and P's and new q's, the values at grid points 22 to 29 are punched out first for the p grid and then for the T grid. pVoint 22 is in first position on the card, on 23 in second, etc. The program now begins over if the prograxmed switch is on run. If changes in the program are necessary", the programmed switch can be set up to stop and the the program Will halt just before beginning again. The program can then be restarted at 0001. Prgrm odifications (1.) Skipping the 1 equation, put (16 8002 0[432) in 0001. (2.) Skipping 5 equation and q calculation., put 20 1529 1524 in 0463., (3.) Skipping F calculation., put 20 1530 0624 into 0428 and 20 1530 1524 into 091-3. (4.) Other revisions: The program is made up of three subroutines — T equation,

-179Tsubroutine, OO$O subroutine., 0650 subroutine, 04&50 In aditio, the Tsubroutine must be entered twice with no p sboutine in between. Thep subroutine must be entered twice with no T subroutine in between., If the q subroutine is used, it must be used immediately after the V subroutine. Notice that all subroutine entries are in groups of two. With the above information the user can make up his own control program to transfer control among the subroutine and to provide for any output that may be desired.

APPENDIX D

8 Q~~~~~~~~~~~~~~ __.3 6 — IZ 0. 6 4 — 0 0 o I I g0 40 80 120 160 TEMPERATURE (0C) Figure 69 Calibration of Copper-Constantan Thermocouples

z 5 I0 Lu QO C reference junc. I ~~~~~~~~~~~~~~Leeds and Northrup a. ~~~~~~~~~~~~~~~~calibration D 4~~~~~~~~~~~~~~~~~~~~ o oo o a: H3 2 1 40 80 12016 TEMPERATURE (0C) Figure 70 Calibration of Chromel P-Alumel Thermocouple

,icn 05 (3) co)