THE UNIVERSITY OF MICHIGAN INDUSTRY PROGRAM OF THE COLLEGE OF ENGINEERING NATURAL CONVECTION INSIDE A CIRCULAR CAVITY Poa To Chu-.,"_~ r I A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in the University of Michigan Department of Chemical and Metallurgical Engineering 1961 December, 1961 IP-545

ACKNOWLEDGEMENT The author wishes to express his gratitude to Professors F. G. Hammitt and J. J. Martin, the co-chairmen of his Doctoral Committee, as well as the members of the Committee, Professors S. W. Churchill, B. A. Galler, V. C. Liu, and E. H. Young, for their guidance and valuable suggestions in completing this work. The author is grateful to Messrs. J. M. Deimen and R. H. Radius, research assistants, for their assistance in constructing the test facility, and in taking the experimental data. To Mr. R. Bare of Argonne National Laboratory, the author wishes to express his sincere thankfulness for devising an analog computer circuit. The analog computer time provided by the Argonne National Laboratory is also acknowledged. The financial support of the Phoenix Project from 1959 to 1961 and a Rackham Scholarship during the year 1960 granted to the author are also acknowledged. The assistance of the Industry Program of the College of Engineering in the final preparation and reproduction of the manuscript is greatly appreciated. ii

TABLE OF CONTENTS Page ii ACKNOWLEDGEMENTS.......................................... LIST OF TABLES............................................. v LIST OF FIGURES.............................................. ABSTRiBACT....................................................... I. INTRODUCTION........................................... II. LITERATURE SURVEY......................................... i lu.-al Convection From a Flat Plate.................. 2. Natural Convection From the Outside and Inside of a Cylinder,............................... 3. Natural Convection Inside the Enclosed Planes........ 4. Natural Convection Inside a Vertical Circular Cylinder....................................... 5. Natural Convection With Inter —al Heat Generation...... III. THE THEORETICAL CONSIDERATIONS........................... 1. General Considerations............................. 2. The Basic Equations.................................. 3. The Solutions for the Boundary Layer Type Flow Regime.................................... o. (a) The Formulism of the Integral Equations.......... (b) The Asymptotic Solutions........................ (c) The Solutions in Numerical Form.................. (d) General Remarks on Computer Program............. vi 1 4 4 6 7 8 11 14 14 1' 25 25 28 37 43 45 61 78 85 94 94 97 111 115 129 4. Solution for Flow Regime with Similarity. ao.oo.. O. (a) The Similarity Transformation.................. (b) The Eigenvalues and the Approximate Solution..... (c) The Similarity Solution Obtained from an Analog Computer.................................. 5. General Discussion of the Theoretical Results......... IV. THE EXPERIMENTAL PROGRAM................................o 1. General Considerations.............................. 2. The Experimental Facility............................ 3. The Experimental Pro *- e.e.................. 4, Observations........................................... V. CONCLUSIONS........................................... iii

TABLE OF CONTENTS (CONT'D) Page NOMENCLATURE................................................ o o 131 APPENDIX I. Derivation of the Boundary Layer Equations......... II. Detail Derivation of the Momentum and Energy Equations from the Boundary Layer Equations for Numerical Computation......................o o III. Brief Review of Lighthill's Solution for Boundary Layer Flow Regime..................... IV. Details of the Similarity Solutions...... o........ V. The Details of Digital and Analog Computer Program~. VI. The Physical Properties of Mercury as Functions of Temperature.............................. 133 139 153 156 164 169 BIBLIOGRAPHY........................................o............ o 176 iv

LIST OF TABLES Table Page I Summary of Computed Results for Boundary Layer Type Flow Regime..................................... 44 II Summary of Experimental Results..............o....... 116 v

LIST OF FIGURES Figure Page 1 Schematic Representation of the Problem Configuration e......................................... 15 2 Similarity Flow Regime With and Without Stagnant Portion.................................. 20 3 Schematic of Physical Variables................... 21 4 Comparison of the Asymptotic Solution to the Exact Solution for Large Prandtl Number.....o......o. 36 5 Consolidated Solutions in Boundary Layer Type Flow Regime,.................oo o 46 6 The Development of Boundary Layer (t1 = 1 x 105)..... 47 7 The Development of Boundary Layer (t1 = 107)o o...... 48 8 Velocity Coefficient, y, as Functions of Axial Distance (t1 = 105)................................... 49 9 Velocity Coefficient, y, as Functions of Axial Distance (t1 = 10)................................... 50 10 Velocity and Temperature Profiles at Various Axial Positions...................................... 51 11 Velocity and Temperature Profiles at Various Axial Positions................................ 52 12 Velocity and Temperature Profiles at Various Axial Positions............................ o... 53 13 Velocity and Temperature Profiles at Various Axial Positions..................................... 54 14 Velocity and Temperature Profiles at Various Axial Positions............................. 55 15 Velocity and Temperature Profiles at Various Axial Positions....................................... 56 16 The Eigenvalue Generating Function With Three Zeros... 70 17 The Solution for the Temperature With Different Eigenvalue.................................... 71 vi

LIST OF FIGURES (CONT'D) Figure Page 18 Comparison of Lighthill's Temperature Solution With the Present Solution.......................... 75 19 Comparison of Lighthill's Velocity Solution With the Present Solution.............................. 76 20 Nusselt Number as a Function of Rayleigh Number Similarity Flow Regime (2/a = 47.5)................... 77 21 Comparison of Temperature Function from Power Series Method and from Analog Computer............ 81 22 Comparison of Stream Function from Power Series Method and from Analog Computer...................... 82 23 Comparison of Axial Velocity from Power Series Method and from Analog Computer....................... 83 24 Solution of G'(T) from Analog Computer................. 84 25 Consolidated Theoretical Results...................... 86 26 Dimensionless Temperature Distributions for Various Prandtl Numbers from a Vertical Heated Flat Plate............................................ 89 27 Dimensionless Velocity Distribution for Various Prandtl Numbers from a Vertical Heated Flat Plate.... 90 28 Schematic of Experimental Set-Up...................... 98 29 Assembly of the Test Section.......................... 100 30 Assembly of the Thermocouple......................... 105 31 Summary of the Test Results........................... 118 32 Wall and Centerline Temperature Profile............... 124 33 Radial Temperature Profile........................... 125 34 Typical Heat Input Distribution....................... 127 Plate 1 Photographs of the Heaters............................ 103 2 Photograph of the Experimental Facility............... 112 vii

ABSTRACT The objective of the present work is the study of the fluid motion and heat transfer in a vertical cylindrical cavity, closed at the lower end and opened to a large reservoir on the top. The temperatures of the walls of the cavity and of the reservoir are maintained constant and uniform with the wall temperature greater than that of the reservoir. The variation of the fluid density with temperature is taken into consideration only as it affects buoyancy. The fluid is otherwise considered uniform and incompressible. With this configuration the motion of the fluid is generated by the density differential within the fluid, and heat is transferred through the walls of the cavity to the reservoir by the convection motion and conductive in the fluid. Based on the assumption of laminar flow, and with the approximation of boundary layer theory, a theoretical analysis has been made. Three possible modes of fluid motion are obtained from physical consideration, i.e., the boundary layer type flow regime, the similarity flow regime, and finally, similarity flow with a stagnant portion at the bottom of the cavity. From the equations of motion, it is shown that the parameters differentiating the flow regimes are the Rayleigh number, the Prandtl number, and the aspect ratio of the cavity. A set of solutions in numerical form has been obtained for the boundary layer type flow regime. This set of solutions covers a viii

range of Ra x a/f (Rayleigh number x the aspect ratio) from 5 x 103 7 to 1 x 10 and a Prandtl number range of 0.02 to 10.0. For the same flow regime, an analytical solution has been obtained for the case of very large Rayleigh number. This solution assumes an arbitrary Prandtl number. Due to the complexity of the governing equations, the solution for the flow regime with similarity has been obtained only for the fluids with large Prandtl number. Two solutions, one by an analytical method and the other by a numerical method, have been included in this paper. These two solutions are in agreement, despite the fact that they are derived from two entirely different methods. An approximate solution for this flow regime has been obtained by Lighthill(25) which gives good agreement with the present solution. However, the solutions from present work are useful in case more accurate results are required. An approximate solution by Lighthill(25) to fill the gap between the solutions for similarity flow regime and boundary layer type flow regime has also been quoted in the present paper for the sake of completeness. An experimental program has been conducted to attempt to validate the results of the analysis. A facility compatible to the physical conditions specified in the theoretical analysis employs water and mercury as test media. The results from the water tests agree reasonably well with the analysis, while the test results from mercury show considerable discrepancy. It is believed that for the facility used in the investigation laminar flow inside the cavity is almost impossible to attain for low Prandtl number liquids (as in the case of mercury). The turbulent ix

motion effectively reduces the rate of heat transfer, axld gives a Nusselt number much lower than that predicted by the theory. However, since no data has been available for natural convection heat transfer of mercury in a closed cavity under the present conditions, the experimental results from this investigation do represent a new contribution. x

CHAPTER I INTRODUCTION Flows which are generated entirely by the action of body forces on fluids with density variations due to heating (or cooling) are referred to as natural - or free - convective flows, heat transfer by natural convection is of both theoretical as well as practical interest. The governing equations of natural convection are non-linear, and are coupled through both temperatures and velocity. Even with some simplifying assumptions, it is still a challenging task to solve this set of simultaneous equations in order to obtain the numerical results for this phenomenon. In engineering practice, the use of this mode of heat transfer through the walls of hollow passages in turbine rotor blades for cooling is one of its (0) many applications. Recently, with the rapid development of nuclear power reactors, the natural convection process has become of even greater importance, because this means of heat transfer appears in some of the (2) many schemes for extracting the heat energy from an atomic pole(2) Very few theoretical investigations have been made for the natural convection heat transfer other than those which are restricted to very simple configurations, such as flat plates and horizontal or vertical cylinders. Experimental data are equally lacking. Water and air are almost the only two fluids being considered as test media to verify the limited theoretical results. Indeed, literature is available for natural convection

-2 - in more complicated geometries and less conventional liquids such as flows in channels, or using liquid metal as test fluid*, but those investigations are often semi-emperical or purely experimental. Since there are so many physical variables involved, and each of them may take a wide range of values, a few semi-emperical results are of dubious worth in the absence of theoretical considerations. Therefore, in order to answer some of the many pending questions concerning natural-convection heat transfer, such as the effect of a confining wall (closed geometry) and of liquid metals (characterized by their low Prandtl number), this phenomenon is analyzed herein. To be more specific, the problem is to study the natural convection and its associated heat transfer in a vertical, cylindrical cavity with fluids of different Prandtl numbers at a given constant wall to reservoir temperature difference. Only the case of laminar flow is analyzed herein. From a scientific point of view, this is reasonable since the laminar analysis will shed some light on the study of the more complicated case of turbulent flow. From an engineering point of view, this is equally desirable even though in the practical application the flow will often be turbulent. Study of the corresponding laminar flow. being more reliable, is essential to bringing out certain broad qualitative details of what happens. Thus, in the following chapters, it can be seen that a general physical analysis is made to distinguish possible modes of fluid flow, the equations of motion, which are compatable to physical conditions but * For instance, References, 3, 4, 5,furnish an excellent example.

-3 - are nontheless reasonably simplified to permit mathematical manipulations, are then set up, solutions for different flow regimes are obtained either numerically or analytically, and the validity of these solutions is discussed. Finally, along with the theoretical analysis, an experimental program is also described. Water is used to represent high Prandtl number fluids and mercury is regarded as a typical low Prandtl number liquid. Some existing experimental data are also compiled herein. Indeed, these experimental data are used to verify the validity of the assumptions made within the scope of the present analysis, and hence, the correctness of the theoretical predictions.

CHAPTER II LITERATURE SURVEY Generally speaking, the case of natural convection heat transfer in an open geometry, such as heat transfer from a vertical plate, or on the outside of a vertical or horizontal circular cylinder, has been extensively studied. Experimental data for various fluids are also available. However, the information of its counterpart - natural convection in a closed space, is comparatively much less. Since many ingenious approaches used in solving the cases with open geometry can also be applied to the solution of the cases in confined space, it is desirable to summarize some of these simpler cases. 1. Natural Convection from a Flat Plate Although it is believed that Lorenz(6) is the pioneer investigator for natural convection from a vertical plate kept at a constant uniform temperature, different from the ambient*, E. Pohlhausen obtained the first correct approximate solution of this case(7), He used the boundary layer simplification of the equations of motion and formally transformed them * Lorenz in 1.881 published his paper on this subject. He assumed that heat was transferred from the plate to a layer adjacent to it by conduction. Thus a density gradient is formed and heat is carried away by the flowing stream. The horizontal velocity in this layer was completely neglected and consequently the equations of motion were simplified. These assumptions were later proven inaccurate by Schmidt experimentally. However, it is quite amazing to find the resemblance between Lorenz's simplification and boundary layer theory which was published by Prandtl somewhat 50 years later.

-5 - into ordinary differential equations from which he obtained a numerical solution for fluid of Prandtl number 0.73. The theoretical results thus obtained were then verified experimentally by the data of Beckman and Schmidt(8) who measured both velocity and temperature distribution and found good agreement.* Later, O. A. Saunders(9) solved the same problem with a different approach. He assumed a polynomial solution for the above mentioned transformed equations and fitted the coefficients of the polynomial to a desired accuracy. A special feature of his solution is the removal of the Prandtl number limitation which appeared in Pohlhausen's solution. Saunders further verified his results by measuring the overall heat transfer coefficient in water and mercury. In both cases, good agreement was seen. It is to be noted that this is one of the very few experiments of natural convection heat transfer from a flat plate in liquid metal. Using the integrated boundary layer equations and assuming a quadratic form for the velocity and temperature profiles, Squire ( * The detailed description of Pohlhausen's solution has never been published. E. Schmidt quoted this solution in one of his papersU ) and stated briefly that the theoretical solution was verified by his experimental data. Recently, in a Boiling Heat Transfer Symposium held at Argonne National Laboratory on September 6, 1961, which both the author and Professor Schmidt attended, the author queried Professor Schmidt about Pohlhausen's solution. According to him, Pohlhausen studied Schmidt's data and estimated the slope of the velocity and temperature functions. From these boundary values, Pohlhausen was able to integrate numerically the ordinary differential equations from the equations of motion. The same difficulty is to be encountered in a later section of this paper. However, the author is fortunate to be able to find the correct values by a trial-and-error computation with the aid of an analog computer, which was not available at Pohlhausen's time.

-6 - obtained an approximate solution of Nusselt number as a function of Prandtl number (Pr) and Grashof number (Gr)o In spite of the simplified method, his solution did check with the known "exact" solutions to a reasonable accuracy. Later, Braun and Heighway(l4) modified Squire's method by introducing more terms, which are functions of Prandtl number, into velocity and temperature profiles. It is indeed true that the introduction of these functions does improve the accuracy of numerical results. However, the simplicity of Squire's method is losto With the aid of the large scale electronic computers, the obtaining of a numerical solution of this problem becomes more feasible. Thus, with boundary layer approximations, Ostrach(ll) and Sparrow(l2) obtained a set of numerical solutions corresponding to fluids of various Prandtl numbers. Hellums and Churchill(13) developed a method of solving the exact equations of motion and obtained. the solutions for different cases including non-uniform plate temperature. 2o Natural Convection from the Outside and Inside of a Cylinder A semi-empirical study of natural convection from a horizontal cylinder was first made by Ro Hermanno(15) His results were satisfactorily checked by Jodlbauer(l7) who measured the velocity and temperature distribution of air surrounding a heated horizontal cylindero Numerical solutions of natural convection inside a cylinder have also been obtained by Hellums(l3) which check with Martini's(l6) data reasonably wello Natural convection from a vertical cylinder is very similar to that of vertical plate if the curvature of the surface is not too big.(10)

-7 - Under this condition, Sparrow and Gregg(18) applied a set of "similarity transformations" to the boundary layer equations and obtained two simultaneous ordinary differential equations which were solved numerically. In this same paper, Sparrow pointed out that the overall heat transfer coefficient could be approximated by considering heat transferred through a thin layer adjacent to the wall by conduction. No experimental data were available for comparison*^ however, K. Pohlhausen*- used another set of transformation to solve the case where the cylinder wall temperature varied linearly with the axis(l9). Due to the simplified geometry, the possible similarity solutions for vertical flat plates and cylinders have been systematically summarized by Yang(20). He started with boundary layer equations and determined the necessary conditions under which the original partial differential, equations can be transformed into ordinary differential equations. From these mathematical requirements, the compatibility of physical boundary conditions can be asserted. This paper has thus established the existence of various similarity transformations with which the original equations are simplified and eventually solved either analytically or numerically. 3. Natural Convection Inside the Enclosed Planes A theoretical prediction of heat transfer and fluid flow between two long flat plates kept at different constant temperature was first made by Batchelor(22) in 1954. The essence of his paper is to determine which of the several different flow regimes occur at a given Rayleigh number * It might be interesting to compare this solution to Lorenz's solution. ** Not E. Pohlhausen mentioned previously.

-8 - (Pr x Gr) and aspect ratio (the ratio of length of the plate to the width of the enclosure, L/d). Batchelor did not completely solve the governing equations of motion, instead, he obtained several polynomials which could approximate the exact solutions at extreme values of Grashof number. Using the stability criteria, he was able to determine the coefficients of the polynomials and the range to which they were applicable. His prediction was then compared with the experimental data by Mull and Reiher(21) Although the predicted values are generally higher, the comparison is satisfactory. Later, Poots(24) modified his approach by employing an orthogonal expansion of the governing equations and developed a technique of nbtaining numerical approximation. Several sets of solution corresponding to different Rayleigh number and aspect ratio were thus obtained. The theoretical results from this paper are comparatively more reliable as verified by the experimental data(23). With the growing computer technique, this method has assumed increasing interest. 4. Free Convection Inside a Vertical Circular Cylinder f25) This subject was first studied by Lighthill 5) in 1.953 and subsequently by Leslie(26) and Liu(27). Lighthill studied the fluid motion inside a vertical circular cylinder, closed at one end and opened to a heat reservoir at the other end. The temperature of the cylinder walls is kept at constant, different from that of the reservoir. Thus, heat is transferred through the tube wall to th te reservoir by combined natural convection and conduction. A special feature of this problem is that,

-9 - due to the presence of the closed end, the fluid motion is made to form an internal boundary by itself. This is quite conceivable if one considers that the fluid adjacent to the solid wall is moving upward by buoyancy while the continuity requirement makes it necessary that fluid in the central portion of the tube has to flow downward. Thus, somewhere in between, there must exist a stationary layer. Indeed, it is shown by the equations of motion that three independent parameters, theRayleigh number, the Prandtl number and the aspect ratio, determine the fluid motion uniquely. Nevertheless, in addition to this, Lighthill predicted that fluid motion is not exclusively of only one type. He predicted that for large values ofRayleigh number (Ra) or small aspect ratio: (2/a), the flow would be of the boundary layer type. This type of flow ceases to exist when the first critical value of Ra x a/2, t1, is reached. Then the motion is of the "similarity type" by which it is meant that the velocity and temperature distributions are similar along the axial distance. The similarity regime again is terminated at the second critical value oftl, below which the fluid motion is stagnant in some portion of the tube. Using a modified Squire's method(l0), Lighthill obtained approximate solutions and determined the appropriate critical values of t1 for different flow regimes. Using the theory of eddy properties (exchange coefficient concept), he estimated the corresponding values for the case of turbulent flow. It should be noticed, however, that a partially linealizing assumption is made in his solution so that the analysis is applicable only to fluids of large Prandtl number. To be exact, the

-10 - Prandtl number of the fluid to which his solution could apply should be infinitely large. Physically, this corresponds to neglecting the inertia of the fluid as compared with the forces due to its viscosity and the external acceleration. In a paper published in 1957, Ostrach(28) extended Lighthill's similarity solution to the case of linear variation of wall temperature. Still later, Leslie(26) used Lighthill's method to solve the problem in the same geometry but with a first order perturbation term which is introduced by considering the external acceleration to be at a small angle to the axis of symmetry. Apparantly, his solution is intended to throw some lights on the "turbine blade" problem(l) where the Coriolis force is exerted on the moving fluid. In a paper by Liu(27) yet unpublished this same problem is attacked by a more sophisticated method. Liu modified Poots' method(24) of orthogonal expansion of the equations of motion and obtained solutions in a limited range. Since the boundary layer approximation is not involved in this case, his solution should be very accurate and more reliable. Experimental work on this subject has been done by Cohen(29), Martin(30) and Cresswell(31). As a whole, the experimental data check with Lighthill's laminar analysis to a reasonable accuracy, the error in the theoretical prediction of turbulent flow is high, however. Since in Lighthill's analysis a parameter of large Prandtl number has been assumed, the experimental work is mainly for conventional fluids. Nevertheless, the experimental data for air which has an average Prandtl number of 0.73 do agree with the theory in spite of this assumption. This fact was pointed out beforehand by Lighthill in a discussion of the applicability of his solution(25).

-11 - The experiments by Hartnett, et al., (32) were conducted in the same geometry except that the condition at the walls being constant flux, Both mercury and water were used as test media and the data were correlated in terms of average Rayleigh number (Gr x Pr) versus average Nussult number. It is reported that corresponding to a given Rayleigh number, Nussult's number for heat transfer in mercury is much smaller than that in water. Furthermore, it has also been reported that the fluid motion in mercury is very unstable for all values of Rayleigh number. The same phenomenon has been observed by A. G. Smirnov in Russia.(3) No rationalization has been given in both reports, however. 5. Natural Convection with Internal Heat Generation Because of the advance in nuclear power reactor technology, heat transfer by natural convection with internal heat generation has assumed increasing interest. The basic equations of this case are very similar to those without internal heat generation, with the exception that a heat source term is included in the energy equation, Theoretical solutions by Hammitt,(2) Hallman,(33) Ostrach,(34) etc., are available in the literature, among which the work of Hammitt is of most interest to the author. Hammitt studied the heat transfer in a completely closed circular cylinder with a known internal heat source (not necessarily constant), the wall temperature of the cylinder being subject to a known axial function. With the boundary layer simplification and the assumption of a large Prandtl number for the fluids, he set up the equations of motion in integral form. Hammitt further modified Lighthill's temperature and

-12 - velocity profiles(25) to satisfy his boundary conditions. The final form of these equations are then programmed for an IBM-650 computer to obtain the numerical results. The solutions thus obtained are verified by an appropriate experiment from which a good agreement has been seen. Since Hammitt's solution involves an arbitrarily known heat source and a very general boundary conditon, Lighthill's solution in the boundary layer flow regime as discussed previously would be a special case of Hammitt's solution if the heat source is assigned everywhere zero inside the cylinder except at the bottom and the wall temperature is assumed to be constant throughout. This fact has been pointed out and numerically checked by Hammitt in his paper. Later, Hammitt and Chu(35) employed the results from the previous work to obtain a numerical method for transient heat transfer by natural convection in the same geometry. No experimental data are available to verify the theoretical results from this paper, however. Also Hammitt and Brower(36) made a somewhat preliminary analysis of low Prandtl number fluid flow in the same geometry. Finally, internal heat source experiments similar to those of the above studies are presently underway in England.* To summarize the above, the phenomenon of natural convection has been quite extensively studied both theoretically and experimentally. Nevertheless, due to the inherent difficulties of mathematics and experimental technique involved, the information so far available is rather fragmental. Several elegant analyses on laminar flow have been published, * Personal communication from Dr. Do Wilki.e, Windscale and Calder Works of UKAECG September 1961.

-13 -but only under very restrictive conditions, such as a particularly simplified geometry, a limited range of Grashof or Prandtl number etc., those analytical results agree fairly well with the experimental data. With the exception of few papers such as Lighthill's(5), the analysis for natural convection in turbulent flow has virtually not been attempted. Certain emperical correlations are in existence, however, due to the absence of basic understanding, the applicability of these correlations is rather limited.

CHAPTER III THE THEORETICAL CONSIDERATIONS 1. General Considerations Figure 1 is a schematic representation of the problem configuration. The coordinate system is symmetrically placed at the closed end of the circular cavity whose length is i and radius is a. The walls of the cavity are kept at a uniform and constant temperature, To. The open end of the cavity is connected to a reservoir which is kept at a uniform constant temperature, T.o T1 is always smaller than To. By this arrangement, heat is transferred through the walls to the fluid inside the cavity and carried to the reservoir by convection and conduction. Radiation heat transfer is neglected throughout the entire analysis. The temperature difference (To-T1) is assumed to be small as compared with the absolute temperature. Furthermore, for sufficiently small temperature difference, it is reasonable to assume that flow is laminar.* The body force, due to gravity, is parallel to the axis of symmetry. As mentioned in the previous section, the fluid motion is not exclusively of one type. Both Lighthill(25) and Batchelor(22) have pointed out that three different flow regimes are possible, namely, a boundary-layer type flow regime, a similarity flow regime and a similarity flow with a stagnant portion. The parameterswhich differentiate one regime from the It appears that the magnitude of the Crashof number determines whether the flow is laminar or turbulent. Indeed, a large temperature difference will tend to make the Grashof number big, this is indeed not the only factor, however. For instance, in a centrifugal field the acceleration is so big that even under moderate temperature difference, the Grashof number can be of order of 1012. Nevertheless, the assumption that T1-To being small is essential as will be seen in following section. -14 -

-15 - x RESERVOIR, T1 it To IR Figure 1. Schematic Representation of the Problem Configuration.

other is Rayleigh's number (based on radius, see Nomenclature) and 2/ao Since each of the following solutions applies to a specific flow pattern, it is of advantage to analyze physically the significance of each flow regime, and its relation to the parameters involvedo* When the product of Ra and, a/~ is sufficiently large, the flow up the sides of the cavity approximates free convection flow up a vertical plate. This flow is of the boundary layer type; and it is well known that the solution thus obtained has given good agreement with experimental data,(8'9) For flow inside a circular cylinder it is also known that the curvature of the wall does not affect the form or the validity of the boundary layer approximation, provided that its thickness is smaller compared to the radius of curvature(1O0) The present case is different from that of the flat plate since the fluid motion occurs in a confined space which makes it necessary that the flux upward along the walls must be balanced by the flux flowing downward at the center, At one extreme, when the boundary layer occupies only a negligible portion of the cavity cross-section and consequently the downward velocity is small, the present solution should be identical to Pohlhausen's solution(8) which assumes no vertical velocity outside the boundary layer. Therefore, for a large value of Ra x a/i, the solution is of the boundary layer type. This yields the first of the three flow regimes with a downward velocity outside the layer taken into consideration. One might expect it to be valid down to a value of Ra x a/Q for which the boundary layer first fills the whole tubeo Actually, however, it should * A similar analysis is also given in Lighthill's paper (25)

-17 - break down for a smaller value of this parameter, roughly that at which there is no longer a maximum in the volume flow of unheated fluid (i.eo, of fluid outside the boundary layer) at the orifice cross-section. It is possible that unheated fluid be drawn into the boundary layer to become partially heated; but not vice-versa. When Ra x a/e is less than this critical value, the boundary layer mixes with the central flow and when a uniform condition is reached, fills the whole tube. The nature of this resulting flow can be more easily perceived at the extreme case when Ra x a/2 is sufficiently small. This is one in which all tendency for the boundary layer to thicken with distance from the closed end has disappeared (ego., similar to ordinary, fully-developed pipe flow). Then the distributions of velocity and temperature are similar at each section of the tube, only the scale increases as the orifice is approached. It will be found from the equations of motion that under these circumstances the scales of both the temperature and(axial) velocity profiles increase linearly along the tubeo The condition that the central temperature should rise from its value T, at the orifice to the value To at the bottom then determines the possible values of the aspect ratio Q/a for this particular flow as certain multiples of Rao This is the smallest value of Ra x a/e for which there is motion all along the tube. For still smaller values, which may be regarded as acheived by extending the tube beyond a hypothetical closed end, the additional length is filled with fluid at rest at the temperature of the walls, the "similarity flow" in the remainder being unchanged. This is the third regime of flow which occurs when Ra x a/e is less than the second critical value described above.

-18 - In the following sections, the solution for the first flow regime is obtained by using a modified Squire's methodo(10) An analytical solution is given for large Ra x a/e when the boundary layer is such that high order terms of boundary layer thickness in a series expansion are negligible. This solution is applicable for fluids of any Prandtl number. For the intermediate range of Rayleigh number x a/e, still beyond the first critical value, the solutions are expressed in numerical form. Again, these solutions are given in a wide range of Prandtl number (10 Pr > 002). For the similarity regime, an analytic solution is obtained only for fluids with large Prandtl number. This simplification is important in order to make the governing equations manageable. It is to be noted that an approximate solution of the same problem under this simplification has been obtained by Lighthill.(25) However, the present solution employes an entirely different method which is believed not only more rigorous mathematically, but also more meaningful physically. A factor of order of 10 separates the first critical value (corresponding to the breakdown of the boundary layer solution) from the second (corresponding to the onset of the similarity solution). The intermediate range must be filled by a solution which permits variation down the tube of velocity and temperature profiles which nevertheless fill the whole tube. Due to the great complexity present, such a solution is very difficult to obtain in a rigorous mathematical manner. However, an approximate solution(25) is known to the author (by Lighthill), and this solution is quoted in order to complete the range of present worko Indeed, it is hoped that this gap can be filled by later endeavors on this problem.

-19 - Figure 2 is a graphical representation of the similarity flow regime with and without the stagnant portion. 2. The Basic Equations The equations of motion of the present problem are similar to the ordinary equations of free convection flow. With the boundary layer approximation, the equations of conservation of mass, heat, and momentum respectively, in steady axisymmetrical flow are* U 0V V o (2.1) aX aR R aT aT,~2T 1 (T au au 1 ap 2u 1 au U - + V R = - f - - - + v(I + - (2.3) 6X OR p OX OR R OR = 0 (2.4) aR Here U and V are unknown velocity functions, U being the axial velocity (parallel to X-axis) and V being the outward radial velocity (along the radius R). T is the temperature, p the pressure which is assumed to be a function of X only as shown by Equation (2.4). This is a consequence of the boundary layer approximation. f is the external acceleration oriented in the negative direction of X. K is the thermal diffusivity, v is kinematic viscosity and p is the density, all in consistent units. X and R are the coordinates as shown in Figure 1. Figure 3 shows schematically the physical quantities. A table of nomenclature is given at the end of this paper. * The derivation of these equations can be found in the Appendix.

-20 - i 'I i I SIMILARITY REGIME WITH STAGNANT PORTION FLOW WITH BOUNDARY LAYER NOT FILLING TUBE Figure 2. Similarity Flow Regime with and without Stagnant Portion

-21 - X RESERVOIR TEMPERATURE T - T, UNIFORM AND CONSTANT T EXTERNAL ACCELERATION U T =TO ALONG THE WALLS, UNIFORM AND CONSTANT "-NR X Figure 3. Schematic of Physical Variables.

-22 - As usually assumed in natural convection flow, only density variation with respect to temperature is taken into consideration. Thus, if the temperature difference To-T1 is small compared to the absolute temperature scale, p1 varies with temperature according to the equation. p = Po [1 + a(T-To)] (2.5) where the subscript "o" signifies the conditions at the wall (R = a), and a is the volumetric thermal expansion coefficient. The condition of no slip at the solid boundary requires U = V = 0 at R = a. Thus, the momentum Equation (2.3) at the wall takes the form 0= f - + v[U + 1 (2.6) pO ax ~iO Ri(R=a2.6) Now Equation (2.5) is combined with Equation (2.6). By this means, the thermal buoyancy force is made apparent. Upon substituting dp/dx from the combined form of Equation (2.5) and Equation (2.6) into Equation (2.3) the momentum equation becomes 8U 2 U 1 fu R=R U - + V - = f(T-TO) V + - =a (2.7) 6X 6R v o) 2 R 6R R=a Equations (2.1), (2.2) and (2.7) are the three governing equations to be solved for three unknowns U, V and T under the conditions at R = a and at X = O; U = V = O and T = To (2.8) and at X = Q, R = 0; T = T1. (2.9)

-23 - It is found convenient to write the governing equations and the boundary conditions in a dimensionless form. The substitutions which reduce simultaneously Equations (2.1), (2.2), (2.7), (2.8) and (2.9) into dimensionless form with the least number of parameters are U = u V = v a T= T - t (2.10) R = ar X = ix Thus, the final form of the governing equations becomes ~u ~v v a +- + - = 0 (2.11) ax 6r r u at+ vt 2t + at (2.12) ax ar 2 r ar l(u u + v u) = -t + [2 + 1 ]r=r (2.13) 1(U V It + v - - ] rl a dx or o r +r Tr r while the conditions in Equation (2.8) become u = v = t = 0 at r = 1 and at x = 0. (2.14) and the condition at the orifice or Equation (2.9) is 4 afa (To-T1) a a R x Gr Pr = a t 2.5 tlIx=l,r=o vki le

-24 - Now, one of the most important results from any solution of the above equations will be the overall heat transfer rate. In this connection, conventional Nusselt number, based on the radius of the cavity, is defined, i.e., ha Nu = and h = Q = A(T1-To) 2Tra (T1-To) Combining these two definitions, one obtains aQ 1 1 T aa Nu = aQ _t 1 1 2 K(7-)R dX K 2ta (Tl-To) 2niK(Tl-To) Q 2jK(Tl-To) J 2 K) a dX o (T -Toy) \)R=a 0 Upon the substitution of dimensionless variables defined in Equation (2.10), the Nusselt number can be written as 1 Nu = - f ( - dx (2.16) tr r=l o Before getting into the details of solving above equations, two conclusions can be readily drawn; (1) the final results must contain the parameters Prandtl number, a, Rayleigh number, Ra, and the aspect ratio i/a, since these are the only parameters in the governing equations and the boundary conditions. (2) From Equation (2.13), it can be seen that zero velocity everywhere inside the cavity (u = o, v = o) for any nonzero t is impossible. This simply means any temperature difference across the cavity, however small, will produce some steady convective

-25 - motion.* This fact is useful in obtaining the similarity solutions in the following section. The complete derivation of the basic equations as well as a discussion of the significance of the omitted terms due to the boundary layer approximation is found in the Appendix. 3. The Solutions for the Boundary Layer Type Flow Regime (a) The Formulism of the Integral Equations In this section an approximate solution for large Ra x a/~ is obtained. As mentioned in the previous sections, this corresponds to boundary layer type flow. Because of the complexity involved in the governing Equations (2.11) through (2.15), an exact solution throughout the fluid is very difficult to obtain. Therefore, an approximate solution of their integrated forms is attempted. Experience has shown that the integral method predicts overall heat transfer rate, being an integrated quantity, with reasonable accuracy (around 3 percent at high Prandtl number and 13 percent at low Prandtl number) (14). However, the errors in skin friction and mass flow rate are usually high.(10,14,37) Since in the present case, the Nusselt number is of main interest, the use of this method to obtain the solutions is justified. The equations of motion integrated over each cross-sectional area, in dimensionless form, are 1 urdr = 0 (3.1) 0 As an example of a problem with the opposite nature, consider natural convection inside a rectangular region where two horizontal sides are maintained at different temperatures. Now the Rayleigh number must exceed a certain value in order to have non-trivial solutions for the velocity.(22)

-26 - 1 2 (3.2) a / rutdr = ()r=l (3.2) o0 2 1 2 1- a ru2dr = - rtdr + 1 (U 2u (3.3) Cra ax J 2 -r )r2 r=l o o The forms of Equations (2.11), (2.12) and (2.13) at the walls and on the axis are 2t 1 it ( + * -a) = o (3.4) Or2 ' r dr r=l at = t 1 at (u x)r=0 ( 6 + )r=o (3 5) 1 (u u a 2u 1 6u o 1 ( )= (-t + [ -2 + - ]r=l (3.6) ( ax r=o r=o r=l (3.6) For this flow regime a solution with a variable boundary layer thickness is needed. Because of this complication, a simple profile is used. No attempt is made to apply the more delicate of the above conditions, namely, those involving the values of the second derivative, as in Equations (3.4) and (3.5). Equations (3.3) and (3.6) are thus combined to eliminate (2)r= Since for this regime of flow, u is independent from r outside the boundary layer, ( ) vanishes. Thus, the momentum equadr r=o tion becomes 1 1 1 ru dr = - rtdr + 1 (1 u u + t) + (u) (3.7) a ax J2 a ax r=o 6r r=l o o The solutions of Equations (3.1), (3.2) and (3.7) together with the boundary condition Equations (2.14) and (2.15) are obtained in the form of polynomials in r. Since from the condition of axisymmetry, the solutions must be the even functions in r. That is, the terms with odd

-27 - powers of r are missing. A consideration of the physical situation as well as test results of Hammitt(2 ) reveals that the temperature function t is a monotonously decreasing function from a constant value t_ at r = 0 to zero at r = 1. The simplest form of the velocity function u, with least energy dissipation is to have one (positive) relative maximum and one (negative) relative minimum over the interval 0, r - 1, being negative constant outside the boundary layer and positive in the vicinity of r = 1, vanishing somewhere between and also at r = 1. In fact, the following profiles due to Lighthill(25) are assumed rtl 0< r < -Y 0 < r < 5 t -it)(i - (^ i p <.^(3.8) u-til ("r)2[ 0 <r<< ( u y{1 - (1 ) 2[1 + 5(r-l)]} < r < (39) 1-P where B, y, 5 are functions of x. Equations (3.8) and (3.9) already satisfy the boundary conditions at r = 0 and r = 1. Now Equations (3~l), (3o2) and (3o7) are used to determine the unknown functions I, y, and 5. It is to be remarked here that P could be physically regarded as the one's complement of the "boundary layer thickness", i.e., 1-P is the boundary layer thickness. The continuity requirement of the temperature and velocity functions suffices to determine the interval of P varies over 0 < P < 1. As will be shown later, the upper limit of 5 is completely determined by a pair of given values a and t1. The substitution of Equation (3.9) into Equation (3.1) yields 86 - 5(5+2p+p2) (3.10) (-p)22(3+2p)

-28 - while the equations of conservation of energy and conservation of momentum after elimination of 6 from Equation (3.10) give the following system of ordinary differential equations,* 1 d [Z7 -1485+1017+4722p2+2622p3+297P4+27 5] a dx 840 (3+2p)2(l1p) = 1 tl(l-P)(3+P) - Y 3(3+4)+352) (3.11) 12 (3+2P)(1-)2 d (tl 45+132+181p2+62p3) _ 1 (3.12) dx 840 3+2P 1-P Since it is assumed that the boundary layer arises from the closed end of the cavity, D = 1 at x = 0, likewise at this point y = 0. Thus x = 0, =, = 0 (3.13) (b) The Asymptotic Solutions It now remains to solve Equations (3.11) and (3o12) under the conditions defined in Equation (3o13). Two difficulties are immediately encountered, namely, (i) the given system of equations is highly non-linear, (ii) the given condition at x = 0 happens to be the singular point of the equationso An analytical method of solving a system of non-linear equations is not generally known, and it seems very unlikely that Equations (3.11) and (3o12) can be solved analytically. Therefore, a numerical method is relied upon for the solutions of the above system. However, a * A complete derivation is given in the Appendix.

-29 - well defined initial point is necessary to initiate the numerical process,* In the present case, the conditions known at x = 0 do not define the first derivatives of B and y and consequently the numerical integration cannot be readily performed. Since the system of Equations (3.11), (3.12) and (3-13) do not satisfy Lipschitz condition,** the solution may exist but uniqueness is not guaranteed. This simply means that through the given point x = 0 there are infinitely many solutions for the given system of equations. Therefore, it is necessary either to successfully remove the singular point from the differential equations or to determine from physical reasoning which one of these infinitely many solutions would yield the desired results and thus to proceed immerical process. It has been found in the present case that an asymptotic integration of the above equations combined with a conclusion drawn from dimensional analysis will determine the desired solution with success. Substituting P = 1 - 5 (3.14) into the above differential system, the following equations can be obtained 1 d (2f) = til - 7I2 (3.15) a dx d (7g) = 1 (3.16) dx * See Milne, "Numerical Integration of Differential Equations." McGraw-Hill Co. See Ince, "Ordinary Differential Equations." Dover Publication, p. 23.

-30 - f = 240-65505+4880~2-13603+14445-956 6 1 1 pf 280(5-22)2 f' df -3000+3600o-450t2-960 3+540 4-117t5+9t6, 1 p df = Pf dS 70(5-2t)3 t2 t2 420-680o+367t2-62 3 84o(5-2) g, dg _-2560+3670t-1664 2+248g3 dt 840(5-2t)2 I1 = 1__ (4-) = =Pi1 12 Z 3(10-10t+352) 1 1 j (5-20) t2 I2 2 and the initial conditions become x = 0, = 0, g = 0 (3.18) When the Rayleigh number x a/i is large, the boundary layer thickness,, is small and the higher order terms of g are neglegible as compared with the leading coefficients. Thus, the polynomials in Equation (3.17) can be approximated by appropriate constants. The differential system Equation (3o15) and Equation (3.16) is then reduced to 1 d (12 ) = t - 6 - (3.19) a dx 35 3 t2 d (1 7) = (3.20) dx 10

-31 - By expanding the differentials in above equations, Equations (3.19) and (3.20) become 12 1 (2 7 _ d) tl - 6_ (3.21) 35 a S dx _2 dx 3 2 1 = Z (3.22) 10 dx from which the following equations are readily obtained. d _ -35atl 3+30y(21a+24) (3.23) dx 3672 dy = 10 (3.24) dx S Now it can be seen that the independent variable x is missing in both equations. As a result of that, an elimination of dx yields a first order equation with variables y and 5 only, i.e., dy -36072 (3.25) d7 [35ct1i3-307(21l+24)] By a transformation of 7 = |3g (3.26) Equation (3.25) becomes d = -105atl@+1800 02+1890ag dS ~ 'h atl-30 0(21a+24) The variables are now separable and the result of this separation is dS = 35atl-30 Q(21a+24) d (3.27) 0Q[-105atl+90 8(20+21a)]

The general solution of Equation (3.27) is inC + ent = - 1 nOG - 120 n[-105atl+90 0(20+21a)] 3 90(20+21a) or 1 4 C~ = [-105at +90 Q(20+21a)]3 and by Equation (3.26) 4 Cy = [-105atl+90 y2 (20+21a)] 20 (3.28) Equation (3.28) clearly shows that the given condition 0 = 0, y = 0 does not suffice to determine the unknown constant C uniquely as predicted previously. In order to determine a particular solution for Equation (3.27) which is physically meaningful as well, the results from a dimensional analysis are invoked. It is well known that in case of boundary layer type natural convection flow, the boundary layer thickness grows according to the fourth root of the distance x, (e xl1/4) while the increase of velocity, 7, is proportionate to x3/4.(1314) Therefore, the ratio of ~3 to 7 is independent of the distance Xo Bearing this in mind, one can easily show that Equation (3.28) does yield this particular solution if one sets C equal to zero. Furthermore, the proportionality constant can also be determined from the same equation, i.e., 7atl 3 (39) 120(1+1.05a) with this explicit relation of 7 and i, the system of Equations (3.23)

-33 - and (3.24) can easily be solved. The solutions are found to be* =3 88(1+105')1/4 x1/4 (3.30) atl = 3.43( tl )1/4 x3/4 (3.31) l.05a The Nusselt number can be computed by combining Equations (3.30) and (358), i.e., 1 1 Nu= f ( t) dx = 2 f dx (3.32) tlt0 r r=l o Thus, the overall average Nusselt number, as from Equations (3.30) and (3.32) is ctl 1/4 Nu = 0.692( t l 4 (3.33) By definition of "asymptotic integration", the accuracy of the results increases as 5-1 increases and approaches exactness when 5 is infinitesimally small. This is equivalent to say that the accuracy increases for an increase in t1. It is therefore natural to investigate the lower bound of t1 below which the error introduced becomes intolerable. Owing to the non-linear nature of the problem, it is almost impossible to compute the error as a function of tl in a rigorous mathematical * A less elegant but /asier way to obtain these solutions is to assume the solutions | = Axl/, 7 = Bx'/4 and determine the constants A and B from Equations (3.23) and (3.24). Indeed, if the postulated solutions are consistent, they should reduce these differential equations into algebraic equations for which A and B are unknowns. However, in doing this, the insight into the differential equations is completely lost.

-34 - manner. Nevertheless, it is not too difficult to estimate the order of magnitude of the error by comparison with a known "exact solution" and thus to determine a lower bound, If the Prandtl number of the fluid is sufficiently large, the left side of Equation (3.11) vanishes and the equation is reduced to an algebraic form. With this simplification, Lighthill(25) obtained a set of solutions by integrating numerically the energy equation, which, in the present notation, can be written as t1 d [3(4- )(42o-680o+367=2-62 3)] (3 dx (10-10~+32) His solutions are carried out within the range of 0.62 > g > 0.1 which corresponds to 5400 < t1 < 2,4 x 106. The lower limit of t1 in this case is determined by the condition of maximum volumetric flow at the orifice which has been discussed in the previous sections. If a - co, the asymptotic solutions Equations ( 3.30), (3.31) and (3.33) approach the limiting values of = (240)1/4 x1/4 (335) tl For the purpose of comparison, the range of the solutions of For the purpose of comparison, the range of the solutions of Equation (3.34) has been extended to < 0.01 which corresponds to t1 > 4.0 x 1011. A trapezoid method is used for this integration and the numerical computation is done by an IBM-704 computer.

-35 - The exact solutions for a -> from Equation (3.34) compare with the asymptotic solutions from Equation (3.37) remarkably well. The deviation is less than 1.2% when t1 is 2.64 x 106 which corresponds to a boundary layer thickness of 0.2, about 5.42% at t1 = 4.04 x 104 or at the boundary layer thickness 0 = 0.3. Figure 4 shows the comparison over the region investigated. The dimensionless volumetric flow rate of the cold fluid at the orifice of the cavity is proportional to the flow rate of the core and the cross-sectional area of the orificeo That is, G- 7(1-)2 = y72 (3.38) where G is the dimensionless flow rate. It has been found* that the right side of Equation (3~38) can be expressed by p3(3+) (3+2) (1-p)3 36(3+4p+3p2) and this quantity ceases to be a maximum when B = 0.38 or 0 = 0.62, correspondingly t1 = 3,400, which is the lower limit for this regime of flow. Substituting in Equation (3.29) with a ->0o, and Equation (3.35) into Equation (3o38), it is obtained that G (1-_)2 1 tl3 (3.39) Equation (3-39) has a maximum at 0 = 0.60 which indicates even for relatively large i, the asymptotic solutions are still reasonably A brief review of Lighthill's solution appears in the Appendix.

POINTS TO THE LET OF I POINTS TO THE LEFT OF THE CROSSES INVOLVE A PHYSICAL IMPOSSIBILITY ASYMPTOTIC SOLUTION ---- "EXACT SOLUTION" I I I I I I I I I I I I I I I I I. f I I I I I I I I ' 1 1 1 I I I I I I I I I - I I I I 1 -' II, I - 1. I I - I I II I - i.0 I a: w _J I Z U) U) z J I w 0 n 9 10 -I ____ ____-__ —^^^ ^1 --- - _-_ _-_ ---........ I I w I>rQo zS ~ - 3 a ~- 0 o=2 IC )2 103 105 106 RAYLEIGH NUMBER x 0/I Figure 4. Comparison of the Asymptotic Solution to the Exact Solution for Large Prandtl Number.

-37 - applicable. However, the corresponding t1 computed from Equation (3o35) is 1,850 which deviates considerably from the actual value of t1 = 3,400 (47.1% error). It is,therefore, concluded that when the Prandtl number is large, asymptotic solutions are good for tl > 104 and the deviation becomes large as tl decreases from this value. The critical value of t1 as predicted by Equation (3.39) is not good. Nevertheless, the critical boundary layer thickness at the orifice is in excellent agreement with that obtained by exact integration. There is no absolute guarantee that this conclusion can be equally well applied to those solutions of finite Prandtl number. However, at least, the comparisons above may serve as a possible guide as to how the solutions vary with respect to the physical parameters. It is indeed true that the exact solutions of Equations (3o15) and (3.16) still rely on numerical integration which is the subject of the next section. (c) The Solutions in Numerical Form In this section, a method is developed to obtain numerical solutions of Equations (3o15), (3.16) and (3.17). An IBM-704 digital computer is used for this purpose. The detailed computer program is included in the Appendix along with the flow diagram, The difficulty due to the presence of a singularity at x = 0 can now be overcome with the aid of the known asymptotic solutions. It has been shown that these solutions are particularly good when ~ is small, and it is also physically obvious that the boundary layer thickness g is a continuous function of x throughout the entire region of investigationo Thus, in the neighborhood of x = 0, regardless of the value of

-38 - tl, 5 must be small as required by this continuity (since at x = 0, l = 0), and consequently the relationship between i, 7, and x is accurately described by Equations (3o30) and (3.31) if x is chosen sufficiently close to the origin. The differential system Equation (3o15) and Equation (3o16) is everywhere continuous except at the origin, x = 0. Therefore, the numerical integration can be initiated in the neighborhood of x = 0 with the values of ~ and y predicted by Equations (3.30) and (3.31). Indeed, high accuracy can be obtained by first taking both a sufficiently small value of x and then a small interval of numerical integration. This is feasible since the tedious repetitive computation is done by a computer, Equation (3.15) and Equation ( 3,16) could be programmed directly for mechanical computation. A sufficiently small initial value of x could be taken and the corresponding values of ~ and 7 computed from the asymptotic solutions Equations (3.30) and (3.31)o The values of the first derivatives d~/dx and dy/dx are then to be computed from Equation (3.15) and Equation ( 3.16) respectively. From this point on, Runge-Kutta or any known method for numerical integration of differential equation can readily be applied, and the differential Equations (3.15) and (3.16) pointwise solved. In the present case, however, direct programming has proven disadvantageous. The independent variable x is missing both in Equations (3.15) and (3.16). If the iteration is on x, the stability of machine computation is greatly reduced. As a matter of fact, in some instance, the error is so big that overflow of the computer was detected

-39 - after only the first few computations. Thus, an amended form of Equations (3.15) and (3.16) is needed. Expanding the differentials in Equations (3.15) and (3.16), one obtains 72f' d_ + 2yf d = a(tlIl-I2) (3.40) dx dx yg, d_ + g d_ = 1 (3.41) dx dx 5 where the coefficients are the same as those defined in Equation (3.17). It is found convenient to write these two equations as follows: Pf, -d + 2 Y Pf dy = aI(tlPI1 - L- PI2) 2 dx dx 52 (3.42) f~ d7 1 7g' d + g d- - (341) dx dx The new coefficients have also been defined in Equation (3.17). From Equations (3.41) and (3.42), d~/dx and dy/dx can be explicitly solved, i.e., d_ = tlgPIl 3 - 7(agPI2+2Pf) dx 2 (gPf -2SPfg') (3.43) dy _= tlg 'PI14 - y(og'PI +Pf) (44) dx y7(2|Pfg'-gPf,)

Dividing Equation (3.44 ) by Equation (3.43), one obtains dy = atlg'P iPi4 - y(ggPi +Pf,) _ d atigPI 3 - 7 (agPI+2Pf) (3.45) One also obtains from Equation (3.43) dx _ 72(gPf t-2Pfg' ) d atlgP1S3 - 7(cgPI2-2Pf) (3.46) and the Nusselt Number can be expressed in terms of (, i.e., from Equation (3.32) 1 11 Nu = 2 dx =2 f e d (3~47) o 0 o o a where _1 is the boundary layer thickness at x = 1. Equations (3.45), (3.46) and (3o47) are now programmed for machine computation. The program is so written that the input data are Prandtl number, a; the product of Rayleigh number and a/i, tl; and the initial point, Po (sufficiently small). Corresponding to these values of P o, a t, the computer starts to compute xo and 0o from the asymptotic solutions. With all these as initial values, the computer proceeds to integrate Equations (3,45) and (3.46) by the Runge-Kutta fifth order method.* The pointwise results are stored in a temporary space for the integral of Equation (3.47). The accuracy of integration of Equations (3.45) and (3.46) is of order of At to the fifth power, while for Equation (3.47) the accuracy is only of second order. The machine prints The equations involved are included in the computer program shown in the Appendix.

-41 - out the values of x, 7 and 7(1-_)2 at each ~. It is terminated when the condition x = 1 is satisfied. If x = 1 happens to be somewhere in between the intervals, the size of At at the last point is adjusted according to a quarterth power extrapolation. During the numerical integration, the values of 2 1 dx At for each 5 are computed and summed. The g dS final total of this quantity is printed out at the end of the program. This is the overall Nusselt number computed under the given conditions. It can be observed that both the denominator and the numerator of Equation (3.45) are the difference of two terms. A simplified numerical computation shows that their numerical values are so close that their difference begins to show only after the fourth or fifth digit. This is particularly true when the Prandtl number is large; -- e.go, accuracy of 12 digits is required when a = 104. Furthermore, the stability of the solution decreases as the Prandtl number increases. This is evidenced by the rapid build-up of the errors incurred in each step. Once an error is introduced, the numerical solution of y converges to a small constant value very close to the initial value of r corresponding to xo, This result is not at all surprising if one examines the orders of magnitude of each term in Equation (3.45). The numerator of Equation (3.45) is of 0(&4) while the denominator is of 0(53) and 5 is of 0(0) in the neighborhood of x = O. Therefore, if an error, c, is introduced in the numerator during the numerical computation, the error simultaneously introduced in the denominator would be of 0(-) which is much larger than Eo is than consequently of 0(t) or 0(0). The same argument is repeated when the next point is computed except that now an error has already been inherited in y since the numerical integration

-42 - employes a slope approximation (or a truncated Taylor Series approximation), ie, Yn+l = Yn + (RAS) n (3.46) where R is a correction of Runge-Kutta method and n = (d o Therefore, it can be concluded that if the required accuracy is not met, the error would build up rapidly, and the solution of y becomes approximately constant depending on the initial value of xoo A similar conclusion can be drawn from the fact that Equation (3.45), being a first order equation, possesses a singular solution of y = 0 (which is physically meaningless). This solution happens to be stable as far as numerical integration is concerned. Therefore, it is only natural that an additional accuracy is required in order to obtain a solution along a relatively unstable patho It is, however, fortunate that this obstacle is not overwhelming. It is well known that if the Prandtl number is big, the overall heat transfer rate is not much affected by the magnitude of this parameter.(8,9,10) This is further evidenced, in the present case, by the asymptotic solution Equation (3-33) that a is hardly important at all if it is big enough. In fact, Lighthill's(25) solution of a - oo was well verified by experiments on water(29) having a maximum Prandtl number less than 10. As pointed out previously, the numerical integration employes a Runge-Kutta fifth order method, and the maximum step size does not exceed 0.01 for all the computations. Therefore, the truncation error should be less than 10-10 for each step of computation. The accumulated error and the stability criterion are very difficult to estimate.* No attempt has been made in this connection. See Milne, ibid.

-43 - Altogether about 30 machine computations have been made(the results are listed in Table I). The range of Prandtl number is from 10.0 to 0.02 and of Rayleigh number is from 108 down to 103 for fluids of larger Prandtl number and from 1010 down to 105 for fluids of low Prandtl number. In each set of computation there is a check of volumetric flow rate at each axial distance x. This device was intended to determine the lower limit of tl; i.e., that condition when a maximum volumetric flow rate would not be located at x = 1. Regardless of the parameters t1 and a, this maximum always shows at the orifice. But the Nusselt number 1/4 becomes more and more independent of tl (NuN% tl for high tl). The failure to obtain such a minimum t1 below which this analysis does not apply is thought to be due to the error introduced in the pointwise integration. However, there is no rigid basis for this argument. For the present cases, the minimum t1 for each a is then taken when ~ > 0.6 at x = 1 as per Equation (3.39). This criterion should give a fair estimate as explained in foregoing section. (d) General Remarks on Computer Program At least, three different methods of programming have been tried. A Runge-Kutta fourth order method has been used occasionally.* The results from this program were compared with the results of Runge-Kutta fifth order method and no significant difference observed. A double precision fifth order program has also been used for integrating equations with high Prandtl number, but no success obtained. On the average, a complete set * See University of Michigan Computing Center, MESS.

-44 - TABLE I SUMMARY OF THE COMPUTED RESULTS FOR BOUNDARY LAYER TYPE FLOW REGIME 1. a = 0.02 3. a= 1.0 ti 5 x 1 x 5 x 8.5 x 1 x 5 x 1 x 1 x 103 4 104 4 10 105 105 6 106 107 Nu 3.7161 3.8311 4.5014 4.8079 4.9410 6.8062 8.0149 14.1259 ti 5 x 103 1x 104 1 x 105 1 x 106 1 x 107 1 x 108 Nu 4.7758 5.5616 9.8277 17.7087 31.7950 56.8675 2. a = 0.1 4. a = 10.0 ti 1 x 103 5 x 103 1x 104 2 x 104 5x 104 1 x 105 1 x 106 1 x 107 1 x 108 Nu 3.6606 4.o806 4.3461 4.9077 5.9621 6.6381 11.6010 20.7800 37.2118 tl Nu 1 x 103 5 x 105 1x 104 1 x 105 1 x 106 3.3396 5.0817 6.1296 11.2995 20.9869

-45 - of computations takes approximately two minutes on an IBM-704 type computer (using Runge-Kutta fifth order), while double precision program requires five to six minutes. Figure 5 is a consolidated representation of Nusselt number as a function of t1 for different Prandtl numbers. The dotted lines are the trace of the asymptotic solution. It can be seen that the deviation increases with decreasing tl. The numerical solution reasonably follows the "one quarter power" law and flattens out when a certain value of t is reached. The mark on each line is the estimated value of t1 corresponding to the break-down of the model. This is given by the criterion of > 0.6 at the orifice. Figures 6, 7, 8 and 9 are the boundary layer thickness and the axial velocity as functions of x. Again the asymptotic solution is shown for comparison. Figures 10 through 15 show the temperature and velocity distributions. 4. Solution for Flow Regime With Similarity (a) The Similarity Transformation In this section a solution of Equations (2.11), (2.12) and (2.13) under the boundary conditions Equations (2.14) and (2.15) is obtained, such that u and t are proportional always to the same two functions of r, the factors of proportionality varying with x. A brief examination shows that the only possible form of variation with x is that u and t are both proportional to x, measured from some origin. In fact, such a transformation linear in x does exist as shown following:

10, z 10 1.0 I I - - - - ASYMPTOTIC SOLUTION NUMERICAL SOLUTION POINTS TO THE LEFT OF 0 MAY INVOLVE A PHYSICAL IMPOSSIBILITY 8 10 t =Ra x a t Figure 5. Consolidated Solutions in Boundary Layer Type Flow Regime.

-47 - 1.0.9.8.7 x w C-) z C) 0 -J x.6.5 tF lRa =10 ' =1.0 a'.I/.o =.02 i.4.3.2.1 L / / fA I n L 0.1.2.3.4.5 BOUNDARY LAYER THICKNESS, C Figure 6. The Development of Boundary Layer (t = 1 x 105)

-48 - 1.0.9.8 7.6 x w 0 z.5 X C, x.4.3.2.1.1.2.3.4.5 BOUNDARY LAYER THICKNESS, C Figure 7. The Development of Boundary Layer (t1 = 107).

1.0 x U z -J 4 L.. t1_ R0 -- ____ 6,, ~-7-~~~' =/ // //, v L I//,/ ------ ^ 7 / ---- — x________ ____ __ __ _ _/I -- i........ ______ ___ --- -- - -- - -— _ il- z — _ _ __ ===='Z//......==.I I -I\0 I.01L I 10 100 VELOCITY COEFFICIENT, y Figure 8. Velocity Coefficient, 7, as Functions of Axial Distance (t1 = 105)

1.0 x a4 U) 4 x 4.01 -__________ ____/'/ /! 10 100 200 - - - / / _ n~~~~~~~z z zz zz~~~ _____ _/ / / _ _ _ _ ___ ___i#2~ I o I VELOCITY COEFFICIENT, Y Figure 9. Velocity Coefficient,y, as Functions of Axial Distance (tl= 107).

-51-,-, I. /, /,,10 X 0.2 X=0.5 /; t 3.0 X 0. - X 0.8 0.6 W 1 2.0 / -X00.4 W IF) gure X =o.5, y:26.20 U)s) U)z ~ ~ ~ ~~Ai W w -.0 0 0.9 DISTANCE FROM CENTERLINE Figure 10. Velocity and Temperature Profiles at Various Axial Positions.

-52 - 7 6 5 N I0 -J w -i 4% 4z U) U) 4 z 0 U) z w 0 El 4 3 2 1.0 w.8 w 4 w.6 U) UJ z w z.4.2 0 -I -2 Figure 11. Velocity and Temperature Profiles at Various Axial Positions.

-55 - C- =0.1 t, - i0 5.0. 4.0 Io C3 w 3.0 x C) C) w 2.0 0 z w 2 1.0 0 X=0.2_ ____X___ / /=0.5 x =0,8 / X =0O.2, Y = 0.1405 x 10 -/000~ -0 X =05, r 0.3028 x 02 // _/- ' —X =0.8, Y=0.4665xl0 - At // //,/'/~ 1.0 0.8 -- w r 0.6 g a. w z IC,) 0 C, 0.2 z 0E I.01 1.0 0.9 0.8 0.7 0.6 DISTANCE FROM CENTER LINE 0.5 0.4 Figure 12. Velocity and Temperature Profiles at Various Axial Positions.

-54 - x =.2, = 34.54 a = 0. I: =10 7 7 x:.5, Y=71.50 II 0.W -J I w W.4.8 () u,,I/ — I Cn X=.8, 101.5. oZ =.8 w z W w I -iiim.4 0 _____0 1.0.9.8.7.6.5.4 DISTANCE FROM CENTERLINE Figure 13. Velocity and Temperature Profiles at Various Axial Positions.

-55 - ':1.0 tl 10" >: 10 -.9.8 DISTANCE FROM CENTERLINE Figure 14. Velocity and Temperature Profiles at Various Axial Positions. C) > 2 5 -X:0.2, Y:O.868x - / I -x: 0.5, X:0.3952 xl2 -1.0.9.8 DISTANCE FROM CENTERLINE Figure 14. Velocity and Temperature Profiles at Various Axial Positions.

DIMENSIONLESS AXIAL VELOCITY, u/Y b ~,~ coN 0o 0 (0 " — ^ (D 00 Po.P Hi - PH 0 (t C) 0) Z P m H( C OH 0 nl DIMENSIONLESS TEMPERATURE, t /t,

-57 - By Equation (2.11), a stream function * exists, such that U = ai r r Or v = 1 _ r ax (4.1) The transformation which possesses the "similarity characteristics" is found to be * = xf(r), t = xg(r) (4.2) It should be noted that the functions f and g in Equation (4.2) are not to be confused with those mentioned in the previous section (which will not appear again). After this transformation, Equations (4.1), (2.12) and (2.13) become u = x f'(r), r v = - 1 f(r) r (4,3) f'g - fg' = rg" + g' 1( ft2 - f +g + 1 1f t 1 + C f r2 r r r where C = -f"'(l) + f"(l) - f'(l) (4.4) (4.6) and the prime denotes differentiation with respect to r. The boundary conditions are transformed into f(l) = f'(l) = g(l) = 0 (4.7) f(0) = 0, g(0) = t1 = constant (4.8) and moreover, from symmetry, the solutions must be even functions of r.

-58 - Indeed, Equation (4.5) is of third order and (4,4) of second order, the five conditions in Equations (4o7) and (4.8) would suffice to determine a particular solution for this differential system had the unknown constant C been absent. Owing to this additional complication, one more physical condition is explicitly utilized; i.e., the velocity u must be bounded everywhere in the region of r, that is to say, | im f'(r)l < M (4o9) r 0 O r where M is a finite constant. The transformed energy and momentum Equations (4.4) and (4.5) respectively are non-linear and of high order. It does not appear that an analytic solution can be obtained without simplification. This situation is further obstructed by the given conditions in Equations (407), (4.8) and (4.9), so that even a numerical solution is very difficult to obtain. As stated before, a well defined initial point is necessary for any type of numerical process. For the present differential system, three conditions are given at r = 0 and another three are given at r = 1o Therefore, if a numerical solution is attempted,three additional conditions must be guessed in order to initiate the numerical process and then the results must be matched to the conditions at the other end. The condition specified by Equation (4.9) is most unsuitable for any numerical computation since it is almost impossible to obtain a limiting value from mechanical calculations. Equations (4.4) and (4.5) are subject to a set of homogeneous boundary conditions. Since the differential system is non-linear, the

-59 - well-developed theory of homogeneous ordinary differential equation and eigenvalue solutions are not quite applicable. However, it would not be surprising that this system of equations would exhibit a behavior similar to that of linear equation (which happens quite often in mathematical physics, such as in some non-linear vibration cases). This is to say, the solutions may exist only under certain values of the parameters a and tj. It will be shown later this is the case. The given system will possess solutions only at discrete values of t1 for a given a. These values of t1 are accordingly the eigenvalues. Since the differential equations are non-linear, the sum of the solutions is no longer a solution. Further consideration will be required to decide which one of these eigenfunctions is the desired solution. In summary, the given differential system Equations (4.4) through (4.9) is inherently difficult to solve by any means. The analytical method is practically unknown for a non-linear system of this kind, The numerical type solution encounters three serious difficulties; i.e.: (i) the "two-point problem" with three unknown quantities at either end. (ii) the convergence requirement at r = 0 and (iii) the proper choice of t1 for each a (such that tj happens to be an eigenvalue of the system) before the initiation of the numerical process. Such a choice of t1 must be completely blind. It must, therefore, be hoped to solve these equations either with the aid of certain simplifications or under very restrictive conditions. This latter approach is chosen. It will be seen in the

-60 - natural convection, these governing equations are coupled together. It has been found fruitful to solve Equation (4.11) in terms of g and then fully utilize this solution of f to determine the unknown function g from Equation (4.10). (The elimination of one variable from these two equations is absolutely hopeless.) Equation (4.11) is linear in f and this fact is particularly useful since its complementary solution is readily obtainable. In fact, the homogeneous part of Equation (4.11) is r2f"' - rf" + f' = 0 (4.12) which is of Euler's form and whose solution is f = c1 + (c2 + c3 lnr) r2 (4.13) From Equation (4.13), the Green's function subject to the conditions of f, i.e., f(l) = f'(l) = f(0) = 0 (4.14) is constructed, and finally the particular solution of Equation (4.11) is obtained in terms of g. The Green's function method to solve a nonhomogeneous equation is well known,* The particular application to this problem is given below with further details in the Appendix. The Green's function of Equation (4.11) is obtained from its complementary solution Equation (4.13) and the given conditions Equation (4.14). This function satisfies the homogeneous Equation (4.12) and the end conditions Equation (4.14); it is continuous throughout the whole interval 0 r l 1 up to its second derivative where a finite discontinuity * See for instance, Ince, "Ordinary Differential Equation" Dover Press.

-61 - following sections that a solution restricted to large Prandtl number (a - oo to be precise) is obtained. It has been demonstrated experimentally that the solutions corresponding to large Prandtl number are applicable, to a reasonable accuracy, to fluids even with Prandtl number of order of one.(2,25) With this restriction the simplified equations possess an analytical solution which illustrates all the particular aspects of the original Equations (4.4) and (4.5) that have been described previously. When a -~o, Equations (4.4) and (4.5) are reduced to following forms; rg" + (f+l)g' - f'g = 0 (4.10) r2f"' - rf" + f' = (g-c)r3 (4.11) Comparing Equations (4.10) and (4.11) with Equations (4.4) and (4.5), it can be seen that the order of both equations is retained. The nonhomogeneous term (g-c)r3in Equation (4.5) appears again in Equation (4.11) and finally, though the non-linear terms are somewhat simplified, the basic non-linear nature is still preserved [see Equation (4.10)]. The analytical solution of Equations (4.10) and (4.11) will appear in the following section. Since a truncation of a power series is involved, the solution thus obtained is checked by a numerical calculation from an analog computer, with the initial information supplied by the analytical solution. An altered form of Equations (4.10) and (4.11) is thus needed for suitable analog computer programming, the details of which will be discussed later. (b) The Eigenvalues and the Approximate Solution The task now is to solve Equations (4.10) and (4.11) subject to the conditions in Equations (4.7) and (4.8) and (4.9). As is typical of

-62 - of - 1 is allowed. With these specifications, the desired Green's r2 function K(r;|) is found to be K(r;S) {1-2 + 2,n t (l-_2)n } r2 0 < r < -= {1 + (2;,n r-l)r2} < r < 1 (4.15) Once the corresponding Green's function is known, the solution of Equation (4.11) is readily obtainable by integrating the Green function and the nonhomogeneous term (g-c)r3, i.e., r f(r) = {l+(2en r-l)r2}{g(~)-c} d 4o + {1-_2 + 2 -n _ (1-2)n r f{g() - C} rA 3d (4.16) The integration of Equation (4.16) yields r 1 1 f(r) = + 2r2n r f gr3dS-r2 f gt3dC 1 11 + f r2gedS + S 2r2g $ned~ - f 2r25n rgid,} r r r 1+ Cr2(l+2n r-r2) (4.17) where g = g() Equation (4.3) and Equation (4.9) are now used to determine the unknown constant C. In this connection, Liebni-tz rule of differentiation with respect to a parameter is required; i.e., for continuous functions F(S;r), a(r), P(r) (where 5 in F is a parameter), the following theory is true.* See Kaplan, "Advanced Calculus," McGraw Hill Book Co.

-635 - a(r) f F(;r)dE = F[3(r),r1'(r) - Fiac r),ric i(r) ar c(r) B(r) + - F(f;r)dS a(r) ar By differentiation of Equation (4o17) with respect to r, one obtains the velocity u from Equation (4o7), namely, 1 1 u - dr x {n r J g(5)53d5 + fJ g(S)>ni dS r dr o r 1 - sn r f g()SdS} + ~ C(l + Qn r-r2)x (4.18) r when r approaches to zero, 1 r 1 u(O,x) = x lim, {in r[ f g33 + f gdl - f g~d+ + -C] r ->oo o o 4 1 r 2 2 + Singgdg - f intgd + I C - T Cr2 Thus, for a convergent u as required by Equation (4.9) the unbounded terms of the above equation at r = 0 must be cancelled out among one another. By equating these terms to zero, the unknown constant C is found to be 1 1 2 C = 4 f gadS - 4 J g3-dS (4.19) o o Thus, the final solution of f, satisfying all the given conditions is obtained. f(r) d + rr J gd + 2r ingd f(r) = L o g 3d| + 2r2Cn r f gd5 + r2 g^d + 2r2 f ngd~ 4f g o r r 1 1 1 1 - 2r2 f gd3d + r4 J g~3d + r2 f gid - r4 f gidj} (4.20) Now, Equation (4.20) can be combined with Equation (4.4) to form a differential-integral equation. However, if this is done, the final equation contains very complicated integral and differential terms,

-64 - and the present limited theory on integral equations is not enough to tackle any problem of this nature. Thus, in order to obtain an explicit solution for f(r) and g(r), Equation (4.20) should be further modified. Assuming the function g(r) can be expressed in a power series, i.e., 00oo g(r) = ann (4.21) the integrals in Equation (4o20) can nrw be carried out. In fact, upon the substitution of Equation (4.21), Equation (4.19) becomes* f(r) = o ( n+2)2(n+4) [(n+2)(1-r2) - 2(1-rn+2)] (4.22) 2 n=o (n+2)2(n+4) and Equation (4.4) is now of the form E n(n-l)anrn- + E nanrnn=o n=o + no nanrnl}{J ( j [(j+2)(l-r2) - 2(1-rJ+2)]} 2 n=o J-0 (j+2)2(j+4) - anrn}{kL (k+2)2 [(k+2)(1-2r2) -(2-(k+4)rk2) n ~o k-o (k+2)2(k+4) 0 (4.23) Before the solution of Equation (4.23) is attempted, several interesting facts can already be seen from it. First, the roots of the indicial equation of Equation (4.23) are 0 and -1 so the difference between them is an integer. This means one of the solutions of Equation (4.4) will possess a logarithmic coefficient. This solution thus is not * The detail of this solution can be seen in the Appendix.

-65 - convergent at r = 0 and is discarded. The second term in the series expansion gives al = 0, and the recruiting formula requires all the coefficients of odd powers of r to be zero. Therefore, the solutions of Equation (4.23) convergent at r = 0 is a series with only even powers of r. This is exactly what is required by the assumption of axisynnmetrical geometry. Indeed, it is clearly shown by Equation (4.22) that if g(r) is an even function, f(r) is also an even function. Therefore, a consistent set of solutions can be obtained after the coefficients in Equation (4.23) are determined. Since all the coefficients of odd powers of r are missing in the series solution, Equations (4.21) and (4.22) can be written as 00 g(r) = Z a2nr2n (4.24) n=o 1 E a2nr 2) 2n+2)] f(r) = - z a2n [(n+l)(l-r2)-(l-r2n+2)] (4.25) 8 n —=o (n+l)2(n+2) Thus, the solutions of the energy Equation (4.10) and the momentum Equation (4,11) with these forms of power series will be free from singularities at r = 0, Now Equation (4.24) and Equation (4.25) are substituted into Equation (4.4) to determine all the coefficients such that the given conditions g(0) = t1 and g(l) = 0 (4.26) are satisfied. [Note that Equation (4.25) will satisfy all the conditions of f(r) and the momentum Equation (4.3) provided that Equation (4.26) is satisfied. The verification can be found in the Appendix.]

-66 - From the first condition in Equation (4.26), the leading coefficient in the series solution of g(r) must be tlo By the second condition in Equation (4.26), the form of g(r) at r = 1 is a0 + a2 + a4 + a6 +... =0 (4.27) Since a0 is non-vanishing, the eventual determination of all the coefficients of g(r) will be expressed in terms of aO (the recurrence formula is not necessarily explicit and in the present problem, it turns out that the relations between these coefficients are extremely complicated). Therefore, it can be seen from Equation (4027) that only certain discrete values of ao, and hence the values for tl, will yield a compatible solution for the differential system Equation (403) and (4.4), namely, those values of a0 which satisfy the algebraic relation Equation (4.27). In principle, there are infinitely many roots for Equation (4.27), but as will be shown later, the number of physically meaningful solutions is very limited. The substitution of Equations (4.24) and (4.25) into Equation (4.4) yields the following relation: 16 E n2a2nr2nl = E j [jr2-2(j+l)r4(j+2)r2J+4] n=o j (j+l)2(j+2) [ E a2n -1] a2j [jr2-(j+l)r4+r2j+4] n=o j=o (j+l)2(j+2) [ na2nr2n-li] 2 Recalling Cauchy's formula for the product of two infinite series* Z an. Z bn = E Z anbn-j no n=o n=o j=o (4.29) See, for instance, Rainville, "Intermediate Differential Equations."

-67 - and shifting the indices, Equation (4.28) becomes 16 7 na2r2n-l = a2(n-1){Z f ja2j r2n-1 n=o n=l j=o (j+l)2(j+2) 00 200 00a2j -Z (n-1)a2(n-l){Z ja }r2n-i - a ) 2a a2 r2n1 a2(n-2){ J= (j+)l)2(j+2) n2 0 (+l)(+2) co co 00 n-2 a ja2(n, 2n+ n (n-2)a2(n-2){o Z aj +}r2n-l + z a2 2(n 2) r2n-l n 22 j2 o (j+l)(j+2 n=2 j=o (n-j-l)2 00 n-2 - E a2jr-a2(n-j-2) r2n (4-30) n=2 j=o (n-j-1)2(n-j) The difficulty in determination of the coefficients in Equation (4.30) is from the fact that this equation is not in the form of a pure power series (due to the nonhomogeneous terms in the momentum equation). The coefficients are related to one another in a very complicated fashion. Thus, an explicit recurrent formula is not available. In order to determine the coefficients for the leading terms in a manageable way, a truncation of the infinite series must be made (up to the sixth power of r). It has been found convenient to extract the leading term from the series in Equation (4.24); i.e., if bo = ao, g(r) = bo(l + b2r2 + br + br6 + ) (4.31) and after collecting the coefficients of various power of r in Equation (4.25), it becomes, f(r) = b- [(60b2 + 40b4 + 27b6...)r2 -(120b2 + 60b4 + 36b6)r4 5760 + 60b2r6 + 20b4r8 + O..] (4.32) Substituting Equations (4331) and (4.32) into Equation (4.4) and equating the coefficients of each power of r to zero, one obtains the following

-68 - equations: (60 - 2)b2 + 40b4 + 27b6 = 0 Co 120b2 + (60 + 4)b4 + 36b6 = 0 co 2b4b2 + 180 Cob2 + 18b6 =0 ( 33) where co = bo/5760 (4.34) Equation (4.33) can be solved for b2, b4, and b6 in terms of co. The condition defined in Equation (4.27) is then used to solve for the values of co. That is, b2, b4, and b6 in terms of co must satisfy the condition 1 + b2 + b4 + b6 = 0 (4,35) Theoretically, Equation (4.31) and Equation (4.32), after being substituted into Equation (4.4) will generate a set of n algebraic equations with n unknowns and a non-zero constant, b0 = t1. These n-equations are homogeneous but non-linear. Thus, the eigenvalues of this system cannot be determined by setting the complimentary determinante to zero. Instead, the boundary condition at r = 1 is used to determine uniquely the numerical values of all the coefficients as demonstrated previously. The exact solution could be obtained by allowing n to go to infinity. This, of course, is practically impossible. In fact, even for a relatively small value of n, say n = 4, the steps involved in solving this set of simultaneous algebraic equations are already forbiddingly complicated. The elimination of b2, b4 and b6 from Equations (4.33) and (4.35) (by Kramer's rule) results in the following function of co, called eigenvalue generating function, i.e.,

69 - 120c3(108+180co)2 + 40c3(72+108co)2 - 120(108+180c))(72+1080c0)c3 + 36co(l08+180co)(8-120c20o+120co)2 - 27c(72+108co)(8-120co+1200Co)2 - l80co(lo8+l80co)(-8+84co-3000co) - 18(8-120co+1200c) (-8+84co-300co) 0 (4.36) The roots of Equation (4.36) can most easily be obtained by plotting its value versus co.* This function is shown graphically in Figure 16, where three zeros, corresponding to co equal to 0.06702, 0.1442 and 0.2241, can be seen. Therefore, when the power series representation of g(r) is truncated up to its seventh term (r to the sixth power), these values of co, and hence bo by Equation (4.34), will give a consistent set of solutions of Equations (4.3) and (4.4). The solutions are given by Equations (4.31) and (4.32) with the coefficients determined by Equations (4.33) and (4.36). Corresponding to co = 0.06702 or t1 = 386, the solutions for t and u are t = 387x(.O0-2J6r2+2.60r4-1.44r6) (4 37) u = -x(8.6-41.7r2+52.3r4-19o2r6) (4.38) [the solution for u is obtained by differentiating f(r) and dividing f'(r) by r as shown in Equation (4.3), u = x f'(r)] r In a similar fashion, corresponding to the other two values of co, there are two sets of solutions for t and u. The possible solutions of temperature are plotted in Figure 17. As mentioned previously, since the energy equation is nonlinear, the sum of the solutions is not a solution. Therefore, it remains to determine A short IBM-704 program has been written for this purpose.

-70 - 1000 t 5 500 %... ui 0.10 0.20 Co Figure 16. The Eigenvalue Generating Function*with Three Zeros. * The function E(Co) is defined in Equation (A.36). 0.30

-71 - 1300 1200 1100 1000 900 800 700 600 500 400 z 300 200 100 0 -100 -200 -300 -400 -500 -600 t, - 9t,, C 2241 ^^*~|,< Jt,:386, Co 0.06702 tl=830, CO =1442 ~C '../~. I I.. 0.1.2.3.4.5.6.7.8.9 1.0 RADIAL DISTANCE, r Figure 17. The Solution for the Temperature with Different Eigenvalues.

-72 - which one of these solutions is the desired one. In this connection, it can be seen from Figure 17 that the only possible solution is the one corresponding to co = 0.06702. The other two solutions have a negative temperature difference within the interval of r, and this violates the second law of thermodynamics. Since it is impossible to heat the fluid in the cavity by the heat input through the wall to a temperature higher than the wall temperature. Knowing the temperature function, the Nusselt number is readily obtainable by Equation (2.16), i.e., 1 o ()r=l dx (4.39) tl O or r=1 and by Equation (4.37), Nu = 0.771. It will be of interest at this point to discuss the validity of the solutions Equations (4.37) and (4.38). The previous development has shown a method from which the solutions of Equations (4.3) and (4.4) can be obtained. By taking a sufficient number of terms in Equation (4.30) the solution can be as accurate as desired, provided that the power series is uniformly convergent throughout the entire interval 0 < r ~ 1. Indeed, the complexity of the computation would simultaneously increase with the increased accuracy. The most serious difficulty encountered is, however, to establish the convergence of the power series representation of the functions f(r) and g(r). The present problem is non-linear and of high order. Therefore, the general theory of the power series method of solving differential equations is not applicable, since it is restricted only to linear equations. Furthermore, the recurrent formula of the coefficients is not in explicit

-73 - form, thus, even the radius of the circle of convergence cannot be estimated. Nevertheless, it is interesting to note that the signs in both Equations (4,37) and (4.38) are alternating and the coefficients for r6 are smaller than for r4. The independent variable r is generally smaller than unity, and when r = 1, the condition in Equation (4.26) requires the power series to be zero. With all these favorable indications for convergence, it could be expected that the solutions thus obtained would converge. However, the rigorous mathematical proof is not readily available. The second difficulty is to justify the truncation of the power series, Should the convergence be rigorously established, the error of the truncation could be estimated. But this cannot be done here. The main concern is thus to estimate the error introduced in finding the eigenvalues due to the truncation of the power serieso It has been verified in the theory of Calculus of Variations* that the truncation has little effect on the values of the first few roots of the eigenvalues generating function. For example, the number of the eigenvalues for the present problem is reduced from infinitely many to three due to the truncation of the power series. The value of the first eigenvalue (ioe,, co = 0.06702) is only a little affected while the third value (i.e., co = 0.2214) may be hopelessly inaccurate as the result of this truncationo It is fortunate, that the present solution requires only the first eigenvalue. Therefore, it is believed that the accuracy of the solution must be very high, Again, there is no convenient method to estimate the possible error because of the complexity of the given differential system. See R. Weinstock, Calculus of Variations, McGraw-Hill Co 1952 See R. Weinstock, Calculus of Variations, McGraw-Hill Co., 1932.

-74 - With the same restriction of fluids of large Prandtl number, Lighthill obtained an approximate solution by Squire's. method.(25) His solution predicted a critical t1 of 311 which is about 20% lower than the value predicted by. the present investigation. The difference in Nusselt number is, however, bigger. Lighthill's solution predicted an overall Nusselt number of 0.364 (at tj = 311), while the present solution shows a value of 0,771. Figure 18 and Figure 19 show the comparison between his solution of temperature and axial velocity and the solution by the present method. His temperature and velocity functions are t = 311x(1-2.091r2+1-54r4-0,454r6) (4.40) u = -x(8.36-37.3r2+36.6r4-7o7r6) (4.41) As can be seen both from the graphs and the equations, the solutions from these two entirely different methods do resemble to each other. For smaller values of t1 than 387 (say ts) which is characteristic of the similarity solution when it fills the whole tube, the temperature function is given by Equation (4.37) with t1 = 387 being replaced by the value ts, and x by x - (l-tl/ts). For x <ts/tl, there is no motion and no variation of temperature (see Figure 2.) and the top of this region plays exactly the same role to the similarity flow above it as does the closed end of the tube in the case discussed above. The Nusselt number is therefore equal to Nu = 0.771 s (4 42) which is a straight line on the Nu vs. t1 diagram shown in Figure 20.

-75 - 400 350 300 250 x *' 200 150 100 50 ~o \___ _ - - ----- LIGHTHILL'S SOLUTION PRESENT SOLUTION I I \a\ n:t,=387 N \ \N \ tv = 311 __ —___^^~~ ) \.2.3.4.5.6.7.8.9 1.0 RADIAL DISTANCE,r Figure 18. Comparison of Lighthill's Temperature Solution(25) with the Present Solution.

-76 - x 0 -J -j 4 4 -II. _ -12 - -.4 3 -.2 -.1 0.1.2.3 A.5.6.7 RADIAL DISTANCE. r Figure 19. Comparison of Lighthill's Velcotiy Solution(25) with the Present Solution..8.9 1.0

cr w z Iw cn C') z -J -J w 0 1.0 0.I EXPERIMENTAL RESULTS BY MARTIN ON AIR (30) (0.6 < c < 0.7)( ----) O EXPERIMENTAL RESULTS BY MARTIN ON GLYCERINE (30)(60 < 0 <7600) ( --- —I) o THEORETICAL RESULTS BY LIGHTHILL 3 (25) ocr --- ( ---) / o THEORETICAL RESULTS BY LIGHTHILL / (25) O -— 1.0( ----) // V THEORETICAL RESULTS BY THE PRESENT/ INVESTIGATOR or ---00 ( ) - / /' I / /n I s 2 4 6 8 4 2 4 6 8 5 2 4 6 8 _____ --- - ~ _n In / / __ _ ____________ _ )6 Iv I \ RAYLEIGH NUMBER Figure 20. Nusselt Number as a Function of Rayleigh Number. Simiarly Flow Regime (e/a = 47.5).

-78 - (c) The Similarity Solutions Obtained From an Analog Computer In this section, Equation (403) and (4.4) are programed for an analog computer to obtain the solutionso As mentioned before, this set of differential equations with Uhe given boundary conditions is not very suitable for any type of numerical computationo However, with the known solutions Equations (4o37) and (4o38) the approximate values for the initial points can be estimated. Thus, the difficulty arising from the "two-point problem" is partially removed since the values of C, g' and f" from Equation (4.37) and Equation (4-38) respectively should be very close to the exact values if the analytical solutions obtained previously would be correct. In the meantime, the solutions from the analog computer can be compared with Equations (4,37) and (4,38) to estimate their accuracy since there are some uncertainties in the power series solution due to the convergence, truncation, etco To employ an analog computer to integrate the differential equations, it is most convenient to transform the equations such that their independent variable will range from zero to infinity. In many cases, the independent variable of the differential equation is the time variable and the dependent variables converge to certain constants as time increaseso The independent variable, r, in the present case, however, varies between O to 1. Therefore, in order to facilitate the programming and to avoid the possible complications in simulation, a transformation of Equations (4,3) and (4o4) deems necessary. With the exception of the nonhomogeneous term in Equation (4.3), the given system of equations is of Euler's type. Thus, the transformation, which will satisfy the requirement of the analog computer

-79 - programming, is r = eT Note when r = 1, T = 0 and when r = 0, T = o The transformation in above equation yields the following differential operators; namely, d 1 d d2 1 ( d2 dr r dT dr2 r2 dT dr2 and (4.43) d3 1 ( d +3 d2 d3 -r, (2 d +3 d - d dr3 rS dTT d- T Thus, calling F(T) and G(T) the transformed form of f(r) and g(r) in the previous section, one obtains the momentum and energy Equations (4.3), (4.4) respectively, in new variable T, i.e., F"' + 4F" + 4F' = (C - G)e-4T (444) G" - FG' + F'G = 0 (4o45 while the original boundary conditions, after this transformation, become F(O) = F(O) = G() = 0 (4 46) F(c) = 0 G(o) = (4~47) and C = F"'(O) + F"(O) (4.48) and the corresponding condition of convergence at r = 0 [Equation (4.9)] is now lim e2T F(r) = 0 (4.49) T - oo

-8o Equations (4.44) and (4.45) are programed for a Pace electronic analog computer. The electric circuit is shown in the Appendix.* The input voltages to the computer are the properly scaled values of F(0), F'(0), F"(0), G(0), G'(0) and the value of C, among which the values for F"(O), G'(0) and C are estimated from the results of the previous section. The criterion of the correctness of these values is to require;the value of F(T) approaches to zero, the values of G(T) and e2F'(T) approach to a steady value, at a fairly large T. In fact, this steady value of G (T) should be the value of t1 in Equation (4.47). The best results from the analog computer, as selected from approximately 100 trials, are summarized in Figures 21, 22, 23 and 24. Those curves are obtained according to the initial values of G'(0) = 470 and C = 170 as compared to the corresponding values of g'(l) = 307 and C = 210 from power series solution. The dotted lines on the graphs show the results from previous section. It can be seen from these graphs the general agreement between the solutions from two entirely different approach. The solutions for functions f(r) and F(T) match to each other remarkably. However, a discrepancy between g(r) and G(T) is shown on Figure 21. Nevertheless, not only the general trend between these two curves are the same, the Nusselt number of 0.706 as computed from these results reasonably checks the result from power series method of 0.771, i.e., Nu = 1 G(o) = G-( = 0.706 2 t1 2 G200T the error of which is less than 10%. * This part of the present research was done in the National Argonne Laboratory during the summer of 1961. The author was then a resident research engineer in their Reactor Engineering Division. The analog computer program is devised by Mr. R. Bare of Applied Mathematics Division of the same laboratory. Mr. Bare also kindly assisted the author to obtain the results shown in the following pages.

450 400 (9 I H? TIME(SEC), T Figure 21. Comparison of Temperature Function from Power Series Method and from Analog Computer.

.8.6.4.2 0 ii: -.2 -.4 -.6 -.8 -1.0 10 TIME (SEC), T Figure 22. Comparison of Stream Function from Power Series Method and from Analog Computer.

30 20 10 0 U. CM 0) -10 r l.O- r \_ — 0 — POWER SERIES SOLUTION. 0.A4. 7 IQ.9\p~~-0.8~~~~~~~~~~~ ~ _ _ _ __ __THE TRANSFORMATION IS 'T a Lnr -: 2F'(Tr) — 13 — POWER SERIES SOLUTION \___ __ __ __ __ ___ _ ^ ___ -- -ANALOG COMPUTER SOLUTION I I1 I -20 -30 -40 -50 C ) 2 3 4 5 6 TIME (SEC), T 7 8 9 10 Figure 23. Comparison of Axial Velocity from Power Series Method and from Analog Computer.

600..-ll. 500.. 200 __ 200 -100 -- 0 1 2 3 4 5 6 7 8 9 10 TIME (SEC), T Figure 24. Solution of G'(t) from Analog Computer.

-85 - As discussed previously, no numerical computation is suitable to evaluate a limiting process, this is clearly demonstrated by the irregularity shown in the right part of Figure 23, When T is sufficiently large, both F'(T) and e-2T are vanishing, this portion of the curve merely shows the residue statics in the computer circuit rather than the limiting value of Equation (4.49). The solution from the computer is very sensitive to the input value of C, a +5% change of its value could drastically change the form of the solutions. When this occurs, the solution of G(7) usually does not converge to a steady value (always decreasing). Also, the values from the divisor for e2T F'(T) rapidly increase and the signal of overflow is soon detected. From the theoretical solution of the previous section, this behavior of the numerical solution can be expected since the value of C actually determines the convergence of the solution at r = 0 or, equivalent to T = o. 5. General Discussion of the Theoretical Results The previous sections give the theoretical results of the laminar natural convection flow inside a circular cavity. The analysis includes the complete solutions for flow of the boundary layer type which is characterized by high values of the parameter Ra x a/2 (called tl) for any Prandtl number and two solutions from different methods for the flow in the similarity regime with fluids of large Prandtl number. All these solutions are summarized in Figure 25, A factor of about ten in the parameter, tl, separates the solutions for boundary layer type flow and the similarity flow,comparing

in2 - r 10 z -SIMILARITY -- BOUNDARY -LAYER - - BOUNDARY-LAYER TYPE FLOW REGIME - I FLOW REGIME TYPE FLOW NOT FILLING THE ENTIRE FOR FLUIDS OF VARIOUS PRANDTL NUMBER. " WITH STAGNANT LENGTH OF THE CAVITY. (SIMILARITY - - -- PORTION. FLOW TOWARD THE CLOSED END.) - SOLUTION SOLUTION OBTAINED BY LIGHTHILL (25) I I - " APPUCABLE ONLY AND APPLICABLE ONLY TO FLUIDS OF I I I TO FLUIDS OF LARGE PRANDTL NUMBER. LIGHTHILL SOLUTION (25)-7 s - LARGE PRANDTL | / NUMBER. — -- - _ _ _ _ -O-.O "0.02_ Z 0.02 { — /- -- -. ----.. --- — -- - -- __ --- ASYMPTOTIC SOLUTION SOt --- -- NUMERICAL SOLUTION POINTS TO THE LEFT OF 0 MAY INVOLVE A PHYSICAL IMPOSSIBILITY PRE / SE N ITX- I IN E S TIG A - T IO NI'II-I-I-I I I —I ---- — II ----- ---,POINTS TO THE RIGHT OF A INDICATE /SOLUTKN BY UGHTHILL (25) THE ONSET OF SIMILARITY FLOW REGIME __/.'"Iii cO I 1.0 Qi I LI 130 o10 IO 10' io' 10" t,= Ra o I t Figure 25. Summary of Theoretical Results.

-87 -fluids of large Prandtl number (Figure 25). The solution for this intermediate range of t1 has not been obtained by the author. However, Lighthill(25) presents an approximate solution for this region, which is also included in Figure 25. Due to the various approximations made in Lighthill's work, his solutions for boundary layer type flow and his "intermediate solution" do not merge at the point of transition. As previously mentioned, more polished solutions are required. Indeed, it would be of very desirable to obtain the solutions both in the similarity regime as well as in this intermediate regime for fluids of arbitrary Prandtl number. From conventional boundary layer theory, it can be shown that Prandtl number measures, as orders of magnitude, the ratio of the thickness of the velocity boundary layer and the thermal boundary layer. However, in obtaining the results for the boundary layer type flow regime for the present analysis, it has been assumed in Equations (3.8) and (3.9) that the thickness of these boundary layers is of the same order of magnitude regardless of the value of Prandtl number. The validity of this assumption depends upon the particular mechansim of natural convection flow The fluid motion in this case is entirely generated by the thermal (differential buoyancy) effect. Thus, it is impossible that the thickness of these two boundary layers should be of significant difference, since in the region where the thermal effect vanishes or changes sign, the driving fQrce for fluid motion and hence the velocity will tend to do likewise. The exact solutions for the flat plate case with fluids of different Prandtl numbers have been worked out by Ostrach,(ll) and Sparrow,(12) etc. The differential equations describing this natural

-88 - convection flow are pointwise solved numerically. From these solutions, no significant difference of the thickness of the thermal and velocity boundary layer has been observed except perhaps in the region of very high Prandtl number where a substantial temperature effect does not penetrate as far into the fluid as a substantial velocity. However, this is for the unbounded case of a semi-infinite flat plate in an infinite fluid, where as in the present case within a cylinder the wall-generated velocity is terminated by an opposite buoyancy effect in the central portion of the vessel. Figure 26 and Figure 27 from Ostrach(ll) show the close relation between these two boundary layer thickness for fluids of various Prandtl number. In spite of the various assumptions made it is believed that the solutions obtained for the similarity flow regime are reasonably accurate, since the power series solution and the solution from the analog computer closely correspond. It is realized that analog computers may not give the solutions of differential equations of this complicated nature with great accuracy. Nevertheless, the process of obtaining the series solutions furnishes the method in case solutions of improved accuracy are needed. It should be pointed out here that the integral method ( or modified Squire's method) of obtaining the approximate solutions in the similarity flow regime is not as desirable as it is when applied to the regime of boundary layer flow.* It has been mentioned earlier in this chapter that the overall heat transfer results obtained from the integral Lighthill(5) used this method to obtain the solution which is quoted in Figure 25.

'*re Id ST1J Paq-3H TlOTxall A V lmoIJ sGxaqmU TIPUeXaI snoT;J'A loj uoTqinqTax:sT a9sq.'exIaduia, ssa1uoTsuamTa -9g anS;T X t p ) =I Al Oil Xi E 9 t7 ~ ___A______ X 7EEEO 1|tEl01 9 tZ ZZ el 81 91 I ZEl 01 I 4F4iF7F11Ti I 9 -r -S oI l 80 8 0'I t 10d0 Jd I. E.,i i iNi _ I I I I I ---— \ I 1 V I -69;

-90 - LI.,U..3 ----- ----.2 I n. \. 8 0 14 1 I 20 22 24 7 Grx 14 ~S=(4 ) x.7 I x D |1L.6 cm If.-.5.4.3.2.1 0O Pr /. ~..72 L/ __1 _. _.... ~ ~'~~ —~~~.^0^^ ------ 10 0 I 2 3 4 5 6 7 GrX 114 r (- 4 ) X Figure 27. Dimensionless Velocity Distributions for Various Prandtl Numbers from a Vertical Heated Flat Plate.

-91 - method should be satisfactory since the Nusselt number is an integrated quantityo However, in the case of the solutions in the similarity flow regime, the independent variables x and r in the temperature and velocity functions are separable, and the Nusselt number as by Equation (2.16) is evaluated by integrating the value of the radial gradient of the temperature function at r = 1 with respect to x. Since the variables are separable, and furthermore, the temperature dependence on x is linear, the effect of the integration is merely to multiply the value of g'(l) (which is of course free from x since the variables are separable) by a numerical constant equal to 1/2. Therefore, a slight change in the assumed form of the temperature function could introduce a large error in Nusselt number because the value of the differential of this assumed function is needed, and differentiation quite often introduces more error than integration, The entire analysis is based on the assumption of laminar flow. It is, therefore, of interest to estimate the magnitude of the parameters for which the basic laminar flow model should break down. There are two factors in this flow regime which tend to create turbulence. First, the radial distribution of shear has a maximum in the fluid. Secondly, the temperature gradient upward is negative. Since the similarities exist between the natural convection from a vertical flat plate and the boundary layer type flow regime of the present problem, the empirical criterion of transition from laminar to turbulent flow is referenced here. Hermann(l0) found that for natural convection both from a vertical flat plate and a horizontal cylinder, transition occurred when a Reynolds number based on the maximum velocity and boundary layer thickness

-92 - was about 300~ Based on this criterion and the asymptotic solutions in the previous section (the exact numerical solutions can indeed be used for this purpose, but this would not give an explicit analytical expression for discussion), the dependence of this transition point on the parameters Ra, Pr and ~/a can be estimated. The Reynolds number, Re, according to the above definition is Re = aU (5.1) Applying the transformation in Equation (2.10) for U together with velocity function of Equation (3.9) to Equation (5.1), one obtains Re = y0(p,r) a * = 7Y0(pr) 1 & (5.2) p.- = ( ~ a where 0(P,r) is the quantity in the square bracket in Equation (3o9). Now from Hermann's criterion of Re = 300, Equation (5.2) becomes 300 = S70(B,) 1 a a and from the asymptotic solutions Equations (3.30) and (3.31), the final expression of Equation (5.2) is 300 = 0(,r) 1 x 3 ca a or x = 22.5 o0(p,r) (5.3) Several interesting results are suggested by Equation (5o3). First, it can be seen that for a given tj and aspect ratio, flow should be laminar for the entire length of the cavity with fluids of high Prandtl number, while under the same condition, the flow of low Prandtl number

-93 -could be turbulent at least for part of the cavity. Generally speaking, the quantity 0(3,r) is decreasing with increasing tl, (see Figures 6, 7, 8, 9) this suggests the turbulent motion would occupy a larger portion of the cavity as the parameter of Ra x a/2 increases.

CHAPTER IV EXPERIMENTAL PROGRAM 1. General Considerations The previous theory has predicted that for the present problem, the overall heat transfer results, expressed in terms of Nusselt number, can be correlated as a function of three dimensionless variables, namely, the Prandtl number, the Rayleigh number and the aspect ratio of the cavity. In arriving at this conclusion, many simplifications and assumptions have been involved. Therefore, it is desirable to conduct experiments to verify these assumptions so that the validity of the theoretical results can be ascertained. The experimental facility for this purpose should be so designed that these three independent variables can be varied individually. in addition to this, the facility should allow reliable measurement of the quantity of total heat input for computing the dependent variable, i..e,, the Nusselt number. Naturally, the geometry of the test section and the physical conditions of these experiments should be compatible to those used in the analysis. Based on these prerequisites, the present experimental facility has been designed. The details are discussed in the following section. The solutions given in the previous sections have covered a wide range of Prandtl number. (The numerical solutions cover a range of Prandtl number from 10.0 to 0.02, and the asymptotic solutions are good for all non-zero Prandtl numbers at large Rayleigh number.) Therefore, -94 -

-95 - the fluids used as test media should also cover a large range of this parameter. The work of Martin, Cohen, et al.(29'30'31) has used glycerine, water and air under the same conditions and in the same geometry (the Prandtl numbers of which are 66.0 to 7,000, 1.8-9o0, and 0.7-6.0 respectively). These tests covered quite completely the range of large Prandtl number fluids. In fact, they were intended to verify Lighthill's solution(25) which is applicable only to fluids of large Prandtl number. In the present work, water and mercury are chosen as the test media. The purpose of the water tests is to develop the experimental techniques with the present facility in a case where comparison with similar data is possible. Mercury has been selected because of its low Prandtl number, which is approximately 0.03 at room temperature. No other fluid available at room temperature has Pr appreciably below 0.73 (Approximately the value of air). However, the low Prandtl number is typical of all the liquid metals, which are of present technological interest, but, with the exception of mercury, only are liquid at high temperature. NaK or Na would considerably extend the range of available data (Pr for NaK is" 0.003), but the complication of the required high temperature is involved. For basic research, mercury is considered sufficient. Two test sections have been used: one with an aspect ratio of 133 and the other with 22.3Due to the equipment limitations both are of the same diameter. * Defined as the ratio of the length to the diameter of the test section.

-96 - Rayleigh number can be varied in two ways after the diameter of the test section is fixed; i.e., variation of the temperature difference between the wall of the cavity and the reservoir (which is the quantity used to define the AT in Ra), or of the average fluid temperature, since the fluid physical properties are functions of temperature. The test range for Ra is 105 to 106 for mercury and 105 to 107 for water. The range is limited by the capacity of the test facility and the desire to obtain laminar flow to match the assumptions used in the analysis. Since mercury vapor is toxic in extremely small concentrations, safety precautions must be taken into consideration. A mercury vapor concentration meter was placed next to the test section to detect any possible mercury vapor leaking out of the system. In addition, the facility is located in a test cell equipped with an automatic ventilation system so that if any mercury vapor should leak, it would be carried out before reaching a toxic concentration. Due to the specified closed geometry, the low magnitudes of velocity, and the opaque nature of one of the fluids, it is almost impossible to measure the fluid velocity, though this would be very desirable. Ingenious devices are available in the literature, such as those by Saunders, Schmidt, etc.,(8,9) but their methods require the use of smoke traces in an air stream or an open geometry. Hammitt(2) measured the velocity of natural convection flow in a closed cylinder by timing the movement of dye particles, but again his method is applicable only to a transparent test section with transparent test fluid. Several devices were considered

-97 - in a preliminary fashion, but none of them has proven practical. Fortunately, however, the local fluid temperatures can be measured to a reasonable accuracy. Therefore, not only a fairly accurate overall heat transfer result can be obtained, but also the radial profiles of the temperature functions are obtainable. 2. The Experimental Facility Figure 28 is the schematic of the experimental facility. The test section (A) is welded to the center of the bottom of the reservoir (B). The entrance to the test section has been carefully machined after welding, so that the disturbance to the fluids due to irregularities of the edge is reduced to a minimum. A cooling coil is mounted inside the reservoir to maintain the reservoir temperature constant. The diameter of the coil is sufficiently large so that the cooling surface is reasonably removed from the test section entrance to assure an approximately uniform temperature of the fluid entering the test section. The spiral coil is fabricated from 1/4" O.D. stainless steel tubing, approximately 6' long. This provides a total heat transfer area of about 57 sq. in. A rough calculation has shown that with this amount of heat transfer surface, the maximum temperature difference between the inlet and the outlet of the coolant would not exceed 5~F even when the facility is operated at maximum heat input. This small temperature difference proves to be very important in stabilizing the reservoir temperature. Both the inlet and outlet of this coil are welded to the wall of the reservoir. Thermocouples are inserted so that the coolant inlet and outlet temperatures can be measured during the experiment.

A. Test Section B. Heat Reservoir C. Movable Thermocouple Probe B D. Cooling Coil E. Coolant Recirculating Pump F. Coolant Heat Exchanger _ G. Insulation M. High Precision Voltmeter N. High Precision Ammeter TO Q THERMOCOUPLE P. Potentiometer for Thermocouples SWITCH FO O w co COOLANT TO FLOW G METER POWER SUPPLY I= rTOQ SYSTEM COOLANT (The regular triangles indicate the position of various thermo- A couple beads). F Figure 28. Schematic of Experimental Set-up.

-99 - The reservoir is constructed from a piece of stainless steel pipe of 6" nominal diameter and 8" length, A matching blind stainless steel flange is welded to one end of the pipe and a matching slip-in flange to the other. A hole is drilled at the center of the blind flange for the insertion of the test section. The thickness of the flange in the neighborhood of this hole is reduce to 1/4" (from the original thickness of 1"). In this way, the test section is secured firmly to the reservoir, and the heated length of the test section is maintained at a maximum. Another blind flange is used as the cover for the reservoir. Two small holes are drilled in this flange: One for filling (plugged during operation), the other is used to hold the movable thermocouple. An "O"-ring seal is used to prevent leakage at this point since the thermocouple rod is a moving pieceo Two test sections were used. Both are of 1-1/2"0 oDo stainless steel seamless tubing of 0.075" wall thickness. The shorter section is 18.0 inches long and the longer one is 30.0 inches. The bottom end of each test section is sealed by a thin stainless steel plate and the open end is welded to the reservoir as described previously, A small hole of "1/16" diameter is drilled in the sealing plate which is welded to a needle valve to provide a drain. Although both the test sections have been carefully machined before being used, the presence of small irregularities on the walls is almost inevitable. However, it is believed that these are so small as to have possibly only negligible effect on the motion of the fluid, since the scale of the velocities inside the cavity is very small. Figure 29 is the assembly of the test section.

-100 - PLUG MOVABLE THERMOCOUPLE HOLDER 11I I ~~I RESERVOIR I ~~$I 6"O.D. S.S.PIPE I COOLANT 8,LONG COOLANT OUTLET INLET I CLNI I I l l I I I 0 ATCOOLIN G 0 I INSULATION 0 STATIONAR COI L ATION VBOX THERMOCOUTP- I S S LE I3/ I\~~~~~~ I"v~~~~~~~ I I 2 D SAIHEATING ITHICKVV T NTAPE I INSULATON - o.075 " DRAIN " SEALING I12"NEEDLE PLATE 0.125" VALVE THICKNESS SCALE iT3/8P Figure 29. Assembly of the Test Section.

-101 - The mathematical model described in Chapter III requires a steady and uniform temperature at the walls of the cavityo Physically, such a condition can only be obtained by applying a wall heat flux which is nonuniform in the axial direction. The required heat flux distribution may be approximated by using separate heaters and suitably adjusting the power input to each heated section. Theoretically, an infinite number of heaters are required but this is indeed physically impossibleo Electrical resistence heaters are used. This proves effective both for construction and in controlo A possible alternative is the use of steam heaterso After a brief consideration, this idea did not appear feasible in the present caseo The rate of heat input to the system would be hard to control with steam. Even though the outside wall temperature were maintained uniform, the rate of condensation of steam (ioeo, heat flux) would vary along the walls so that the inside wall temperature would not be uniform. Also, a suitable steam supply was not readily available. There are five electrical heating coils for the long test section and three for the short oneo Each of these heaters is formed by winding resistance tapes around the insulation which is wrapped on the tube wall. The heating tape is 1/161 wide and its thickness varies from 0.003" to 0.005I? depending on the designed capacity of the particular heater. The glass cloth and mica electrical insulation is used around the tube wall. The mica is sandwiched between two layers of glass cloth. The total thickness of this insulation is less than 1/32", and throughout the

-102 - whole experiment, no burnout has been detected. The spacing between two windings of the coil is kept less than 1/32". Thus, in view of the relatively high thermal conductivity of stainless steel (i.e., compared to the electrical insulation) and the thinness of the electrical insulation, a fairly uniform temperature of the tube wall could be attained. The heater lengths are 4", 6", and 8' for the short section and 3 ", 3 ", 6 ", 8 ", and lO "respectively for the long one(listing from the top of the tube). The sharp increase in the length of the heaters is due to the distribution of the temperature and heat transfer along the wall as shown by previous theory, i.e., both the change of temperature and the heat transfer is smaller toward the bottom of the tube. Therefore, even with a comparatively longer heater, the lower end of the tube can still be kept at uniform temperature. There is no heater across the bottom of the test section for either tube, although the boundary conditions of the basic problem do require a uniform temperature for all the walls and some heat transfer through this portion is inevitable. There is no convenient means to attach a heater to this wall and also the amount of heat transferred from this portion of the tube is small as may be seen both from the theoretical consideration and the physically small heat transfer area.* A portion of the heater construction is shown in Plate 1. Six thermocouples were welded along the outside walls of the short test section (which was constructed first). It was found, however, their readings were * Note the definition of Nusselt Number in Equation (2.16) does not include the amount of heat transferred through this wall.

-103 - *~:"'~~~~:~!ilii'~ ~~~ A-.:~:i:: ~::;i_~::. A-:D -~~~~~~~~~~~~~..^ 4 CB^1-"^"- i ~ *^^ ^sdif -.,. SU s _+//W ~-~c Plate 1. Photographs of the Heaters.

-104 - affected by the local high temperature of the heating tape to be of significant assistance. Therefore, for the longer test section, the thermocouples along the outside vertical wall are eliminated with the exception of locations at the test section entrance and the bottom of the sealing plate. This gives a rough check on the movable thermocouple at the top and the bottom of the test section and an indication of temperature at a point which it cannot reach. However, major reliance is placed for the entire investigation on local temperature measurements from the movable thermocouple. Figure 30 is the assembly of the movable thermocouple. It is composed of a 1/8" O.D. stainless steel tubing which shields iron-constantan wires of approximately 0.008 diameter, insulated by magnesium oxide powder. The assembly is about 5 ft. long. A 3/16" O.D. by 4'-6" stainless steel tube houses the smaller tube to give the required rigidity and guidance to the assembly. An eccentric three-pronged circular spider with a diameter of 1.35" (which is the I.D. of the test section) is located at the bottom to maintain the required spacing between wall and thermocouple rod. The lower end of the thermocouple rod is bent to a right angle. A bead is formed at the tip which is 0.3375" from the centerline of the rod (this distance is equal to one quarter of the inside diameter of the test section). At about 1/2" above the bend it is located by the spider. The eccentric center of the spider, and hence the center of the thermocouple rod, is located half way between test section centerline and the wall. The remaining portion of the thermocouple rod is housed in the 5/16" tube and is held at the both ends of the tube in two bearing sleeves so that it

-105 - TOP OF RESERVOIR HOLDING DEVI TO POTENTIOMETER / - a$tS. THERMOCOUPLE I CONDUCTOR C HEX BOLT T BORED TOA - P.G. DISKS(2B A 8RPAPER DISK BETWEEN - DRILL ROD I OD SS TUBING I+ I HEX BOLT BOREDTO — " " " P tLE XX J% L | TOP ASSEMBLE BOTTOM ASSM ICE Figure 30. Assembly of the Thermocouple.

can rotate freely with respect to the housing tube, but can only be moved vertically together with the housing tube. On the top of the upper end of the thermocouple is a dial which indicates the angle of thermocouple turning. Along the wall of the housing tube are marks of half inch interval. Therefore, the exact location of the thermocouple bead can always be known. The housing tube is then secured in a specially made coupling screwed on the cover of the reservoir. Three "O"-rings in series are used to prevent vapor leaking. These allow rotation and vertical movement of thermocouple. A three-legged supporting device is erected on the top flange to provide a rigid support for the required vertical gliding motion of the assembly. The centers of the spider and of the coupling are carefully aligned so that the centerlines of the thermocouple and of the test section are parallel. Since this assembly, in spite of its length, is supported at both ends, its vertical motion is positive and only vertical motion is possible. By this arrangement, the thermocouple bead can move along a circular path of a diameter equal to the radius of the test section as well as along the axis. Since axial symmetry in the flow is assumed, the temperature measurements on the half circle of the locus of the bead, defined between the wall to the center of the test section, should be the radial temperature profile of the fluid at a known distance from the bottom of the cavity. Furthermore, since two sets of mutually symmetric readings are obtainable by rotating the bead to a full circle, the assumption of axisymmetry can be checked.

-107 - The presence of this movable thermocouple probe in the fluid would certainly affect the fluid motion. However, it is believed that this effect is small and can be neglected. Since the cross-sectional area of the probe occupies less than 3% of the total area of the cavity (the legs of the spider are made of 1/16" stainless steel rod), it seems reasonable that it should not cause a great change in the fluid motion. Indeed, since the probe is off center, it would change the symmetry of the motion. In this connection, the experimental work of Hammitt is cited. He used an off centered thermocouple probe, much like that here described, to measure the temperature field of natural convection in an entirely closed circular test section with internal heat generation. As a check on the possible effect of the probe, a second "dummy" probe was inserted in a series of tests, with all the physical conditions kept the same. The temperature profiles were then measured in the same facility with two probes in place. The additional probe served only to preserve the geometrical symmetry in these tests. Several sets of comparable data thus obtained have shown no noticeable discrepancy. Therefore, it can be concluded that this temperature measuring device does not distort the physical phenomenon to a detectable degree. Since the spider is about 1/2" above the thermocouple bead and cannot come out of the test section because of its diameter. The temperature reading at the inlet of the test section cannot be measured by the same device. A fixed thermocouple, similar in design to the movable thermocouple, is used at this position, one end being attached to a

-108 - coupling welded to the wall of the reservoir. In order to avoid possible confusion, these two thermocouples are made of the same material. The whole unit is mounted vertically. The top section of this structure (including the reservoir) is enclosed in an insulating box, about two inches of glasswool insulation being used on each side. Around the test section itself and the drain valve, are wrapped with a double 4 inch layer of commercial thermal insulation. Three layers of aluminum foil are sandwiched between the insulation layers to reduce radiation heat loss to a minimum. Three thermocouples are imbedded in the insulations to measure its temperatures. It has been found that heat loss is tolerable in most of the tests, since the temperature difference between the ambient and the surface of the insulation has never been more than 1.0~F. The estimated maximum heat loss corresponds to about 5% of the heat input. All the thermocouples are led through a manual multi-point switch to a Leeds and Northrup potentiometer so that temperature readings are taken one at a time. However, the delay between readings is not sufficient to cause severe inconvenience. The cold junction is immersed in an ice bath. Individual variable transformers are used for power supply to each heater, since each heater assumes a different heat input to the test section, the design capacity of these being selected according to the approximate heater requirements; i.e., the transformers for the top heaters are much bigger than those for the lower ones. The largest transformer (for the topmost heater) has a maximum capacity of 1.12 kw (8.0 amp x 140 volts) while the smallest (for the lowest heater) only 0.1 kw. The maximum total heat input is approximately 2.8 kw, but under no circumstance could this be realized due to the distribution of the required heat input.

-109 - The output from the power transformers is led through a control panel equipped with voltmeters, outlets, switches, indicating lights, etc. A voltmeter is connected to each power transformer to indicate the approximate voltage across the corresponding heater (there are thus five such voltmeters). A high precision Weston voltmeter and ammeter are plugged into the sockets on the board so that they can be connected to any of the heater circuits. Thus, precise readings of the power input to each heater are made one at a time. The regular voltmeters merely indicate roughly the condition of the heaters during the test. Since the capacitance and inductance of the heaters is negligible, the power factor is very close to one, total power input is accurately calculated from the product of current and voltage. Conventional safety devices are provided for the protection of the precision meters and transformers, etc. The cooling system is of the recirculating type although initially a much simpler device using tap water directly was employed. However, this simple coil presented disadvantages as explained below so that it was necessary that the design be changed. The temperature of the tap water is usually well below the temperature of the reservoir. This big temperature difference between the wall of the cooling coil and the reservoir fluid appeared to cause instabilities in the fluid entering the test section. In some instance, it is almost impossible to achieve a desired reservoir temperature. Also, the flow rate of tap water in the laboratory has been found quite unsteady, and consequently, to maintain a steady temperature of the reservoir is very difficult. In some cases, this reservoir temperature is so unsteady. as

demonstrated by the rapid fluctuation, that a reliable measurement is impossible. This is more pronounced in mercury tests than it is in water tests. By adjusting the flow rate of the tap water in the auxiliary heat exchanger of the redesigned system, the temperature of the recirculated coolant can be controlled and the undesired large temperature difference eliminated. Because of the thermal inertia (or the time constant) of the heat exchanger, a large part of the tap water flow rate fluctuation is damped out and not passed to the secondary coolant which directly affects the steady state of the reservoir temperature. There is yet one more advantage for this type of cooling system, Since the properties of the test fluid are functions of temperature, the Rayleigh number, and hence the parameter tl, of the tests can be varied at a constant power input. In other words, it is possible to hold the power constant (and thus, approximately the temperature difference between the cavity wall and reservoir remains unchanged), and to change the parameter t1 by adjusting the coolant temperature via tap water flow-rate. This capability of the facility proves helpful in attaining the desired data. The cooling system is composed of a cooling coil as mentioned previously, a small centrifugal pump, a tubular single pass heat exchanger, a rotormeter, and a small glass jug used as the reservoir of the recirculated coolant. The coolant leaving the coil is passed to the heat exchanger where its excess heat is carried away. It then flows into the glass jug

-111 - the outlet of which is connected to the centrifugal pumpo It then is delivered to the coil again thus to complete the cycling process. The flow rate of the coolant (ordinary tap water) in the shell side is indicated by the rotormeter which is ranging from 0.05 to lo0 gpm, in accordance with the designed capacity of the total heat input to the test section. The power to the driving motor is adjustable so the speed of the pump can be controlled, but this device has never been utilized since the control of the tap water flow rate is easier and effective enough. Plate 2 is a photograph of the facility. The arrangement of components is the same as those shown in the schematic drawing, Figure 28, which can be used as the referenceo 3. The Experimental Procedure As mentioned before, water and mercury have been selected as test media and two test sections of different length have been constructed) Water tests were conducted only in the short test section (the 18" section) while the mercury tests were carried out in both sectionso With the exception that the recirculating cooling system was used only with the mercury tests the experimental facility and general procedure of the tests have been essentially the same throughout the entire programo After the test section and all the auxiliaries were assembled, the heat reservoir as well as the test section was thoroughly washed with acetone to remove grease and other dirt. The movable thermocouple and the top cover flange were carefully aligned with the centerline of the test

-112 - k --- — i I: i j:: i i i 1a1i i r Iii f T I i i.il-j 1 ii::::~ia x: i 1 ii i I i i 1:: 1 iI I, _:::i:::ii:I i I1 i `.; Plate 2. Photograph of the Experimental Facility.

-113 - section, and then the cover is assembled to the reservoir. The whole assembly, erected on the upright structure, was carefully leveled. After the fluid is poured into the system, the hole on the top flange is plugged. Distilled water and triply-distilled mercury were used. The instrument was carefully calibrated before used. To start a test, the lowest heater is turned on first, then the recirculating pump is started (this is, of course, eliminated for water tests), with zero coolant flow rate. "Coolant flow rate" refers to the flow rate of tap water in the heat exchanger. The recirculated coolant is at a constant flow rate for all the tests. After an hour or so, when the whole system is warmed up, the subsequent heaters are successively turned on. In the meantime, the movable thermocouple is positioned to measure the wall temperatures at frequent time intervals. In accordance with the thermocouple readings, the heat input to each section is adjusted until an approximately uniform wall temperature is reached. This adjustment takes 4 to 8 hours depending on the total heat input. Generally speaking, it takes a longer time to adjust the heaters for larger Rayleigh number runs than it is for the smaller. For the mercury tests, the mercury vapor detector is used as a safety device. Since the operating temperature is relatively low and the system is well sealed, the vapor detector has never registered any substantial reading. The maximum temperature for the mercury tests is approximately 350~F, but the reservoir temperature is only 170~F. The corresponding saturated vapor pressures are 9.65mm and 0.03 respectively. Except toward the entrance of the test section, where a mixing of cold fluid entering the orifice and the hot fluid leaving the cavity is

-114 - inevitable, the wall temperature readings are quite steady for most of the water tests. The unsteady runs are those of large Rayleigh number for which the flow is believed to be very turbulent. For the mercury runs, almost all of the wall temperatures are unsteady for all the tests. Thus, the wall temperature is considered uniform and constant if the reading of the movable thermocouple is fluctuating within the range of the scale of the potentiometer (Leeds and Northrup type-8662) about the desired value. This is approximately equivalent to a fluctuation of + 5~F. After the desired time-mean wall temperature is reached, the power input to each heater is measured using the precision voltmeter and ammeter, and the total heat input calculated. The insulation temperatures are checked in order to estimate the heat loss. Based on the readings from the thermocouples embedded in the inner and outer wall of the insulation, the heat loss is estimated to be less than 10 BTU/hour. This is small compared to the total heat input (which is of the order of 200 BTU/hour so that the relative maximum heat loss is about 5%). Although the temperatures of the coolant inlet and outlet are measured, the overall heat balance of the system has not been attempted. In the water tests, where direct cooling system was used, the measurements of the coolant flow rate were not accurate enough that a meaningful heat balance could be made. Since no flow meter is incorporated with the recirculating system, the recirculated coolant flow rate is an unknown. However, this quantity is of no real interest since only the net heat transfer to the system through the test section wall is of main concern. The tubings

-115 - of the cooling device are made of thick-walled tygon material. They are therefore reasonably remote to the small change of the room temperature which might be expected to have an effect on the steady state in the heat reservoir. When equilibrium is attained, the wall temperature and the centerline temperature are measured at several axial positions with the movable thermocouple. A radial temperature profile at a fixed axial position is also made. Due to the unsteadiness of the temperature near the throat (the entrance of the test section), this axial position has always been taken near the middle point of the test section. Altogether 47 tests have been made in this fashion: 14 with water, 14 with mercury in the shorter test section and 19 with mercury in the longer test section. Since the major interest is in the heat transfer rate, the radial temperature profile measurements have been omitted in some of the testso(Table II is a summary of the test results.) 4. Observations As shown by the theoretical treatment, the results of this experiment should be correlatable in terms of Rayleigh number and Nusselt number with Prandtl number and the aspect ratio of the cavity as parameters. The data from the present investigation can then be grouped into three categories, i.e., water in the short test section (Prandtl number form 2.0 to 6.0, I/a = 26.6), mercury in the short test section (Prandtl number from 0.015 to 0.02 with same aspect ratio), and mercury in the long test section (same Prandtl number with aspect ratio equal to 45.5). The temperature differential between reservoir and wall of the cavity is used to define the AT in the Rayleigh and Nusselt numbers.

TABLE II SUMMARY OF EXPERIMENTAL DATA Water Runs (Q/a Pr x GR = 26.6) Nu Mercury Runs (2/a = 26.6) Pr x GR x 10 5 Nu 4.01 x 5.1 x 5.95 x 6.61 x 6.92 x 7.8 x 9.11 x 9.52 x 1.032x 1.04 x 1.1 x 1.45 x 1.81 x 2.01 x 106 106 106 106 106 106 106 106 107 107 107 107 107 107 6.51 7.32 7.54 7.52 8.50 8.91 8.02 10.10 10.50 10.01 10.02 9.35 10.05 10.07 1.2 1.22 1.32 1.55 2.52 2.55 2.74 3.42 4.40 5.60 6.80 6.85 8.02 8.72 0.251 0.222 0.140 0.223 0.20 0.245 0.272 0.270 0.285 0.310 0.29 0.25 o.308 0.201 Mercury Runs (2/a = 44.5) Pr x GR x 105 Nu 1.80 1.98 2.10 2.30 2.42 2.90 3.33 3.65 4.05 4.93 5.65 5.95 6.40 7.55 7.96 8.62 9.60 9.62 10.92 0.192) 0.278 0.202 0.120 0.142 0.115 o.141 0.122 0.124 0.144 0.105 0.126 0.155 o.14o 0.124 0.161 0.142 0.155 (unsteady runs, zero coolant)

-117 - Since both these temperatures are not quite uniform, an arithmetic mean wall temperature and the time mean temperature reading from the thermocouple at the throat are used for this purposeo Since the physical properties of the fluids are functions of temperature, their values used in computing these dimensionless numbers are taken at the mean temperature between the reservoir and the cavity wallo Furthermore, since the Prandtl number for each test of the same medium also varies due to this temperature dependence, the final correlations from this experiment are actually expressed within a small spectrum of variation of this parameter. The range of this parameter for each group of tests can be seen in the prior section. Values of the various properties of water and mercury are taken from Heat Transmission and from Liquid Metal Handbook,Reference 46 and 42 respectively, Some of the applicable properties of mercury as functions of temperature are included in the Appendix. Figure 31 is a summary of the experimental results. The points on the top of the graph are the data from the water tests, This includes the data from the present investigation and also that by Martin and Cohen(29) in a facility which is generally similar except that the aspect ratio is 4705. The four points in the middle of the graph are the mercury data for a flat plate experiment taken by Saunders (9) Since when the Rayleigh number is sufficiently high, the results from the flat plate analysis should converge to the present theoretical results, these data points are also included for comparison, The points in the lower portion of the figure are the mercury data from the present experiments. Two different symbols represent the data from tests with different aspect ratios (26.6 and 45.5), The data at the left portion of Figure 31 are

102 10 I I I I I I ill I I I I I I I I I 10 THEORETICAL SOLUTION (LAMINAR ANALYSIS) -- - - X MARTIN'S DATA FOR AIR, a.=0.70, I/a =47.5 0 MARTIN'S DATA FOR WATER, 9.0> 0'>1.8, /a =47.5 _ * SAUNDER'S DATA FOR MERCURY, v. 0=.02, FLAT PLATE CASE PRESENT INVESTIGATION DATA FOR WATER 6.0> 0>2.0 V,,,, MERCURY (t/o =26.6) I I I IIIII~ I~ I. | -II...,I:.II'...I.....-. I- -a = 1.0, 1/ 1=26.6 A,,,,,, (/a =45,5).0/ 6 F I I I 111111 1|1|1 --- —-— 1-1-L_______ — _g,g 3 > __""'_ I; 'l t/o:26.6': i,? *)~~~ ~Ij ___~ r~~~ _ 0.015<0' <0. 02 __ _ z I H I 1.0 9, '. - I\ ~Oq I "A *A^ _A _ t,/a 45.5 A A I - 0.1 lo' lo' 10o o108 Ra = GR x PR Figure 31. The Summary of Experimental Data.

-119 - test results from air by Martin and Cohen(29). No test has been conducted in the similarity flow regime in the present investigation, however. The dotted lines are the results of theoretical considerations. As can be seen from this graph, the values from the theoretical predictions are higher than the experimental results in all the cases, and the deviations increase as the Prandtl number of the fluid decreases. It is believed that this discrepancy is due to the fact that the flow inside the cavity under the experimental conditions is actually turbulent while the theoretical analysis is based on the assumption of laminar flow. The water data of the present investigation agree reasonably well with those of Martin and Cohen, which indicates that possible errors introduced into the experiments are not significant. As shown by the theoretical curves, smaller aspect ratio will increase the Nusselt number. The present data, with an aspect ratio of 26.6, lie above the data by Martin and Cohen, whose facility has an aspect ratio of 47.5. The experimental results from mercury deviate much more from the theoretical consideration than those from the water tests. With the exception of three points, the aspect ratio still plays an important role. The test results from the short section are well above the results from the long test section. The three higher points which are the exceptions mentioned correspond to zero coolant rate. The reservoir temperature was thus very close to the wall temperature so that the magnitude of AT was considerably reduced, giving rise to a comparatively small Rayleigh number and a large Nusselt number. Since experimental conditions for these tests are considerably different from the rest and from those required in the analysis, it is very hard to make qualitative comparisons or to draw any quantitative conclusions.

-120 - It has been previously mentioned that both the fluid temperature and the wall temperature are time-averaged values and considerable fuctuations occur. This same phenomenon has been observed by Hartnett(4) Hammitt(29) and Martin(30). For both test fluids used in the present investigation, it has been observed that the largest fluctuations occur at the entrance of the test section and are reduced as the closed end is approached. The oscillation is much larger in the mercury than in the water runs. This agrees with the observations of the previous investigators(32). Hartnett measured the frequency and amplitude of these oscillations in his constant flux experiments(32), but no conclusion could be made from these measurements, other than that these oscillations were of non-linear type. No time-dependent measurements of these oscillations have been made in the present investigation. However, a test in mercury was made with zero coolant flow rate and zero heat input to all the heaters except the lowest one, which drew a total heat input of approximately 50 watts. Under these conditions, it was hoped that the disturbance from the mixing of the cold and warm fluid at the entrance of the test section would be eliminated. However, the oscillation of the fluid temperature was still observed in this case. Therefore, it can be concluded that natural convection flow with low Prandtl number fluids is exceedingly unstable under almost any condition. The present analysis, as well as the analysis by Lighthill, using an integral method, can not predict this observed instability. Rather, a stability analysis is required of the type used by Chandrasekhar and Elbert(37), who treated the simpler situation of a thin layer of fluid heated from below.

-121 - The calculations for this simple geometry predicted a cellular type of motion for normal Prandtl number fluids such as air, water and oil when the temperature differential between the bottom and the top surface exceeds a critical value. For a fluid of low Prandtl number, such as mercury, the analysis of Chendrasekar and Elbert(37) predicted a situation of "over-stability", in which local fluid temperatures oscillate, and the simple cellular pattern predicted for water is no longer valid. Such predictions have been verified experimentally by Fultz and Nakagawa(38) working with a thin layer of fluid. A similar situation is very likely to occur in the more complicated unstable flow geometry treated in the present experiment. For the water tests, the reservoir temperature is quite stable except in the proximity of the entrance of the tube, where a fluctuation of the order of 5~F has been observed for the entire range of tests. A test was made with the cover of the reservoir and the movable thermocouple removed, and it was observed visually that periodic discharge of hot fluid from the cavity existed. This shows that a steady boundary layer with upward velocity around the cavity and a steady downward velocity in the core which are assumed by the theoretical analysis may not exist at all. Instead, a very complicated mixing process may have taken place in this region. The highest temperature difference between the reservoir and the cavity is about 80~F, and at this temperature, the fluid temperature in the cavity is about 205~F. Since the system is not pressurized, this is the highest Rayleigh number water test possible with the present facility (Ra = 2.01 x 107, Nu = 10.7).

-122 - For the mercury tests, the reservoir temperature is very unsteady, especially in the vicinity of the test section entrance where the temperature fluctuation is usually of the order of 20~F. Although it has not been visually observed, the periodic discharge of warm fluid must have been in existence. This is evidenced by the periodic variation of the thermocouple readings. The period of this temperature shift ranges from approximately 1 minute to 20 seconds as timed by a stop-watch. Generally speaking, the period is longer and the fluctuation smaller for small Rayleigh number tests and vice versa for large Rayleigh number runs. The temperature of the fluid in the reservoir is very sensitive to the coolant flow rate, and it is very difficult to maintain a temperature differential between the fluid in the reservoir and the fluid in the test section less than about 60~F. The movable thermocouple measurements show that this temperature differential exists only within the top two or three inches of the test section. The rest of the tube is in a more or less stable condition, though a temperature fluctuation of 5 to 10~F has still been observed for all the mercury tests. The sharp response of the reservoir temperature to coolant flow rate suggests the existence of excellent heat transfer in the reservoir, Since the reservoir temperature is fluctuating, it can be assumed that the fluid motion in the reservoir is far from negligible. Therefore, it is probable that the fluid motion in the test section is also affected by the turbulence of the entering cold fluid causing some of the additional instability in the test section, In some tests, attempts were made to reduce the tap water coolant to bring the reservoir temperature up to as high as 170~F to reduce the temperature differential. It was hoped that turbulence in reservoir

-123 - might thus be decreased. This reduced the temperature fluctuation of the reservoir fluid only slightly and the oscillations at the throat of the test section still persisted. The theoretical model employed in the analysis assumed a constant core temperature extending from the top of the test section to the bottom and equal to the reservoir temperature. However, the actual observation is quite different from the assumption. In most cases (both water and mercury tests), the core temperature is very different from the temperature of the fluid in the reservoir but close to the wall temperature, although in the water tests, the radial temperature differential between the centerline and the wall of the cavity is larger than it is for mercury tests. But by no means can be said that the core temperature is even relatively close to the reservoir temperature. It appears that the cold fluid of the reservoir is rapidly mixed with the warm fluid within a very short distance after entering the test section (2 to 3 inches). This greatly reduces the effective temperature differential in the cavity and, hence, the total heat transfer through the wall. Nevertheless, the Rayleigh number and the Nusselt number are computed based on the temperature differential between the reservoir and the wall, which does not actually represent the effective temperature differential as far as heat transfer is concerned. Figure 32 shows a typical axial temperature profiles of the wall and the core. Figure 33 shows a typical radial temperature profile.

-124 - 30 28 26 24 U) w I UJ z I. 0 or 0 z I.rl U) Z 0 a\ 22 20 18 16 14 12 -— WALL TEMPERATURE - - — CENTERLINE TEMPERATURE o BEYOND THIS POINT THE VARIATION OF POTENTIOMETER READING IS TOO BIG FOR ANY RELAIBLE MEASUREMENT I I~.........: 10 8 6 4 2 n 200 220 240 TEMPERATURE,OF 260 280 Figure 32. Wall and Centerline Temperature Profile. (Mercury Run No. L-M-2 Reservoir Temp. = 162~F Power Input = 383 watts)

-125 - 1.0 -J 4 w: 0.9 D: I0.8 0 0. -- -1 0.2 0.3 0.4 0.5 RADIAL DISTANCE, INCHES Figure 33. Radial Temperature Profile. (Mercury Run No. L-M-3 Reservoir Temp. = 172~F Power Input = 353 watts Profile Taken 10" from the Closed End)

-126 - Due to this rapid mixing process at the throat of the test section and if the wall temperature is maintained constant, only the wall near the throat is effective for heat transfer. Figure 34 shows a typical heat transfer distribution for the short test section. As can be seen, about 85% of the total heat transfer is through the top four inches of the test section (from a total length of 18 inches). This situation is even more pronounced with the long test section. From the temperature measurements and the power input data, it is indicated that the heat transfer rate is very low for almost the entire length of the tube except for the top few inches. This reduction of the heat transfer rate below theoretical expectations is due to the mixing effect of the turbulence in the flow near the inlet of the cavity which reduces the effective wall to centerline temperature differential very considerably, Also, for the entire tube the cross streams in the turbulent flow field reduce the effective density gradient. Within the range of the mercury tests, it is believed that the entire fluid is in turbulent motion. Due to the inherent instability of the motion analogous to Chendrasekhar's prediction, the onset of turbulent flow could occur at a very small Rayleigh number (as evidenced by the test described previously), the fluid motion is even more perturbed by the entering cool fluid from the reservoir which may be the cause of the effective mixing process at the throat. There may be another mechanism which affects the fluid motion in the mercury tests. It is recalled that the test section is wound with resistence wires as heaters. When a current is passed through the heaters, it induces a weak magnetic field inside the test section (the current

-127 - 18 NO.I HEATER (4" LONG, 139.5 WATTS TOTAL POWER INPUT) 15 o') Q 3 0 z 10 z IF 5 5 NO.2 HEATER (6" LONG 21.7 WATTS TOTAL POWER INPUT) + I NO.3 HEATER (8" LONG 6.4 WATTS TOTAL POWER INPUT) O L 0 10 POWER INPUT, WATTS/UNIT 20 TEST 30 35 SECTION LENGTH Figure 34. Typical Heat Input Distribution. (Mercury Run No. S-M-1 Reservoir Temp. = 74~F Power Input = 167.6 watts Wall Temp. = 110~F

-128 - through the top heater is about 8 amperes and there are 50 windings). Thus, a "pinch-effect" on the mercury is exerted giving an axially non-uniform radial pressure gradient. Although this gradient is small, so are the driving forces for the flow. Thus, there may be some non-negligible perturbation. Although similar heat input distribution has been shown in the water tests, the temperature oscillation is much less in the lower portion of the cavity. This could suggest a locally laminar flow. In the present facility, visual determination of the fluid motion is impossible, especially in the lower portion of the test section. However, Hammitt(2) observed a partially turbulent and partially laminar flow pattern in his test facility. In general, the Nusselt numbers for the water tests are much higher than for mercury. This is predicted by the theoretical analysis but to a much lesser extent than actually observed (Figure 51). The difference is ascribed to the presence of turbulence especially near the cavity inlet as has been discussed above. It is also noted (Figure 31) that Saunder's data from the flat plate (9) do not quite agree with the boundary layer type analysis suggesting that an analysis without consideration of stability criteria is not suitable for natural convection flow with fluids of low Prandtl number.

CHAPTER V CONCLUSIONS In Chapter III several methods have been developed to obtain the theoretical solutions corresponding to different laminar flow regimes. Within the accuracy of the approximate methods, e.g., the integral method for the solutions of boundary layer type flow regime, the truncation of power series to obtain the solution for flow with similarity, etc., these solutions should give fairly reliable solutions to the equations. In Chapter IV, a compatible experimental program has been described. The aim of this experimental program is, of course, intended to verify the theoretical results. In spite of the fact that the experimental facility employed has been carefully designed, due to the physical difficulties in fulfilling the conditions specified by the theoretical considerations, the data thus obtained can be expected only to approximate the theoretical results. The final results of the entire study can be seen in Figure 31 where both the theoretical and experimental results are included. The theoretical solutions agree with various closely related and already known solutions fairly well. However, the experimental results from the present facility show a considerable discrepancy with the theoretical results. By comparing the results from the present investigation with the existing data, it is believed that experimental error is not the major cause of this discrepancy. From the various experimental observations, it is believed that the assumption of laminar flow does not hold in the present test range. -129 -

-130 - This is particularly true for the fluids with small Prandtl number. The effect of turbulent mixing particularly in the entrance region gives a choking effect which considerably reduces the heat transfer, giving rise to a much smaller Nusselt number than predicted by the laminar analysis. An empirical estimate of the onset of turbulent flow has also been given at the end of Chapter III, which could serve as a possible guide to the roles played by the different parameters. From the results of the present investigation, as well as the results from previous studies, it is concluded that an analysis which does not consider the turbulence and instability of fluid motion in this particular geometry would be of little worth as far as engineering practice is concerned.

NOMENCLATURE A. Dimensional Quantities (in consistent units) X the axial coordinate R the radial coordinate ~ the height of the cavity a the radius of the cavity d the diameter of the cavity U the axial velocity of the fluid motion V the radial velocity of the fluid motion K thermal conductivity K thermal diffusivity the viscosity v the kinematic viscosity p the density c the heat capacity a the volumetric thermal expansion coefficient Q the rate of heat transfer f the gravitational acceleration p the pressure T the temperature, T1 = the reservoir temperature, To = the wall temperature of the cavity, AT = To - T1 T the time elapsed on the analog computer B. Dimensionless Quantities u the axial velocity v the radial velocity -131 -

-132 - x the axial coordinate 1r the radial coordinate t the dimensionless temperature differential based on wall temperature = af(T -T )a /vxC; t1 = afa AT/vxe, At in Ra and Gr should be AT. Pr, a the Prandtl number = cp/K Gr the Grashof number = afp2x3At/[2 Ra the Rayleigh number = Gr x Pr = afa3At/vK Nu the Nusselt number = 2jTa(T1-TO)/a Re the Reynolds number = aUUp/J (the definition of Gr, Ra, Nu and Re in the text are always based on a, the radius of the cavity),,5 ythe dimensionless function of x, defined in the text i5 the boundary layer thickness = 1-j (see Section 3 of Chapter III) or the parameter appears in the Green Function (see Section 4 of Chapter III) f,g the functions of 5 in Equation (3.17) or the radial dependence of velocity and temperature in Section 4 of Chapter III I1, I2 PfyP, the functions defined in Equation (3o17) G the dimensionless volumetric flow rate at the throat of the cavity or the transformed g(r). See Chapter III, Section 4(c), C the dimensionless integration constants in Section 3 of Chapter III, or the non-homogeneous term is the momentum equation in Section 4 of Chapter III ao,al, a2 * * blb2,b3,. -,o, the coefficients of the infinite series rfo0 subscript refers to the values of functions at the wall of the cavity, or subscript used to designate a given point (e.g., o0) E quantity of possible error in -Lnuerical computation

APPENDIX I. Derivation of the Boundary Layer Equations There are several ways to derive Equations (2.1, 2.2, and 2~3). The most rigorous method is, of course, to start from the Navier-Stokes equations in cylindrical coordinates and then to examine the order of magnitude of each term. Based on boundary layers theory, the terms of lower order may be dropped so that: (1) the entire momentum equation for radial velocity, V, vanishes (2) some stress terms in the momentum equation for axial velocity, U, are neglected (3) the terms for axial conduction and energy dissipation are dropped from the energy equation. Thus, the desired equations of motion are obtained. However, the above involves much argument based on physical intuition and the whole process is very lengthy. Therefore, the alternative is to employ a simpler method as the following. Consider an annular ring of a circular disc with an infinitesimal thickness dX (Figure 1-A) The width of the ring is dR which is R distant from the origin R = O. Temperature, T, depends on R and X alone. R dR I —.-.Z-7T.......dX Figure 1-A -133 -

Under the condition of steady state and axisymmetry, the time independent conservation equations can be written for this system. Only the temperature variation of density is taken into consideration, the fluid is otherwise incompressible and uniform. With these assumptions, the equation of conservation of mass gives z (X TRdRU -dX + C (2tRd V -Vct = ~ Upon expansion of the derivatives and rearrangement, one obtains bU t bV t = (2.1) which is one of the desired equations. In a similar fashion, the energy equation (2.2) is obtained, i.e., the change of convective heat transfer plus the change of conductive heat transfer = 0 (for steady state), or a7 (wR"dIT U CTc)ds t cap Rdz Vf cT) ri (1.1A) a 3LKzrR - C- )7c1 - o Expanding the derivatives and collecting terms, one obtains U 6 VbV T '7T -7 71 T — '-) + U v 's ' +, — and by Equation (2.1) f t T _ _i 7t t T\) (2.2) + VR

It should be noticed that Equation (1.1A) does not fully describe the energy relations in the element. As mentioned before, the axial conduction term and dissipation energy terms have been neglected. Physically, this merely means that the convective heat transfer in the axial direction is of a higher order of magnitude than conductive heat transfer in the same direction and the quantity of heat generated due to viscous friction of the fluid is negligible as compared with the overall heat transfer. Bo-th ass- iilptions are reasonable and conceivable. The momentum equation in the axial direction is composed of the terms involving the time rate of change of momentum and the terms for the various forces acting on the fluid. By Newton's law, these two should be equal. Thus, Y (zTrRluR. u)dX -, - ( C-i4Rcx v fv u)dR ( R = Time rate of change of momentum in axial direction~ and -%IRdRdT f f -th rf RdRc actXi th fl Rd. t ). R = C 1 3A) the forces acting on the fluid.

-136 - In Equation (1.3A), again the gradient in the axial direction is neglected along with various stress terms. The negative sign for the external acceleration term is due to its orientation which is in the negative Xdirection. Equating (1.2A and 1.3A) and applying Equation (l.lA), the desired momentum equation is obtained, i.e., U - v R = - -f eb p " ( RI,R ) (2.3) The pressure gradient of the fluid in the radial direction can be neglected from a conventional order of magnitude boundary layer analysis. Therefore, the pressure is assumed as a function of X only and the partial differential 6p in Equation (2.3) can be now changed to an ordinary differential dp 6X dX This assumption also leads to Equation (2.4) in the text. Using the same assumptions, Equations (3.1, 3.2 and 3.3) can be derived. Thus, the assumptions of steady state and incompressibility iinpij that the net mass flow through the entire disk of a thickness dX is zero; ioe., 5J [- RdRU = 0 (1.4A) and by Equation (2.10), the dimensionless form of Equation (1.4A) is a — d — =0 (3.1) Jo If the dissipation energy and axial conduction terms are negligible, the net convective heat transfer across the entire disk is equal to the heat

-137 - conduction through the wall, i.e., b 2R d Uf CT rj (1.5A) 0 Again by Equation (2.10), Equation (1.5A) becomes x f tdr = (>a b ) 'l (3.2) The derivation of the momentum equation in integral form is a little more involved, however. First, consider the time rate of change of momentum. Since the dependence on the radial velocity vanishes in the integral as shown by Equation (1.4A), this inertia term is merely Io (aRd Ru-n) There are at least three force terms; the external acceleration, the axial pressure gradient and the shearing force at the wall (the rest of the forces being neglected). These three forces are, respectively, Z- fR I aRRdR.d d a t 4-) R a C X (1.6A) o o 0 0 All the integrations are, naturally, with respect to R. Equation (1.6A) can now be rewritten as Forces =-Z1RdR (f+) ) 4R )d d (-l.A) 0 Recall that previously the pressure gradient has been determined in terms of stress at the wall, where velocity terms are zero. Thus +-I I i a(v U +I- - (2.6) 0 oR R R R

-138 - and the density variation with respect to temperature follows 7J 4 C ~ T - 7o3 (2.5) Combining Equations (1.7A, 2.5 and 2.6) with the inertial term above the integral form of the momentum equation is obtained. i.e., ca fo(2RrUm )urx =-r.)Rf(f the mm)tux.+q(a)o 0-S R+-L)- (V. \c_ Thus, in dimensional form, the momentum equation is > | R D> dR = -;if ( T- To) R2 d R -b Z ( -e(R) Changing the variables into dimensionless form by Equation (2.10), the final form of the momentum equation is 0 '0 /L LJ d ii- rL. t (I ri- + u, (3-3)~ T' ~ ~ ~ ~ ~ I -Z b (r iO2 I~~~-I

II. Detail Derivation of the Momentum and Energy Equations from the Boundary Layer Equations for Numerical Computation It has been shown previously that the equations of motion in integral form are, r L d r=o (3.1) AfL*(3 drz. = ( )(3.2) e =1 fL443 (3.8) A= \ rc L Lk xt^^ )?C rL '<^^ // L^( and the velocity and temperature profiles are assumed to be: f - (3.8) It is now desired to show the detail derivation of Equations (3.43 and 3.44) from these equations. It should be noticed here that three unknown functions, 8,y, and p are present in Equations (3.8 and 3.9), all of them being functions of x only. Equations (3.1, 3.2 and 3.3) are mutually -139 -

-140 - independent. Therefore, they should suffice to determine the unknown functions uniquely. It can further be seen from Equations (3.8 and 3.9), that r = p is the position of the adjoint of velocity and temperature functions. Substituting Equation (3.9) into Equation (3.1), one obtains j P- Lj. - -il -' ' j( P3)' f ((nL-n^/I)AU =, P'I E= (2.1A) Let I n~~-^^. =,-p) - Equation (2.1A) thus becomes >0)+*3X+;y 20 _(1-0) | (I t Sp - d-p ) -^)+ I f, -Fd, ) - 2 - (3t )- -- ) -3 - C-p)(32f), S _ 5(3.t2(~') (3 10) CI-~~C3+2r)

-141 - It follows upon substitution of Equations (3.8) and (3.9) into the momentum Equation (3-7) that; o S LZzd r = (L A + (L IAd iv ~0 o JjrxLt; d L- It - zZ AI I+ ) 7. - Lf. PLj + j{-2t^( ) N-l)Ul a(P t- )l2(N- )jldh P =-1t - fZt f-.( 2_ r ) [ -, ] rL-I ( PI- 'B c —, )l] - I-',I.~- = fl. 1i4^ 'L 2. 4i-[rl~~-*^ (f~ {,-r) T '2 1 I 'Tl - ) t [ -P ) E id r Again let I-P drL= (I-P)dt /rL- = (C-p ) ^-P l = \-P)(- >) r., - p,; -= / > O o 0 - - jt W j(o-uPr) ip j'^{ y\f-P)(-J c(2. A)

-142 - The integrand in Equation (2.2A) is split into two parts and carried out separately, i.e. c*)c ji (,, r 1a' l+ S 1,t)~ ( crI ) 4 7 -Ii =i, zc\-r21P- t 7 ((I-P) >-, dd o 0 -F0 () ) - tI = -,-e c -p) t-f o I-~T(,-i)- If 5,-P)( - o =-3 l-P) 1 5+5t-T —) 3-3f1'+^ ) = -L l~-p) T '5 C3+p-y(\- P) ( 3-3 p+ 7 = (\-^[-G- l3+p)+3o ^o-2 2

-143 - Recall 5(35rzft (3 ) (3+lf-(\FPy2 I- ^,(,4() 5( 3+L(B4 45A)(l -tIoJ3) 5 (3+23+(z\) (5+3~)? = - ~o I ~ R) ( 8 st 4 h _ 30 1 O llo (3t52t) S4o (3C^^2(1- ) or p _-5P(1^ (5+zl)(3t2F ) -2O(3>tS zkl r) (1 II~ ID r) (I-1) t3t125) s(3t-2. + )Z( 5-3 ) (\-P P;-P)i- 4 -5=( \- 5)Et st1)L3+2) = =s-' -5 - 4 ) ) - (4 5 -; 1-3 '-1 7 - t i;? 4-rIt 5) =1 - srot?4 l4>+35lz2=t 13 —i58 4 -S448/35 -3o 33+2.p>t r )-(,8to) (\-g,,~3+ (p) = -20 t(l+ lllr-l5s3X-plo5p3-'-zo - ~ ') = - 148O - 4bD 3 t 5ce r-+ 1oo r + 1 44c r 3t 4 c35 = 1 (Sft?^)(^- ~-^10 k\ 44 r, r4 -I 2- -+5 p q 7 -t 3 1 5o l t 125 Bo f - 4- p S +-qt P Substituting these quantities into the definition of P, it is obtained.

-144 - I,7 5. - I -1 ( I c -t 4- 4 - r -> -t * i - P - - 7 f ll oC3 (V O-i f, I r /- < 2 t 3,,. 35~~5$ 38713+ k40r t 3,/1 / + - C, I r i j, 'Z - -~~ ~, 15 4 '~~~'a I ~~~ -f C- +0 z A;. "', I~~~~ I- -L tY ) IZ~ u — 3 p7 f + (-0t 2- 434d,2- 3 f a 7 0 4 — 2 /1 r-Ip ' 4-D (It50-^ le -4 0 ) The momentum equation can now be written as 0 cr (I) - - t.c-l I I. T - dY I x k' " (Ba ) bJ a *aP [L- = - -Z f~-P Np-+T O - )-L - -P L i _, - { ' I _. y -9^ -.) - (43 +- ^ 95 (3+ -2 1 )- - L_^zpO~l-f^ j C-512.)( I- f ) I \ dr LZ Ax I d

-145 - and 0 \ - S)-i1o-._ / -.. P.4 ~-lf.t)-P_) ~-lf)t (tJ~-), -~ =c jl(-P~i~4|~t1-P)550t- 27( f?)(. -1) t 'IV- -3'" K)' KJd) 2. Collecting the appropriate terms gives, I )+4* 7T + -(,3t ) ) 3L 4) (+ t 2 3 j c\-^)^ | \^)(-P) +- 2S\ ()-P,)l(5+5l3)-,^ t 32o-p^ Ib (5+0) + S -P

I —,. r rl, lr I- 4 _ -~ t4 -.H'Ij II 11 I) C c 1 (I r4 c - (-I ^I + lcd r. rf 4 | $ I - -_ V -1 t 1 -- I1 -t -,4 4- 4< -1 -1 --I-1d r 1 --- II r4 -1 N 11 CA rs CO_ (41 t l-t0,I t~ t.') fri C(34 ri N - C) rif I — I-) Co? f. 4 rI rl C /-\ CL00 O0 '+s rl nl 00 rc^0 tl r* tO 1I 0 rU sV 5i l^ 4 - 0 -ri1 0 6o in rc r(Ct fr< rt ln N r< 0 N <Cf -tcl_ r I dj~ +P iS /-9 C-.A r^ 1-1 r\ (VI c2 -(ri 4 t0 1 4~-; I I 11 - 11 II I -I >~ I -1I I I It r IZ-1F rbz -Ib 16 CI cr -( II '4 o/ -Icl II f II II C-11~~

-147 - and finally, _ 3) { C, 1 4 5-t i o~ P+4-1 1 3t z 1-62 2: 4 3 I 2 TX( 84-o0 __+zpZ_ _ () *( 1c - - y, y (3t 43t3 (2.3A) 1 l 1+ ) ()-y3 ) Likewise, the energy equation is utilized in a similar fashion, i.e., -6- p -I t (-) (1-) =? _,V, 4^+ a -3. ["a(i-)(-i f i-1)tl 0 I - i r -^^^ t ^-\)(-r - rw(^ Cl-pc ( -Yb-p ) -^t -g - P -01(i-^'t3 = {.(hi ) -p)( (h-V a * ^-1)^^ a; i

rL" I '41 m', -— ' -- — ' ---- 0 Q1_,- - ' r-' - - 1 * ~-v N I- rJ.^ c^_ c^ <"',-i~~, I t i- d coH~~~c t_ U- Q2 c - i - L - -I Ct /^+~~~~~~~~~~~~~~~~1 ri- ^ ct^-^'? ^r~~~~~L 'f- C O- ' - '|rt ^ c- %-) () CA.4-t- I -. I r II cI> 11 II II ii II- ' I+I /,-, |..- rNl —. L I i L c^- c^ <^- ^- ~ 0 0 ^|^> ------ ^1^~~~~~n~~~~~~~ 3^~~~ z~'. -t c ^^ro I ~ ~ ~ ~ ~ ~ d d -I 02_.._ -I r/ 1 0 0 >~~~~~~~~~~~~~~~~ II II II II II II II II II I ii ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~I.-.-..

-149 - 'I?) ~t - = [- r-^> bn/L =l -kL kk Equating these two parts together, one obtains 3L T 4S A^^!^^1 g l: - = (2.4A) by l 4l34o -A) 1_ Since the r-variation is taken care of, the final Equations (2.3A) and (2.4A) are ordinary differential equations, namely A.T -1454o p(.2 )fez - 1 1 i is O-p), _y3(,t4^^5+:.-) (3.11) d - t 45 s 1^P-tlm i n 67 >p3 -, 1 J _(3.12) 1 ty \ -yc3^^ J 1-0 Subject to the conditions

-150 - which are obtained from Equations (2.14), (3.8) and (3.9). From Equation (3.11) 5 t CpT )( I 4' Le t ^y(3-4(2+3p>(3 Let -I4 q 5 & i r 4 7 1 t4-21 -., j t 4' 1 f 2- -1 485 n (O I- ) 47 - i(,-^, );+Z'^'z ^ )+2l ( i-2) 7-( 2t-9 4o0(32-. o f.1 9 )::4 oo - 6 5s5o t~ 48~0o -3c I- S t - 44 t L p - I ~ ~ ^ o (254-2W'- 4 53 7o ^C5-^)3 Y(3.-) I '?, 1 7 8r4t, = 4 ~0- oq"S0 -1 G43 t4 -y3 (J1^ ~ ~~~~ ~.. r -Z.4 ~,

-151 - I, = Z —t (4 — ) =- F;L = 2. 3( o-1o%+3-S ) ( 5- z -%) 'I L ' I go Therefore, Equations (3.11) and (3.12) can be written as ky dXI =B'^ G- rA id?F d v I E, ~ (J., d< - x H? 2 (: ti+,PI -YPZ C-) tl51, 11-T, ^ i) or as, d1 cl1-.11 "14 = (j l iI.) J r+ t ) X -> 1

-- ro') I NJ LC\ H H >0 -I VA!< -4 -Oz n- -C —, y)f rl-~ >0 r<k I ri >Q -t >0 >0 I >0 >0 I yp -r 4:s -1 l 173 — l C (-j N ')O P4 II -6 (y-L qln3 (NZ 4-., (N.H 4-o (q) I b A >0 9fJ It p\r -4K r>0 Q/,/ 4iD Dt~ls~,4_, CN-`C — >If q (>< >I c1 -rI D3 >0 -Vb I-0 -I II >9=1 <s^= f\ -N 4 --I rW~ b I -RT lo 11 r\ Ir~ cA C>0 1 I-E b (N >0 II lu T -5 II -- u

IIIo Brief Review of Lighthill's Solution for Boundary Layer Flow Regime(25) If the Prandtl number of the fluid is very large, then o-1 approaches zero. With this condition, Equation (3.7) is reduced to an algebraic form, ie., | A = I ( L e t ( L) (301A) Now from Equation (3.1) and Equation (3.9), 5 is solved for as previously: (2.3A) left s or (c1 ) (I (3.1 Due to the simplification shown in Equation (3.1A), Equation is reduced to the solution of 7 in terms of P. Equating the ide terms of Equation (2.3A) to zero, one obtains t, 3(%+4-ft3F) O o) y (3 ) (_3_) -_(l-) (3.2A) v (, 3+4[t3 3^) " Equation (3.2A), combined with Equation (3.11), gives the energy equation for large Prandtl number. td (-I3i(3+-p )(45tI3 -Z-)l al2-3 I7 I d l 3oI(+ 4r>3 ) J I (303A) or if 3 p~[() - - ai ~(\3tZ) (4.'3 iFs 131 F - p then Equation 3A) can be writn then Equation (3.3A) can be written as 3d-Ft ) L= 1I (3.4A) -153 -

-15 4 - When the right side of Equation (3o4A) is integrated from x = 0 to x = 1, the limits of integration of the left side are from = 1 to l = 1, where p1 is the value of P at the orifice. Therefore p, I= t, O -p F() (3.5A) Now, if tl is given one can solve for p1. Recall the average Nusselt number is given by D~= ili-i.^ LN [(2.16) and by Equation (3.8) and Equation (3.4A) Nsu-= Z\ T e s tJ aDee d= -tF(P ) fI (306A) 0 lI Thus, the final solution of this problem with large Prandtl number depends on the integration of Equation (3.5A). It does not appear, however, that an analytical expression exists for this integrand. The numerical solution was carried out by Lighthill in his original paper. Later, the author programmed this integrand and obtained the same results which can be seen in Figure 4 of the text. It can be seen that when t1 is large (and consequently p is small) Equation (3.2A) converges to Y t 3 (3.7A) a result which can be obtained from Equations (3.35) and (3o36)o This is, of course, only natural since Lighthill's solution is just a special case of the general solution shown before.

-155 -Since from Equation (3.9), 7 is the dimensionless core velocity, the dimensionless volumetric flow rate at the orifice is proportional to 2 y, which by Equation (3.2A) is 0(3+)(3+ 23)( 1 -p)3 s(3c+4 p t 3(3 ) This function has a relative maximum when 0 = 0.38 as mentioned before.

IV. Details of the Similarity Solutions It has been shown previously that with the assumption of sufficiently large Prandtl number, the equations of motion (2.21), (2.22) and (2.23) possess a similarity transformation and the resulting equations are g '+h(f+iY- -f O= o (4.10) "'- L" + =(- C) L (4.11) and the boundary conditions are (+) = f') =() (o) =o0 (4.1A) (g) =o;3 (o)=t t (4.2A) I ' V f l 1 9 (4.9) It has also been proposed to solve these equations by a combined Green's Function and power series method. In this section these details are carried outo The Green's function method to solve a linear ordinary differential equation involves to solve the homogeneous part of the given equation, to construct the appropriate Green's function and, finally, the integration of the Green's function over the given interval of the differential equation. It appears that between Equation (4,10) and Equation (4.11), the latter one is more suitable to Green's function solution since its homogeneous part is linear and boundary conditions in Equation (4.1A) are all homogeneous. -156 -

-157 - In fact, the solution of the homogeneous part of Equation (4.10) f --. t t =0 is readily obtainable. It is of Euler's type and the general solution is -() = CI ( Cz + C3 /L) /L (4.3A) the Green's function accordingly assumes the form e ar al.o qa(Ce i )Ea o (44A (4.4A) l b, I bt b, Quar ) < c -X I There are altogether six unknown quantities in Equation (4o4A)o To determine these quantities, three conditions in Equation (4.1A) are utilized and the rest of the necessary conditions are from the basic properties of the Green's function. The Green's function for a linear ordinary differential equation of third order should possess the following three properties; (i) G(S;~) is continuous over the entire interval 0 ~ g ~1l (ii) Gr(S;S) is continuous over the entire interval 0 - <, 1. (iii) Grr(r;S) has a finite discontinuity at r = ~. This discontinuity is determined by the coefficient of the highest derivative of the given equation, i.eo, jCj) - iC tt ) =-_ (4.5A) Therefore, these six independent conditions suffice to determine the Green's function in Equation (4.4A) uniquely.

-158 - Note that [o -O 0.no A L-~O FrLL t-r L VO fL — 0,4 fL p-nf L —>O -M O. y%2 =0-, I\ Tl for all n> 0o Thus, by f(o) =; f'1) = o=, cl, =0 b, = bz and by.;. (R'; ) _ ~ Ca \d3I1,I-x.) ). C (-) +[L'-,+nRx L) bZ o <-LC< By (i), (ii), (iii), respectively a, t A(3 P\ ) = (-\ + It ^ - I Q ) b2 cz tt +)ct-2I-)C = 4( t~-) b -k7 3, " + N-t )C -4 ( +l+ -%) 6 = — V -S (4.6A) Applying Kramer's Rule, the linear system (4o6A) is solved. The results are 41=4O c = (- z C3 - _ ~t 4-~. I b,= -4 k'=_

-159 - Thus, the desired Green's function is [7(f={ L-% t t 5t;t \7 t(4.15) It can be seen that when r = 0, the first equation in the set (4.15) of is identically zero; so is the second when r = 1. In a rather straight forward but tedious fashion, it can be verified that Equation (4.15) does satisfy all the prescribed conditions, as does the homogeneous part of Equation (4.11). Now the solution of Equation (4.11) can be obtained in terms of the unknown function g(r), and the constant C, i.e., f(A,-) = [ ^ (^ [,X) C] (4.7A) 0o After a lengthyintegration, Equation (4.7A) is found to be |2 ~L Ar At s| H t -AJ it s L 3 " { ()0 l I It^ ^ f ZA t i 1 -I Z CfO+2 A -A ') (4.17) i nAs mentioned before, the unknown constant, C, is now to be determined by Equation (4.9). In doing this Liebnitz Rule of differentiation is invoked. As the result of this elimination, Equation (4.20) is obtained, i.e., \ i | 0 L'ta SI' ( I atoi i NA' + 'd tz.i%4p d J (4.20) o 0 0

-16o Now, the solution of function g(r) is assumed in power series form, s0 (N') = X/: A -O(4.21) such that g(r) satisfies simultaneously the conditions in Equation (4.2A) By employing following integrals g 2,ll' = ---- -t + Cl and )?^T 0~ the final solution of Equation (4.11) is i) ~ t | d L -)n i -;4Z [ )? A ) ( I - 3 (4.22) and for the reason discussed before,. g(r) is an even function, the even solution of f(r). is thus fX > — 7[(itl2- - - ) (4.23) It is of interest to varify that (4.23) is the correct solution of Equation (4,11) since the method of obtaining this form of solution is so indirect, and the steps of computation so involved. Equation (4.23) has to satisfy the following conditions: (i) fic)= fli)=f(o) =o (4.1A) (ii) t iL'8_R fAtfi=^ C- ) I' (4.11) (iii) Clf'CO-rF'^j) 5ince t~,)=o bk (4.2A) and (iv) 1 f'(r) be bounded at r = 0. r

Since O I I, [o h,00 0 ( L l 1 < 2 % a f}^I O Co I, L "L)( l)r I 'Ot )-O - ~J " ~''''= 8, = o ( +I, c) L (- )( I-4 ) - ( -V)n- 4)7 =o Thus, the conditions in (i) are all satisfied. Define ~ ( lt-l ).(I2-z) I. then f(N)= Z S i, ult-+ h — _ )(Cr -.) -( -/u ) f (A)=2 s, lc1)(iiU2R-4z )- 2A- l+4 I S0 - \Z -F.~,~, rll1111n.-1i?~~ 2 ) 2 — + 3 )2- 2n+3-p. )J II ) I 2 5 [ ( ) I ( Z 4 )(2 ) (2 -- ) -./L The left side of Equation (4o11) is, by substitution, The left side of Equation (4o11) is, by substitution, v <r3 xr1+3 2 n -Z4(D — ) -t (Xt4) (' + 3) -t. ) /L -( 4 ) ( 2 R-\2 3)) t n- (21-Al 4)( 2 1, ) A2L'3 + (+>)C,(? 4~ r -.[ zt- (2'-. 4)/ 'L2M3 ' }

-162 - 1= z \154 " - ( CS2I1 [ 3;C I 5 — 1. (Y+))/L t 8 (h I) (ot ) / L I CvA t l) h 21+1+),-, I +1 1~: 1i2 1 3 - /2 (i)4i)(-t2) IL h1o * 3,() = n-((fL> n^o 2.C Lnc - tI))-r-) (4.8A) Since by Equation (4o19), 1 = D U (20u 1, Y - 4- c ( - ) J- -1.A'Z '= 4 L u1 ' 2 -Q~l - = 7 / -.. ---. c,~ ' )~ii,~. (4.9A) Substituting Equation (4,9A) into (4o8A), Equation (4,11) is obtained. Thus the item (ii) is verified. Now -"(1 )- f"'() -9n -3zD1 -4-n g(ht l) CVt -) Q2n I)0.o CtL-) ( A I ) ( o +l +) t0o h=-_D -( ti l ) (1 ti) + 2...t ) (t,,,) l - 2.?2. 1I —).. j. _t- (c i \I)(C+2L) 7'2Ju

-163 - Since.to ^cK)-o b Ca<. -2 m O u f ^ C= + (n-f () 1t, ta,,,tt) and item (iii) is verified. Item (iv) is easily verified since -f'( S) n t% 2,t ) 2 J) t_ ' ( _S_ rS(l-n)(l-z \ ))-2 >t(2,rtLK[)l2 - ZJ 7 21 L Y yt=O h If the series of g(r) is convergent (which is assumed) this series is also convergent.

V. The Details of Digital and Analog Computer Program

-165 - FLOW DIAGRAM FOR RUNGE-KUTTA FIFTH ORDER PROGRAM 1. Correspondence between the notations used in machine and in the text: machine text machine text X,U PF Pf Y,V 7 PFP Pf' Z x G g (Ref. NU Nu GP g' Eq. (3.17) PR a PI1 P11 in text) GR tI PI2 P12 2. EQ.(1): - 120 1 -- (3.29) 120(1+ 1.05a) EQ. (2): x - -7atl 4.. - From Eq. (3.30) 1600(1 + 1. 05 ) EQ.(3): Nu- 70tl t3 — From Eq. (3.33) 600(1+ 1.05a) 3. Al 0, B1 0 A2 H 1/3, B2 J(l)/ A3 2H/5, B3 (6(2) + 4J(1))/25 A4 H, B4 (15J(3) - 12J(2) + ())/4 A5 2H/3, B5 (8J(4) - 50J(3) + 90J(2) + 6J(1))/85 A6 4H/5, B6 (8J(4) + 10J(3) + 36J(2) + 6J(1))/75 These are known as Runge-Kutta Equations 4. Eq. (3.45) and Eq. (46) are written as d - 3 _ NA + NB + NC. V dg DA + DB + DC U dx = V2(G. PFP - 2.U. PF GP).V2 dt DC + DA + DB The terms NA, NB,... etc are corresponding to those appearing in the text. 5. J(N) H * d K(N) - H dx dt M(N) - 2 K(N) U Those are known as Runge-Kutta correction terms. 6. YH - (23J(1) + 125J(3) - 8a3(5) + 125J(6))/192 ZH - (23K(1) + 125K(3) - 81K(5) + 125K(6))/192 HUH - (23M(1) + 125M(3) - 81M(5) + 125M(6))/192 Flow Diagram of the Digital Computer Program.

-166 - co'ltMPILE MlAD, E::.:ECUTE, PUi-1CH.- lEC' TM..ET, tF -----—,L --- —-----— fi'G_ _ _ER-N _______ --- —------— _. __.___________ -— _________________.______ I HTEGER rt _____REPO FEPE' FQR1PT IN.GR5PR5U ____________________________________ DI[MENHSION J,:,::,':::' L C 6::, r,1,::, ) READ READ FORMAIT INGRPR, _1 V.'ECTOR 'tLUES IN = i F12. 2,2F12.4*$ PF F 2400. +U* (:-6550. +U* 48L. +l:*-1..: 3'C0+U*: 1 44. -9'. *I1.', ):,.: __ --- —------------ - --- ------------------- --- 1:28':. * C25. +U* C-20. +4. *U)-l: PF'F - ' -3000.+UII* (.3-00. +I'I11 5-4 r+l.I*, -. +!, (,'If540.=I +U. 1-' 7.4. + ___G - (420. +!:-6i80.i +_iI,*:.3_,,7. -f,2. *U:":,:-C 40L. *::'5 -2.-:,:-_1:,.5_ ) F' 1 4= s4. - ). *12. F 3 = '- -.P +GR*PI G + iU: -;'. *FF+PFP+3. **G* - 2: ------—. *-PR*GRPFI 1 - GU --- * i. —:.. -. G. *F' F+F'FF' +..F.3FR'.. + F'P*I 2::, -:CU+ *IIU:,* 9. *PR*GR*P I 1 *G:G~)..:. *:..PF'F+FFP+3_. *FR*G*p2):. ______________________________1=C9.3_____U*_*_______..P ___*________-__* _,F+PFF,+.F3.*PRp r.. — -~,::,.': U e"-U e Ut ) +,::.. * F'-F.-*,F:_.F; I- -:'-.'.::.'.'.'4. '.- -+'G, u F' F. +-F:F-F. + 3 +- F'R + i' - P- [ 2:,:, 1.*: C,-:-:-.,-. 3 H = 0.001 X= U PRINTt F ORM IfT T ITLE PR GRS, Z:.::, ':',,.. -::.' + 1,.' '*. ' E --- —----- TOR Y._'RLUES TITLE = $i 1 H1,-i 0, i45HSOiLUTIONS FOR -BOINDARY LYEL'..E 1F.R T-'-F' —'E FLOW.. REGIME/1 H O, S1 0, 4HPR =FC.. '4, S:8, 4HGR =E1 0. 4 I.1 HO Si1,:' 21 HIiEN2INLES, 12, 1 3HD IMESEIIONLESS, S1E 1 13HD I MEN'IONLE3lSS 3S1 2, 13HOIMENSIONLESS./1H,S 15, 1 4HA.:IAL-DISTANCE,- 1 ';, 14HB.L. TH 4ICK t iE- SS14, 3H.ELOCIT' Y, S 4 1 4,14HCORE FLOW RATE 1. HO.1, SS22,F7.4, -5S 1 8-, F7.4 2Es25.7*; START THRlOUGH ALPHA, FOIR N1, 1, N. 13.. TFNRtS-FER TO L CNH' L(I), A i 0 B = TRRNF-:ER TO r:ALC L:C2) P H.3. B -..'1':,.3. TRANS-FER TO CALC L _LC3) A -'. h * 5. B - = k t.'.. 1 +;*J 1 ) 25. TRHNS'-FER TO CALC B =,..1 4.'- 6 1...,.:" +- '1.J.::'2* +4..,:: 1 '.'.:...14. TRF-.I-kNFiF.F Ti CLCL LC5) A *H 3.H, - = 0 '': (3)-'. +"0.,:' ~. +.'.J, '1:., )',81 TF.RHNSFER TO CRLC L.(6:) H 4 *H...-..5 B _.. - '-',::::. *J4 + 1 C 3) +36S. *.J, 2-). +. *J:" 1:,.. L C IU.;:+R F =:2'400.+..I*C —,550.+U*C4::880.L +LI*,-1 3i,-,+iJ*,:C 44.-9. U.)'./' 1,:' 280. (C25. +1_. *-20. +4. 1_'.:, ':,- - PFP -=. -30010. +L..I* C3i;-.00 +ULI* -4i50+U*:-9.0. +LI* C5410. +U* -1 1 7. + 1 + -. I _1:.. -' 70. +:125 + -1.,:'- I. + 4. *1 *._ (5. * -2. _., ) 1G =- 42 0. +ULI* C-':,0. +I*C 3.e17. -+.2. * U ) ):). / C/8 40. * C5. -2. *1_I) GF____' C-2560. +UI*,..3670. +_*'-1 '-,664. +248. *U::'::,)__ 1.......840.*~,:5. -2. * U::.C5. -2. * U) — ~F --- —-----,l 1.':4.-U)/1i2. F'I2 = C3.-1 *::10. +U*C-1 O.+3.*U:U)))/5.-2. U Program of Runge-Kutta Fifth Order Method for the Solutions in the Boundary Layer Type Flow Regime

-167 - NA F' R*G3R*GP*PF'I 1 *U*UI*U*U " H B == F-V' R * GPPPI2*L I.. _______ i -i = FF......F...................__................. l B - - V* F' R * G PI F' I ' _ I.I D F' RGR*Gi*F'I 1 *u li11*U.................... -..G*PI2..... rD - -2. *. F'F ~.-'H..-*. (5*F...,'"i =i" * F- *... F'::.::.,:(F+ GDC +:Df i+DBE::.....................,:: r::, - 2.. C. 1.. i:,*::K C r.H)..H..... I........ H..':.: -. ',.,..: ' '. iJ.: Ht.. ': C: + D C i + [1 E:.. '+~ ' = C+H 2.".-..JCJ _. —..,..-..=-.y+<23, *J<1^ 2.5..*J. I 2 C3..1..J.5. INU +z~ ~ ~~~~~~~=:t..Y C i ';" - Zt +: 2, 1+12 *1 j ' t 1 ' + 1 22. * *:,:.:- 1.. " '-' ': ~ -.- - -, -- I 0........... -.........-iti.......'''' ', f i iRHEHEFEF.T Z, GI1,0 P.. 5 N-. U-L 1 C'.1 P.. 7.5 EI ID IF CONDITIO AL....,.,.,.,...,. H, -. ' E..Ei P I T FIF.. UTi. Z,.' ',:1 -. -..,.:*, ' _~___yECTOF.: Vi:.iLLES. OUT. —:.,,1H. '22. F7.. F7 4 ' -"2:'"' I.,HEHIEVE.ER Z.L. 1.0 j TF. ANiSFEF. T ST —,RT FP RI.T FORF,:H T HEFtTTC..,l,.PR...-iU ''',,! E: T F.:,.,i P-t L i 5 E '-. I _.- T T -: -. i H..,1. H R i:-i I,:" F-: 1' 'I t E: E R -'F 7'. t ' 1H.1 t i'-:,Hr:F4 —'NDTL 'F::ME'R = F7 '. I 1 6H1;, ',s '' E L - '. r. i:iT i;,. Fi j.,r,-::',T.L-" R E A:. '. _ _ _~.... ED CF- P IR ARR

4EMII -?dF dF d Ito -.._~ -/i. -' r j eO 10 oe — e- C-lfA C u/./o The Circuit for Analog Computer Program, )-oo 0 7Q J To X-AVy/s Of -/'00 / -<a 8irr DOO

VI. The Physical Properties of Mercury as Functions of Temperature -16 9 -

0.36 0.34 I LL () LL LU. 25 -i w x I — 0.32 0.30 0.28 0.26 I~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ I~ F~~~~~~~~~~~..... ///~~~~~~~ /'O /~~ /.O/ / Z/ zzzzzzzzz4zzzrrrrr, I I 0.24 0.22 0.2n W. f- 0 40 80 120 160 200 240 280 320 3 TEMPERATURE, OF Figure 2A. Thermal Diffusivity of Mercury, K, Ft2/Hr. vs. Temperature, OF. 60 400

THERMAL CONDUCTIVITY, BTU/FTa-HR- (F/FT) 0 (D H (~~~0 - i i i H' ___ _ (D, C) 0 0 C+ i9y 1-Jo._ o '0 0 0 m @~~~~~~~~~~~~~~~~~~~~~ I C) ~ ~ ~ ~ ~ ~ ~ ~ ~ X c —i to " os, m __ _ _ __ _ __ ______ g 0H ~ - _ _ H Cl- - 1 CO l-t 0 0 04- H ____ ___ (D o c+ ~0 C (0i

0.025 0.020 b z 0315 I 0.015 _ _ _ X 0.010 -. 0 100 200 300 400 TEMPERATURE OF Figure 4A. Prandtl Number of Mercury vs Temperature, OF.

13.58 13'54 13.50 _ 13.46 U 13.38 13.34 13.30 13.26 13.22 50 70 90 110 130 150 170 190 210 230 250 270 290 TEMPERATURE, OF Figure 5A. Density of Mercury, p, Gm/cc. vs Temperature, OF.

5.0 J 4.0 >3.0 Z 2. 0 _00 200 __00 400 500 TEMPERATURE, ~F Figure 6A. Kinematic Viscosity of Mercury, v, Ft2/Hr. vs Temperature, OF. I -- I

3.34 3.32 0 x LL, 3.30 0D -a m 3.28 t4 w I o 3.26 w 0a cn 3.24 3.22 3.20 - ----------- - \ --- —--— ~N -- ^- ------ I \ --- —\ I --- Y --- ----- -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - - \ --- —- — ~~~~~~~~~~~'", _ _ _ _ _ _ _ _ ___ — \,^ _ —~~~~~~~~~\, I Il 0 40 80 120 160 200 240 280 320 TEMPERATURE, OF Figure 7A. Heat Capacity of Mercury, c, Btu/lb - of vs Temperature, ~F. 360 400

BIBLIOGRAPHY 1. General Discussion on Heat Transfer, Seco IV, Inst. of Mech. Engineering, 1951. 2. Hammitt, F. G. Free Convection Heat Transfer in a Closed Cylindrical Tube with Internal Heat Generation. Ph.D. Thesis, Nuclear Engineering Dept., The University of Michigan, 1957. (Also available as Industry Program Report IP-259, College of Engineering, The University of Michigan.) 3o Smirnov, A. G. "Free Heat Convection of Mercury in Closed Cylindrical Tube with Internal Heat Generation." Zh. tekh. Fiz., 27, No. 10, (October, 1957) 2373-80. (In Russian) English translation in Soviet Physics - Technical Physics, 2, No. 10, (October, 1957) 2209-15. 4. Hartnett, J. P. and Welsh, W. E. Experimental Studies of Free Convection Heat Transfer in a VerticalTube with Uniform Wall Flux. ASME Paper 56-A-113, 1956. 5. Siegel, R. and Norris, R. H. Tests of Free Convection in a Partially Enclosed Space Between Two Heated Vertical Plates. ASME Paper 56-SA-5, 1956. 6. Lorenz, Lo "The Conductivity of Metals for Heat and Electricity." An;ialen d. Physik, u. Chemie, 13, (1881) 582. 7. Pohlhausen, E. Z. angew. Math. Mech., 1, (1921) 115. 8. Schmidt, E. and Beckmann,W. Tech. Thermodynamics, 1, (1930) 341-9 and 391-406. 9. Saunders, 0. A. "Natural Convection in Liquids." Proc. Roy. Soc. (London), Series A, 172, No. 948, (July, 1939) 55-71. 10o Goldstein, S. Modern Developments in Fluid Dynamics. 2, Oxford Press, 1938. 11. Ostrach, S. An Analysis of Laminar Free-Convection Flow and Heat Transfer About a Flat Plate Parallel to the Direction of the Generating Body Force. NACA TN 2635, February, 1952. 12. Sparrow, Eo M. and Gregg, Jo L. Details of Exact Low Prandtl Number Boundary-Layer Solutions for Forced and Free Convection. NASA Memo 2-27-59E, February, 1959 13. Hellums, J. Do Finite Difference Computation of Natural Convection Heat Transfer. Ph.Do Thesis, Ch. Engr. Dept., The University of Michigan, 1960o (Also available as Industry Program Report IP-461, College of Engineering, The University of Michigan.)

-177 - 14. Braun, W. H. and Heighway, J. E. An Integral Method for Natural Convection Flow at High and Low:Prandtl Numbers. NASA TN-D-292, 1960. 15. Herman, R. "Warmeubergang bei freier Konvektion." Physik Zeitshr., 33, (1932) 425. 16. Martini, W. R. and Churchill, S. W. "Natural Convection Inside a Horizontal Cylinder." J. Am. Inst. of Chem, Engr., 6, (1960) 251. 17. Jodlbauer, K. Forsch. Gebiete Ingenieurw., 4, (1933) 157-172. 18. Sparrow, E. M. and Gregg, J. L. "Laminar Free Convection Heat Transfer from Vertical Cylinder." Trans. ASME, 78, (1956) 1823. 19. Polhlausen, K. "The Laminar Free-Convective Heat Transfer from the Outside Surface of a Vertical Circular Cylinder." J. of Aero. Sci., (June, 1958) 357-360. 20. Yang, K. T. Possible Similarity Solutions for Laminar Free Convection on Vertical Plates and Cylinders. ASME Paper 59-A-86, 1959. 21. Jakob, M. Heat Transfer. 1st ed., 1, New York: Wiley and Sons, (1949) 532-542. 22. Batchelor, G. K. "Heat Transfer by Free Convection Across a Closed Cavity Between Vertical Boundaries at Different Temperatures." Quart. App. Math., 12, No. 3, October, 1954. 23. Mull, W. and Reiher, H. Gesundh.-Ing., 28,(1), (1930) 1-26. 24. Poots, G. "Heat Transfer by Laminar Free Convection in Enclosed Plane Gas Layers." Quart. J. Mecho and Appo Math., 11, Part 3, 1958. 25. Lighthill, M. J. "Theoretical Consideration on Free Convection in Tubes." Quart. Jo Mecho and App. Math., 6, Part 4, 1953. 26. Leslie, F. M. "Free Convection in Tilted Open Thermosyphon." J. of Fluid Mech., 7, Part 1, (January, 1960) 115. 27. Liu, V. C. On Laminar Free Convection Flows in Cavities. Abstract published in Bulletin Am. Phys. Soc. Paper read before the Divisional Meeting of the Division of Fluid Dynamics held in Ann Arbor, Michigan, November, 1959. 28. Ostrach, S. On the Stagnation of Natural Convection Flows in ClosedEnd Tubes. ASME Paper 57-SA-2, 1957. 29. Martin, B. W. and Cohen, H. "Heat Transfer by Free Convection in an Open Thermosyphon Tube." J. of Appo Physics (London), 5, March, 1954.

-178 - 30. Martin, B. W. "Free Convection in an Open Thermosyphon with Special Reference to Turbulent Flowo" Proceedings of the Royal Physics, 5, (March, 1957) 91-95. 31. Martin, B. W. and Cresswell, D. J. "Influence of Coriolis Forces on Heat Transfer in an Open Thermosyphon." The Engineer (London), (December 27, 1957) 926-930. 32. Larsen, F. W. and Hartnett, J. P. Effect of Aspect Ratio and Tube Orientation on Free Convection Heat Transfer to Water and Mercury in Enclosed Circular Tubes. ASME Paper 60-SA-21, 1960. 33. Hallman, T. M. Combined Forced and Free Laminar Convection Heat Transfer in Vertical Tubes with Uniform Internal Heat Transfer. ASME Paper 55-A-73, 1955. 34. Ostrach, S. Laminar Natural Convection Flow and Heat Transfer of Fluids With and Without Heat Source in Channels With Constant Wall Temperature. NACA TN-2863, December, 1952. 35. Hammitt, F. G. and Chu, P. T. Transient Internal Convection Heating and Cooling of Closed, Vertical Cylindrical Vessels. Accepted for publication by Am. Nuclear Society, available as Industry Program Report IP-452, College of Engineering, The University of Michigan. 36. Hammitt, F. G. and Brower, E. M. ASME Paper No. 58-A-212, 1958. 37. Chandrasekhar, S. and Elbert, D. D. "The Stability of a Layer of Fluid Heated Below and Subject to Coriolis Forces." Proc. Roy. Soc. (London), 198, A231, 1955. 38. Flutz, D. and Nakawaga, Y. "Experiments on Over-Stable Thermal Convection in Mercury." Proc. Roy. Soco (London), 211, A231, 1955. General References 39. Perry, Chemical Engineers' Handbook. McGraw-Hill, 1950. 40. Ince, Ordinary Differential Equations. Dover Press. 41. Rainville, Intermediate Differential Equation. John Wiley and Sons. 42. Katz, D. L., et al., Liquid Metal Handbook, AEC. 1950. 43. Milne, Numerical Solutions of Differential Equations. McGraw-Hill, 1953. 44. Kaplan, Advanced Calculus. Wesley Addison Co.

-179 -45o Schlichting, Boundary Layer Theoryo 2nd Edition, Pregramon Press, 1960o 46. McAdams, Heat Transmissiono 3rd Edition, McGraw-Hill, 1954o

UNIVERSITY OF MICHIGAN 3 9011111115 02829 5171 3 9015 02829 5171