THE UN I VE RS IT Y OF MI CHI GAN COLLEGE OF ENGINEERING Department of Mechanical Engineering Heat Transfer and Thermodynamics Laboratory Technical Report No. 1 TRANSIENT, LAMINAR, FREE-CONVECTION HEAT AND MASS TRANSFER IN CLOSED, PARTIALLY FILLED, LIQUID CONTAINERS us-sein Z.;:4 7 r~''~~~~~~'' -'.:''. A..,a:,.-',,.;... ORA Project 04268 under contract with: NATIONAL AERONAUTICS AND SPACE ADMINISTRATION GEORGE Co MARSHALL SPACE FLIGHT CENTER CONTRACT NO. NAS8-825 HUNTSVILLE, ALABAMA administered througho OFFICE OF RESEARCH ADMINISTRATION AUN ARBOR January 1964

LvIkvtc

TABLE OF CONTENTS Page NOMENCLATURE v ABSTRACT vii I. INTRODUCTION II. REVIEW OF THE LITERATURE 1 Ao Analytical Studies 3 B. Experimental Studies 8 III THEORETICAL ANALYSIS 13 A. Formulation 13 1. Statement and physics of the problem 13 2. General model 14 3. Simplified model 18 Bo The Solution. 22 C. Finite-Difference Approximations 23 D. Stability and Convergence 27 IV. RESULTS 34 Ao The Effect of Gride Size 37 B. Comparison with Experiment 38 V. CONCLUSIONS 4o ACKNOWLEDGMENT 41 REFERENCES 42 APPENDIX. COMPUTER PROGRAM 45 iii

NOMENCLATURE a half the width of the container b the initial height of the liquid Cp constant pressure specific heat, Btu/lbm~F Gr Grashof number = g(Ts-TO)a /v2 g the acceleration of gravity, ft/sec2 hfg latent heat of evaporation or condensation Btu/lbm k thermal conductivity, Btu/hr ft ~F p pressure Pr Prandtl. number = 7/OC Ra Rayleigh number (GrPr) (q/A)w heat flux at the walls of the tank per unit area, Btu/hr/ft2 T temp R t time, sec u x-component of the velocity, ft/sec v y-component of the velocity, ft/sec U dimensionless x-component of the velocity V dimensionless y-component of the velocity x axial distance, ft X dimensionless x y transverse, or normal distance measured from center line, ft Y dimensionless y v

a2 a2* a2+2 W -- + b 2 ax2 y2 thermal diffusivity, ft2/sec 3 coefficient of thermal expansion p density, lbm/ft3 [ viscosity, lbm/ft/sec v kinematic viscosity, ft2/sec, dissipation function for two-dimensional incompressible flow is given by 4= 4 + + a T dimensionless time @ dimensionless temperature stream function h an eigenvalue Subscripts c cold wall g vapor h hot wall s saturation or liquid surface i,j denotes position in the space grid o denotes initial conditions w wall Superscript n denotes the time level vi

ABSTRACT The two-dimensional, laminar, transient, natural convection heat and mass transfer in a closed rectangular container with a free surface is studiedo The x-momentum, the y-momentum, and the continuity equations are coupled to obtain the vorticity transport equation, which is integrated numerically using the finite-difference approximation0 The problem of stability, which is associated with the difference equations, is studiedo The convergence of the solution of the difference equations so that of the differential equation is examined. Results of the calculations for different boundary conditions are given0 Comparison of the theoretical results with two experimental measurements from the literature indicate qualitative agreement vii

I. INTRODUCTION This report presents the initial results of a program of research dealing with transient, free convection inside closed containers. The problem reported is the first phase of a program directed at an understanding and prediction of fluid behavior in closed containers subject to various disturbances, including those of ambient heat flux, change in wall temperature, pressurization, change in gravitational field, and liquid discharge~ The objectives of the general program are the prediction of transient velocity profiles, temperature stratification, and pressure histories in such containers. The present report treats laminar transient free convection in a twodimensional container having a liquid-vapor interface. Transient processes are introduced by a sudden increase in the temperature of the container wall and by a sudden application of heat flux to the external surfaces of the container. For a system at constant pressure and at normal gravity the resulting transient flow patterns and temperature stratification are computed. The problem is formulated from the complete Navier-Stokes equations coupled with those from the First Law of Thermodynamics and the Conservation of Mass. Boundary-layer approximations are not made since the geometry and boundary conditions invalidate complete boundary-layer flow and require additional momentum considerations not usually found in boundary-layer calculations. Reduction of the governing equations is done by numerical procedures using 1

an IBM 7090 digital computer. Special attention is given to the problem of convergence and stability of the numerical solution. Combined with the experimental program presently under way, these results will aid in the establishment of criteria for the onset of transition flow in the container, the start of surface boiling, and will be useful in formulating the solution for the case of turbulent convection.

II. REVIEW OF TBE LITERATURE, In the past, considerable effort has been given to the study of natural convection heat and mass transfer. Many problems have been solved under different conditions of geometry and arrangements. These studies have been of an analytical as well as an experimental nature. An extensive survey of literature showed that the problem of natural convection in closed-end tubes with a free surface has not been considered analytically, however, some problems with related geometries have been considered. The extrapolation of the results of these studies to other geometry can give misleading results owing to the complicated nature of transient free-convection phenomena. A. ANALYTICAL STUDIES The case of a vertical element immersed in an infinite fluid initially at rest has received the most attention of many investigators. The time25 steady laminar flow equations were first solved by Pohlhausen 5 for air. The experimental results of Schmidt and Beckman2 are in good agreement with Pohlhausen's solution. Later, Ostrach20 solved the same problem using numerical methods with high-speed digital computer for different values of Prandtl number ranging from 0.01: to 1000. The transient free convection from vertical flat plates, with and without appreciable thermal capacity and. variable fluid properties, has been studied by different investigators for different boundary conditions~ Leitzke considered the steady-state natural convection between two

parallel, infinite flat plates oriented in the direction of body force in which one plate is heated and the other is cooled uniformly. The measured temperature distribution across the fluid is in good agreement with the theory. A generalization of the same problem was carried on by Ostrach,' in which the plates are maintained at constant temperatures not necessarily equal, and the effect of heat sources and frictional heating was included. As anticipated, the effect of heat sources and viscous heating increases the temperatures and the velocities between the plates. The transient free convection in a duct formed by two infinite parallel plates with arbitrary time variations in the wall temperature and the heat generation was studied by Zeiberg and Mueller~57 The two-dimensional steady-state convection in. a long rectangle, of which the two long sides are vertical-boundaries held at different temperature and the two horizontal boundaries either insulated or have linear temperature distribution, was considered by Batchelorl1 He did not solve for the velocity or temperature distribution, but he considered the determination of the rate of heat transfer between the two vertical boundaries and the type of different flow regimes that occur for a given value of Rayleigh s number and aspect ratio oFor Rayleigh s numbers less than 103, Batchelor uses a power series expansion in terms of Rayleigh's number Ra for the dimensionless temperature 0 and the stream function. On substitution of the power series in the governing differential equations and equating coefficients of the like powers of Ra, the problem is reduced to the solution of a series of linear partial differential equations0 The Nusselt number, defined as

NU= (q/A)w/K(Th-Tc) is estimated to be of the order -8 2 NU = /d + 10 Ra where d is the distance between the plates and ~ the height of the duct. For the case of ~/d + oo, he argues that for the regions not near the ends the temperature and the stream function take their asymptotic values, which are given by the solution of two ilnfinite parallel plates, one heated and the other cooled. For infinite values of Ra, he postulates that there is an isothermal core having constant vorticityo He found that the governing equations for the general. case could not be represented by a polynomial of small degree nor could they be handled- by the Oseen type of linearization. 26 Poots solved the same problem handled by Batcheloro He obtained a numerical solution based on the-use- of-orthogonal polynomials for the solution of the governing differential equations. Following Batchelor, the stream. function and the nondimensional- temperature were assumed to be represented by the complete double series of orthogonal functions 00 00 0= 7-~ m Anm sin nitx sin, =my and 00 00 = I i BnmXn(x)Ym (Y)

where the A and B are constants which were evaluated numerically. The governing differential equations were-reduced to two coupled algebraic equations to be solved simultaneously. The functions Xn(x) were chosen to satisfy the fourth order Sturm-Liouville system and the orthogonality property, The method of solution is tedious, and the calculations are practically impossible for RayleighCs numbers greater than 104 and aspect ratios greater than 4. The results for Rayleigh's number 104 and aspect ratio unity are shown in Figo 14. There is an isothermal core having uniform temperature; there is also a constant vorticity in the core. In the region between the core and the cavity walls, the temperature and the stream function oscillate once near the core and then tend smoothly to their appropriate values at the cavity wall, The boundary layer is continuous between the wall and the core. Lighthill examined natural convection flows generated by large centrifugal forces in a tube closed at one end and open at the other end to an infinite reservoir, where the tube walls are maintained at a constant temperature. Such a situation exists in cooling gas-turbine blades. He predicted that one of the following three- regimes may exist, depending upon the product of Grashof number and the radius-to-length ratio of the tube. The assumed flow regimes are1 Similarity flowo For small values of this product, i.e., for large values of length-to-radius ratio for a given Grashof number, the boundary layer fills the tube. The velocity and temperature profiles are fully developed. He predicted that for this type of flow, the velocity and temperature distribution are similar at each section of the tube, only their scale is 6

increasing as the orifice is approachedo Assuming that the velocity and temperature vary linearly along the tube, he concluded that there is an aspect ratio for which the temperature changes from its value at the orifice to the value at the bottom, Extending the tube beyond the length determined by the above ratio, the additional length is filled with fluid at rest at the wall's temperature. 2. Boundary-layer type of flow. For high values of the product of Grashof number and radius-to-length ratio, i.e., for short tubes, the flow is of boundary-layer type. In the extreme case, when the boundary layer fills a negligible portion of the tube area, the flow approximates the freeconvection flow up a flat plateo 30 Non-similarity regime. This is the type of flow predicted to exist for values of length-to-radius ratio which lie between the values corresponding to the first and the third case. The boundary layer fills a large portion of the tube section. He used the Squire technique to solve the first and the third caseso 9 Hammitt considered the case of a closed vertical cylinder with internal heat generation. He used the Lighthill technique, modified to account for the heat sources. The agreement between the calculated and measured values of Nusselt's number is not good. This is probably due to some of the inevitable assumptions which are made: (1) small inertia forces compared to buoyancy and shear, (2) radial extend of the temperature and velocity boundary layers is the same, (3) the boundary-layer approximations apply. The first assumption is valid for large Prandtl numbers; the second is valid for 7

Prandtl numbers near unity~ The disadvantage of this method of solution is that it is in capable of detailed examination of the end conditions. Following Lighthill, Ostrach and Thornton considered a geometrically similar case with a linear wall temperature. In Ostrach's paper, as well as in Lighthill's paper, attention was given to the stagnation of natural convection flows at the closed end. The same problem considered by Lighthill was solved by Levy,17 using integral methods. He assume that the upward flow consists of a layer of thickness 6 near the wall, the remainder of the tube being filled with cold fluid flowing downward. He assumes three regimes of flow similar to those postulated by Lighthillo If the tube length. is less than or equal to a length Q, the stagnation region does not exist and the upflow convective layer increases with x0 For axial distance x > Ill 6 reaches a constant value d. and such a flow occurs for Q1 < x < 220 For x > ~2, there is a stagnation region at the closed section of the tubeo 28 Romonov, also using integral techniques, solved the same problem considered by Lighthill and Levyo His calculations agree with those of Lighthill for infinite Prandtl number, but they differ considerably for Prandtl numbers near unity0 The measured and the calculated temperatures are in a good agreement for different wall temperatures0 Bo EXPERIMENTAL STUDIES A large number of experimental studies have involved flat plates immersed in an infinite fluid at rest and either heated or cooled. In general. there has been good agreement between theory and experiment~ Considerable exper8

imental work has been done in the field of natural convection in tubes and enclosures. This work has concerned specialized applications and particular configurations. Most of the experiments were in connection with cooling gasturbine blades and nuclear reactor applications. Although most, if not all, of these experiments are not applicable to this study, they will be helpful in indicating the general trend and the type of experimental equipment required. Probably the most comprehensive experimental studies of natural con19 vection in thermosyphons are those conducted by Martin in an attempt to check the theoretical work of Lighthill. His results agree qualitatively, although the measured heat-transfer coefficients are nearly twice as large as those predicted by Lighthill. The three regimes predicted by Lighthill were identified from measurements of heat transfer rates. The heat transfer rate was greatest for large values of theproduct of Grashof number and the radius-to-length ratio. The rate was highest at the bottom of the tube, which indicates that there is boundary-layer type of flowo At small values of the product, theheat transfer varied linearly from the orifice to zero at the bottom of the tube, from which Martin concluded that there is a similarity regime. A region of instabilities, characterized by nonsinusoidal oscillatory flow, occurred between the above two steady regimes. Siegel and Norris3 shed some light on the oscillatory flow mentioned by Martin by exploring the air flow patterns in the space between two heated wide plates, closed at the bottom, open at the top, and insulated at the sides. For spacing of 0~28 of the plate height, the flow pattern was symmetric, with 9

upward flowing boundary layers near each plate surface and downflow in betweeno When the spacing was reduced to 0o21 of the height, the flow pattern became asymmetric, with half the cross section occupied by upward (near one plate) and the half near the other plate occupied by downward flow. For smaller spacings, the asymmetric pattern presisted with periodic non-sinusoidal reversal in flow direction and temperature fluctuations, 4 Curren and Zalabak conducted an experimental investigation to determine the effect of length-to-diameter ratio of closed end coolant passages on natural-convection water cooling of gas-turbines~ They reported no significant difference in the heat transfer for the different length-to-diameter ratios investigated ranging from 5~1 to 2505~01o For the largest length-todiameter ratio, 25o5: l, the boundary layer fills 87% of the tube cross section. The visual studies of Sparrow and Kauffman36 of free convection of water in a narrow vertical enclosure, cooled at the top through a copper surface and open at the bottom to a heated reservoir, revealed that the flow pattern is not steady. No region of the enclosure is permanently a region of upflow or of downflowo The size of the various upflow and downflow regions varied along the length of the enclosure at a given time~ The number and size of upflow and downf low regions also varied with time. However, end effects were observed, and a continuous downflow took place in a 3/4inch band adjacent to both, wallso Generally, the dominating character of the flow was instabilityo Hartnett, et ali,o'5 studied the free-convection heat transfer for 10

the geometry postulated by Lighthill, but with a constant heat flux at the tube wall and water and mercury as working fluids. The effect of inclining the tube was also investigated. Temperature oscillations of the same nature as those reported by Martin and Siegel and Norris were observed. Contrary to the results reported by Curren and Zalabak, the heat transfer was considerably influenced by length-to-radius ratio. A decrease in length-to-radius ratio from 22.5 to 15 results in approximately 100 percent increase in the Nusselts numbers. Skipper, et al., studies the natural-convection flow pattern in viscous oil in rectangular tanks heated at the center by vertical coil heater. The flow pattern consisted of a narrow chimney of hot oil rising vertically around and above the heater surface, and a horizontal layer of hot oil at the free surface separated by a sharp vertical gradient from the remaining cold oil below. The hot oil layer had a small vertical temperature gradient, with maximum temperature at the top. The hot oil layer at the surface became increasingly thick with continued heating. The hot oil was found to flow downward at the walls of the tank, while there were suggestions of circulating currents at the side of the rising chimney. The flow pattern shown suggests that a vortex was formed at the free surface near the center line where the hot rising chimney is bifurcated and spread horizontally along the 6 surface. Eichhorn observed similar vortices. These vorticies were formed at the free surface of water near the walls of a cylindrical tube 2 inches diameter and 5-inches long, uniformly heated at the walls and open at the top. 11

Some of the literature pertaining to the problem at hand will be mentioned later in this study where its application is more important. 12

III. THEORETICAL ANALYSIS Ao FORMULATION 1. Statement and Physics of the Problem A closed container is partially filled with a cryogenic liquid. Initially the liquid and the vapor whici fills the ullage volume are assumed to be in equilibrium at the temperature To, the stauration temperature corresponding to the initial pressure Po. Previous tests indicate that boiling due to heat leaks from the ambient is a surface phenomenon only; vapor bubble formation is confined to very near the surface; the rest of the liquid does not contain vapor bubbles, and a small change in the ullage pressure is sufficient to cause boiling to cease.13 From this initial condition, the wall of the vessel is assumed to undergo a step change in temperature or is assumed to be subjected to a uniform heat flux. Simultaneously, the pressure in the ullage volume is changed to Pa, which may be equal to or greater than PoO The measurements reported in Refs. 3 and 7 indicate that the interface temperature rises very rapidly to Ts, the saturation temperature corresponding to the pressure Ps in the ullage space for these conditions. Buoyant forces, caused by density variations in the liquid, set up natural convection currents. The liquid-vapor system tends to adjust to the new non-equilibrium condition by transferring mass and energy across the interface. The controlling factor is the interface temperature, and each region transfers heat independently of the bulk temperature of the other region. An imbalance of the 13

heat transfer at the liquld-vapor interface is counterbalanced by a phase change, ioe.,o either transient evaporation or condensation occurs at the interface. In this analysis, the problem is formulated in its general form, taking in.to consideration the interaction between the liqulid and the vapor phase. Later the problem is simplified, making it -tractable for analysis without seriously impairi".ng its utilityo 2. General Model Consider a rectangular two-dimensional tank, of width a and height h. The initial. heigh.t of the liquid is b, and the depth of the vapor is co The Vapor origin of'the coordinate system is (o) = b taken at'the middle of the tank. bottom. with x positive in the direction of the I Liquid l d. The g-level is sufficiently gx 1A g x; high, so that the effect of surface tension can be neglected, The location of _.,.........y j tthe liquid-vapor interface at any one time is given'by x X(t)o The two--- 2a......... dimensilonal'transient conditions will be considered. The differential equations governing the velocity and temperature distribution i n the liqulid and vapor regions areo

The momentum equations: The x-momentum equation.: au, u a p at x 6y ax Fa rau rau ae l aT2( al< +h) + 6u 2 - - _ _ + JJ dyk dVJ + (1) a6x L ax 3 6axy a y ay 6 The y-momentumr equation)v+ a6v 6v = p+ The energy equation: t 6ux Y 6at x v y -~a (0~- at i~ (3 where, O dissipation function (aee omenlature) The continauty equation~ ae:p:+ ea(;l +O (4) at ax "a

Initial conditions: T(x,y,o) = Tg(x,y,o) = To (5) u(x,y,o) = v(x,y,o) = ug(x,yo) = vg(x,yo) = 0 (6) Boundary conditions. (a) Velocity boundary conditions. Assuming the no-slip condition to prevail at the tank walls, the following boundary conditions are obtained: u(o,y,t) = u(x, ~ a,t) = O (7) v(o,y,t) = v(x, ~ a,t) = O (8) ug(h,y,t)= -%/(Pg.A) (9) ug(x, + a,t) vg(x, ~ a,t) = vg(h,y,t) = 0 (10) Assuming zero shear stress at the liquid-vapor interface, the interfacial boundary conditions could be stated as: u(X,yt) = ug(X,y,t) = dX (11) From t edt -- (x z,*t)'~Z (xzt) = (12) 6x' = 6x

axis, it is assumed that the y-component of the velocity vector is zero at the axis of the tank and that the x-component of the velocity is an even function of y, ioe.o v(xot) 0 (14) ay, (b) Theiirmal boundary' conditions~ The bottom and the upper surface of the taSn are either insulated, subjected to a uniform heat flux, or kept at a constant temperatureO The walls are either subjected to a uniform heat flux or kept at a constant temperature" Owing to symmetry, the temperature will be an even function. of Yo (I.) aIT (xot) -' (xsot) = O (15) (i)aT, (!~Y~t~) " (2) ao a (oyt) = (hqy t) 0 a6X bo k T o't) k oTg (hyt~) (q/A)w a~:x a = - (5) - at k x at (xS + at) = (q/A)w bo T(x, ~ a,t) = Tg(x, ~ at) = Tw (17) -17

(4) -a. T(X,y,t) = Tg(XIyt) T b - (b,y,t) - Tg (b,y,t) = 0 (18) Any case considered in this work is a combination of the thermal boundary condition, given by Eq. (15), and one of each of the other three boundary conditions given by Eqs0 (16), (17), and (18). Applying the conservation of energy across the interface P hfg(d - - (X t - k (Tg (X In addition to the boundary conditions stated above, we have the equation of state P = f(p,T) for the vapor region. Also, an overall energy and mass balance, considering the system composed of the vapor and liquid regions, is automatically satisfied by any exact solution; it provides a means of checking assumptions or simplifications introduced later in this analysis. As another check on the solution, the net rate of fluid flow across any section. of the a liquid tank is equal to zero, i.e., u dy - O -a 3. Simplified Model Due to the complexity of the problem, only the liquid region is considered, for which a simplified model is adapted. In reality, the amount of evaporation or condensation is small; therefore, the interfacial displacement is neglected, and the interface is assumed to be always at x = b, The pressure of the vapor is considered to be constant; consequently the interface temperature is constanto Constant fluid properties 1, Cp, k and p are assumed. Density varia18

tions are allowed in the x-momentum equation only. The pressure terms and the dissipation function in the energy equation are neglected. Initially the pressure throughout the tank is the hydrostatic pressure. The variations in pressure and density caused by the fluid motion and temperature gradients from initial values are expected to be small. Therefore the following assumptions are introduced~ P =Po +P' and P = Po + P' where Po is the initial hydrostatic pressure, P' the change in pressure from the initial value, po initial density, and p: the change in density. Since the pressure variations are small, the density changes due to pressure are negligible, and density variations are caused mainly by temperature changes. Then pt can be closely approximated by Ip = p o Po (T) (22) where 3 is the coefficient of thermal expansion. Similarly, =P -Po g + (23) ax 6x -ap PP (24) by by 19

Substituting the above expressions in the momentum equations: The x-momentumn = + - a - Po Pg(T-To) - + a a(25) The y-momentum: + u ax+ v ) - +A + (26) The continuity equation~ + V= o. (27) ax ay The energy equation: — +3 U +% +,-t (28) at ax by \x2 The substitutions necessary to reduce the above differential equations to nondimensionalized forms are: u = 2 v V a a ~vb a2 T-To 4 Q t. T x = bX, y = aY (29) 20

S Substituting Eqo (29) in the differential equations, differentiating the x-momentum equation with respect to y and the y-momentum equation with respect to x, combining the two equations, and introducing the stream function, the following equations are obtainedO p2 6aX2 aY22 aYaXa2 aY2 aXaYb2X2 aY2 Pr + Pra a a aa2 - a 2 — a +22 a 2,, (29a) a L Kb2Pb ax2 +y2) Y2?b2 ax2 ay and _a aa 8 +2G (30) aT'aY ax ax aY b2 6X2 Y2, where the stream function' is de:fined by U a= e al condtion The'initial condition:s o(X,Y,o) = (31.) v(x,Y,^O) 0 (32) and the boundary conditions are. 21.

r1(X,L,T) = 6 (XT) x,T) = o (33) Vt(O,YT) = T'(O,YT) = (OYT) = o (34) V(1,YT) = i(1,YoT) 2a(1,Y,T) = 0 (35) 4(X,o0T) = X(XOT) =;~(XO-T) = (36) x(O,Y,T) -(X,0,T) = 0 (37) =x 6;(X1 XT) = Gr*, or - (X,,T) = Ow (38) (1,Y,T) - s, or (,Y,) = (39) where: Gr* = (a/b)Pr gga4(q/A)wl(kv2) (40) Qs = b(Ts-To)Pr ga3/v2 = (a/b).Pr.Gr (41) b From the above system of equations and boundary conditions, clearly U, V, and 0 are functions of x,y,T, a/b, Pr, Gr, and Gr*o B, THE SOLUTION The nonlinear, fourth order system of Eqo (29a) and (30) is not amenable to mathematical treatment using classical methods. Furthermore, Ostromov24 and Batchelorl found that neither the successive approximation method nor the series expansion is suitable for handling such a system of equations for any 22

arbitrary Pr and Gr, Therefore it was decided to employ a numerical method for the solution; the finite-difference approximation was chosen for the work in this report. The advantage of this method over other numerical methods is that the boundary-layer assumption is no longer necessary~ Instead, the full. Navier-Stokes equations can be solved. The results obtained by this method can be used to check the validity of the simplifying assumptions made by others. Co FINITE-DIFFERENCE APPROXIMATIONS The fourth order, nonlinear system of Eqp (29a) and (30) is transformed to a second order, nonlinear system as follows: Let w. ab X a+- (42) b XaX 6y2 Upon substitution in Eq. (29a), the following equation is obtained: a. +V Pr vPr+ (43) a +u a V2 a2 (44) - + X+ va -Y b2X2 bY2 By definition, w is twice the negative of the vorticity, and Eq. (43) is known, as the vorticity transport equations For convenience U and V are substituted for at/bY and - b4/bX. By this transformation, the system of Eq. (30) and (31) is reduced to three second-order, nonlinear equations which are 23

easier to handle by difference methods compared to the two nonlinear, fourth order Eqso (29a) and (30)o The differential Eqs. (42) (43), and (44) are approximated by a system of difference equations. Each first order term in Eqs. (43) and (44) can be expressed in forward, backward, or central differences. The choice of the type of difference approximation is determined by whether implicit or explicit difference method is desired and by the stability requirements, to be discussed later. From the beginning of this analysis, explicit methods were felt to be more suitable for handling Eqso (43) and (44), because the principal interest is in transient phenomena for small physical time0 This implies that, for better accuracy, small time steps should be used from the beginning of computation0 Desinberre5 showed that explicit methods are superior to implicit methods as far as accuracy is concerned; he also found that the accuracy requirement restricts the size of the time increment, in the case of implicit methods, to values smaller than required by stability considerations0 When explicit methods are used, the nonlinear terms U(Cw/6X), V(~w/6Y), U(6Q/6X), and v(o/6.Y) have to be approximated by forward difference when the sign of the velocity is negative, and by backward difference when the velocity is positive.35 Central differences are used if implicit schemes are desired. Since the velocity components change sign in the domain ofinterest, four different sets of finite-difference equations are required, and the appropriate set of equations at each location or nodal point is determined accordingly0 The four sets of difference equations corresponding to differential Eqs0 (43) and (44) are: 24

(i) U >O V>, Wji i Ji -1Wi. Wi -i J -Wi 1 WTi,.J + Ui j W+i j..... AT?' AX 1 AY Pr i-'i j1 2 P Wi+l,-2wi J+Wi lj Wi +.-2wi'j +wi 2j 12 j_ __I_ _ _ __w,,__ _ _ _ _ _ Pr + PrI4 -Pr w j + 31 AY b2 (AX)2 ()245 Q' 8 j 33 1i + j 1 + U.0 1X + j, 1 (AX), (A Y) 2 ~ij-~iJ + Ui, ~-2 @ i, j-i-,j- (46a) (ii) U O V < O Wa i ji - 2i W. 2~i i+9 +z -'gUr. 1 - 3 Wji j1li j Pr i + Pra W -2i.+W +i-l,+j Wi Jw i -2w +Wi j_A Ty j + Ui, e i- j +: i +, 0 i j +L i, A i- AX 1 AY a2 ~ql 28 e o e 1 j+i, j +1 i 2 j i n j 8 a2.. e, +-.,.~ 1+, r20... AT2x()2 (AY)2 (46b) am ~i~~l~j-25ij0ll i jl2i.Q~j

Ii -d<'-d P~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~v ~~~~~~~~~~~~~~~~~~~~~~~~' Hl ~~ A ~'~ IA sl P HJ- 0 -F 0.a d 1,r FJ —~'. IA LH~~~~~~ [~~~I IA c-i L~~~~~~~~ ~u F- H-I, ~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~'2 %1~~~~~~ H''+ - + I' R... Jo HH ro m PO > ~ -. mR N ) HC- J. C. i ~+ II- + tH CJ +'H hi + + +: - (D ( D H" >1C-J. FjH 4 _. C-_. C —J. C-J. H C-a C-J_ R.) I I I > ~ —~ ~ Ho "-~ ~ H. I'o ~. H- + I~ ~. ~. 4 - H+ ~~~~~~~~~~~~~~~~~~C-J. +' M-p I::>,-J c.. I+:C-J. ~ +~~~~~~4 + + + H'- H' H + + > +:\ I:> - C-a. H' Hi H-1. C-iI H' HC +>< c. > <I'-. H oh~~~~~~~~~~~~~~~~~~~~r tc__,. _. C-. CI>. C -J. C-J. C-. -. -~~~1- e_, F _. H' I H. ~~~j ~ - Ci; -F C-iC-J. ~ ~ ~ ~ ~ ~ ~ ~~~+C-J. O C ll ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~~'~. I-'i_.oM ~~~~~~~~~~~n~~~~~~~~I~~~~~~~~~ I i>- O H H' C~~~~~~. I~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~:-J'- ~ "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~J,.

The prime refers to the value of the variables at the present time step, and the unprimed variables correspond to those at the previous time step~ The two-dimensional space domain of interest is bounded by 0 < X _ 1, -1 < Y < 1o Only the region to the right of the x-axis is considered because of geometrical symmetry with respect to this axis.. The integers i and j denote the x and y position respectively, with. i = 2, j = 2 as the origin, X = (i-2), Ax and Y (j-2)o Do STABILITY AND CONVERGENCE The method used in this section follows that of Refs. 27 and 12, The theoretical background for the method may be found in Refo 27, which presents an excellent discussion, of the subjects Stability is a necessary condition for the solution of the difference equations to converge to the solution of the differential equations as the size of the increments Ax, Ay, and At tend to zero. In other words, it is important to identify the behavior of the approximation error as the size of the increments goes to zero and to establish the restrictions to be imposed on the difference equations in order that the error will tend to zero as the increments Ax, Ay, and At go to zero9 This requires that there be no unlimited amplification of the error as the computing cycles become infinite in the limit9 The Von Neumann method of stability analysis is used in. this problem9 This method is strictly valid for constant-coefficient difference equations while the equations involved have variable coefficients eT and V; however, 27

very good approximation for the stability criterion is obtained by treating the variable coefficients as constants throughout the analysis and then letting the variable coefficients take on their most adverse values in determining the restriction on the time increment size. The fundamentals of the application of this method of stability analysis are best illustrated by the problem in hand~ The method will be applied for the case U > O, V > 00 The stability criterion holds irrespective of the sign of U or V. The absolute value of U and V must be used in the equations or inequalities describing the stability criterion, and the appropriate difference equations must be used according to the sign of U and V, as described earlier, The solution of the difference equations can be written as a Fourier series, the form of which is 27 w,(n)j (n) Z i(K1X+K2Y) (47) K_ -K2;,(n). X XP (n) i(K1X+K2Y) K1 K2 where K1 and K2 are integers, n is a superscript denoting the nth time period, and. and p are functions of K1 and K20 Substituting the system of Eqso (47) and (48) into Eqso (45a) and (46a), the following equations are obtained after some algebraic manipulations~ 28

Ka j K(n+lK S(n)eKAX -iK2AY a iKAX + a eiK2AY, ),L? l+a2 e + a3 e + a4 e' + as e K1 K + a6e (n+l) i(K1X+K2Y) = -IKAe X -iK2AY iKIAX iK 2A Sz. - ~ (+l2 e + c e K1 K2 e (KIX+K2Y) 0 From the above equations, it is concluded that the difference equations are satisfied if (n+l) (na) -iK AX -iK2Y KA Y iK2A2AY a a2 e + a3 e + a4 e + a e + a (49) and ( n + l ) - i. i),iX - iK pAY i.W AX 2 = l2 e + c e - + c4 e + 5 e where - a1 = l_ 6aA2m X a 1 4 a Pr Pr + IUJ + AT 2a2 = 2 (Pr + ) a3 = (Y)2 +i AY 29

a Pr a4 - b2 () ATX) AT a5 Pr (y) c2 = 2__ _i I _V IQ (b ax2 + / a. Jv5ca3 AT No definition has been given. to a6 since it has no effect on this analysis. The system. of Eqso (49) and (50) are of -the form~ (n+-L) __ (n) ~ __ (n) _ _ (~"") (') i n Kg~ ) a"i2 () (K aK2) (51) =" al "' (tKaK2) 822 2 (K1,K2) (52) In the matrix notation, the above equalities can be written as (n+l1) 2 a (n) ( a0 a (a 30

The quantity between the first parentheses on the right-hand side of Eqo (53) is called the amplification matrix. The Von Neumann condition necessary for stability is: I1 maxl < 1, where k max is the largest eigenvalue of the amplification matrix. The eigenvalues are given by rall-. a 12 I 0 a21 a22-x Substituting the values of all, a2p oo etco in the above determinant and solving for X. we get 1 al + a2 e + a e + a4 e + aY e (54) 2 =:= C1 + C2 e-iK1AX + C3 -i e + c eiKAX + eiK (55) The coefficients al, a2,aoCl9c20o.etcO all are positive except a, and cl. which may be positive or negative. The largest absolute values of Xn and X2 occur when all the terms in Eqso (54) and (55) are real, ioe.o when K1AX K7AX = K2AY = 2e thenn, kX max ai a2 + a3 + a4 + a5 (56) ~2 ma = C1 + C2 c3 + c4 + c5 (57) Substituting the values of al, a2,oo0, CioooCs5 in X max x~ = xa= 1 (58) max max( 31

Therefore, we can conclude that % max will not exceed unity and will not impose any stability restrictions. If there be any restrictions, they are to prevent the minimum value of k from becoming less than -1o The minimum of the eigenvalues occurs when K1AX = K2AY = i and is given by = a - a, a2 - a3 a4 - a5 m.in h2 min = cL - C2 c3 C- c C or min 1 - 2 AQ Pri +2 + L + I i (59) min (AY)T A A(5 X 1 - 2 AT + (A* Ui' I + I ~ Ij (60) h2 min I' i. -. 2 2 (AY),a Ai AY Therefore, for IlX < 1 the following inequalities should be satisfied: AT(a2 + + IiJ+ <23) 1 (61a) An(2 (~) ~AX AY — Equations (6lamb) are requisite for stability~ For values of Prandtl number less than unity, inequality Eqo (61a) is more restrictive and, therefore, should be used~ For higher values of Prandtl number, inequality Eqo (61lb) must be usedo 32

The numerical solution is carried as follows: 1. The temperature distribution is first calculated using Eqo (45); 20 The advanced values of temperature Q' are used in Eqo (46) to calculate wt; 3~ The stream function is calculated at each time step, using Eq. (42). The solution of Eqg (42) is done numerically, using successive row relaxation followed by successive column relaxation36 33

IVo RESULTS Calculations have been carried out for the case of a container with an insulated bottom whose walls are subjected to a uniform heat flux, and the liquid surface is maintained either at the initial temperature T,, at a temperature, Ts, higher than To, or is adiabatic. The last case is expected to approximate the transient convection in the liquid when exposed to a nonpressurized gas. Different levels of heat flux were considered, including 10 and 1000 Btu/hr/ft2o For these calculations, the fluid properties chosen were those of liquid nitrogen initially at atmospheric pressure, the initial saturation temperature was 140~Ro The fluid properties were evaluated at a temperature equal to the average of the initial temperature and the liquid surface temperature. The values of the liquid properties were taken from Ref 38 and are summarized in Table Io The height of the liquid b, is 1 ft, and the width of the conTABLE I FLUID PROPERTIES Thermal diffusivity U, ft2/sec 8~62 x 10-7 Thermal conductivity K, Btu/hr/ft2o~R 0-0775 Kinematic viscosity v, ft2/sec 1. 68 x 10-6 Coefficient of thermal expansion, B, ~R- 1 l533 x 10-3 Prandtl number - lo 91

tainer is 1/2 fto The flow pattern for the case of a constant wall heat flux of 10 Btu/hr/ft2 and liquid surface maintained at initial temperature To is shown in Figso 1, 2 and 3, using a 21 x 21 grid corresponding to AX = AY = 0~05~ These results show the streamline pattern at different time levels —50 sec, 2 min, and 3~6 min. respectively. An interesting streamline pattern is observed at the free surfaceo For a short time after the beginning of heating, the boundary layer rising along the container walls turns smoothly changing its direction from upward to downward flow (Fig. l)o The downwardmoving particles near the rising boundary layer reverse direction and join the upward flow, thus giving rise to the vortex near the free surface0 The fluid away from the edge of the boundary layer flows nearly to the bottom of the container, where it joins the upward-moving fluid0 For greater times following the introduction of the transient, the streamlines near the free surface show the presence of fluid oscillations (Figs~ 2 and 3)~ These oscillations first form near the wall, their amplitude grows with time, and they move towards the centerline of the container~ The calculations show that this oscillatory phenomenon is repeated with time in the sequence described (Figs0 2,349, and 5)0 The effect of increasing the level of heat flux on the flow pattern is clearly shown in Fig0 4, which shows the flow pattern obtained after heating for 51 sec with a wall flux of 1000Btu/hr/ft2 and the surface temperature maintained at the initial temperature Too The comparison between Figso 1 and 4, which correspond essentially to the same time, shows that the oscillatory 35

streamline pattern develops earlier for high heating rates than the low heating rates. Except for the latter effect, the flow pattern at the higher heating rates has the same characteristics as that for the low heating rate. The magnitude of the velocities, of course, is higher for the higher heat flux The streamlines for the case of a constant wall heat flux of 1000 Btu/hr/ft2 and adiabatic interface are shown in Fig. 7. Figures 8 and 9 show the streamlines and the isothermals, respectively, obtained at 1000 Btu/hr/ft2 wall heat flux and the interface temperature maintained at 1600R, which corresponds to the saturation temperature for a pressure of 45 psiao The axial temperature gradient in the boundary layer is negligible for about 70 percent of the container height at a time of 48 sec. At the upper portion of the container, near the free surface, the axial temperature gradients are considerably greater. On the other hand, the transversal temperature gradient is greater at the lower portion of the-tank and becomes smaller near the free surfaceo The temperature distribution exhibits the same character in all the cases analyzed. These phenomena can be explained as follows~ for small times, the fluid near the container walls flows upward in a thin boundary layer. In its upward movement, the hot fluid entrains some of the cold fluid at the edge of the boundary layer. This heated boundary layer is discharged at and just below the free surface, where its transverse velocity is highest. To satisfy continuity, the heated fluid which is discharged at the free surface causes the colder fluid to move downward, thus producing a series of horizontal isotherms.o With time these 36

isotherms penetrate further below the free surfaceo At the lower portion of the container, the transversal temperature gradient is very high near the wall and negligible in the remainder of the container. It is smaller near the free surface, where the boundary-layer flow is discharged. The thermal boundary layer fills the entire cross section in this region~ A0 THE EFFECT OF GRID SIZE Calculations were made for the case of constant wall heat flux of 1b Btu/h~r/ft2 with the fluid surface at its initial temperature, using 21 x 21, 16 x 16 and 11 x 11 grids. The streamline patterns are shown in Figs~ 1-4, 6. The oscillatory pattern is obtained in each case, but it is on a smaller scale for the 11. x 11 grid. The axial velocity distribution obtained at a given time is plotted for various axial positions in Fig~ lO1 The values obtained using 16 x 16 and 21 x 21 are in close agreement. Figure 10 shows that the difference between the values obtained using different grids is greater near the wall, and at the center line; otherwise they are in substantial agreement Figure 11 shows the dimensionless wall temperature at a location x = ~6 plotted against dimensionless timeo The deviation is greatest for small times0 The difference decreases with time and is practically negligible for dimensionless time of 0o003, which corresponds to about 3~6 mino Figure 12 presents the velocity at different locations as a function of the grid size0 Examination of this figure reveals that the velocities change monotonously with the grid size0 Their values are expected to converge 37

to the true solution for a number of divisions of the order of 30 x 30 This is seen from the extrapolation shown as dotted lines in the figures. The rate of convergence is slower near the boundary and at the center than at intermediate positions. The computations were done on the IBM 7090 digital computer at the Computing Center of The University of Michigano The machine time required to complete the calculations and to print the results for UVT,w, and'V every 20 steps is 1l7 sec per time step for the 11 x 11 grid~ The corresponding times for the 1.6 x 16 and the 21 x 21 grids are 8 and 11 sec, respectively, for printing the results every five stepso Bo COMPARISON'WITH EXPERIMENT An investigation of the literature on natural, convection showed that few present results are applicable to those reported here. These include 6 7 the work of Eichhorn and SO Ko Fenster, et al, o there is other literature concerning flow phenomena in free convection, but, except those cited, none included systems with. a free surfaceO Eichhorn conducted visual, studies of the natural convection laminar flow of water using an electrically heated cylinder 2-inches diameter and 5-inches longo His results are shown in Figo 130 The magnitude of the heat flux was not giveno From the discussion it is concluded that the results represent the unsteady state. Figures 13a and 13b show the flow pattern. observed at high:heating rate; Fig~ 13c shows that obtained at low heating rates~ At low heating rates; the streamlines assume a damped-wave shape~ At high heating

rates, annular vortices repeatedly form near the free surface, roll up until a certain size is reached, whereupon they move away from the cylinder, and another vortex vegins to formo The comparison of Figs. 1-8 with the results of Fig. 13 shows that the shape of streamlines observed agrees with that obtained from the theoretical solutions presented hereo Fenster, et al,, studied experimentally the transient phenomena associated with the pressurization of liquid nitrogen initially boiling at constant heat fluxo Although the initial temperature distribution agrees with that assumed in this analysis,, the initial velocity distribution does not. At any instant after pressurization, there was no difference between the temperature at the tank centerline and midway between the centerline and the wall at an axial location below 0O6 of the liquid height~ These results indicate that the isotherms are horizontal in the core of the tank, which agrees with the calculated isotherms shown. in Figo 9o Temperature oscillations of the type shown in Fig. 9 were also obtained 26 by Poots for the steady state solution of the temperature distribution in a two-dimensional closed cavityo In this case, the walls are kept at a constant temperature, one wall hotter than the other, and the upper and lower surfaces assume a linear temperature distributionO These results are shown in Fig. 14o 39

V. CONCLUSIONS The equations describing the transient, two-dimensional, laminar, natural convection in a rectangular closed container having a vapor-liquid interface have been solved using explicit finite-difference approximations. The velocity and temperature distributions indicate that, for small time periods, the flow is of a boundary-layer type, except near the bottom and the liquid-vapor interface. No indications of numerical instability were encountered. The size of the time increment was restricted to small values by stability considerations based on the method of Von Neumanno However, a small time increment is desirable in order to obtain accurate transient results at small time. On the other hand, if steady state results are desired, the number of computations using the explicit method will be large for fine grids and for systems having a high heat flux.o The machine time may be of the order of 2 hours in this 12 case, and the use of implicit methods would be superior. However, the application of implicit-difference methods to this problem and to problems with other geometry is being continued. The results obtained using the present method agree qualitatively with related cases reported in the literature. 40

ACKNOWLEDGMENT The authors wish to thank Profso V. S. Arpaci and Wo J. Yang for their valuable suggestions and assistance. The assistance of Professor R.CoFo Bartels, Director of the Computing Center, The University of Michigan, during the initial stages of this work is gratefully acknowledged. 41

REFERENCES 1o Batchelor, Quart. Applied Math. 12, ppo 209-203, 1954o 2~ E. Schmidt and W. Beckmann, Forsch Ing-WesO 1, 391 (1930)o 53~ Jo Ao Clark, Go J. Van Wylen and So Ko Fenster, "Transient Phenomena Associated with the Pressurized Discharge of a Cryogenic Liquid from a Closed Container " Advances in Cryogenic Engineering, 5, ppo 467-480!. 4o A. No Curren and C. Fo ZaLabak, "Effect of Closed End Coolant Passages on Natural Convection Water Cooling of Gas Turbine Bladeso" NACA RM R55J 18-ao 5 Dusinberre, "A Note on the "Implicit" Method for Finite-Difference Heat Transfer Calculations)~ ASME Trans, Jo Heat Transfer, Series C, pp. 9495, Feb. 1961K 6. Ro Eichhorn,'"Flow Visualization and Velocity Measurement in Natural Convection with Tellurium Die Method." ASME Tran.s, J, Heat Transfer, 83_, Series C, ppo 379-381o 7. S. K. Fenster, Go J0 Van Wylen and J. Ao Cl.ark, "Transient Phenomena Assoicated with the Pressurization of Liquid Nitrogen Boiling at Constant Heat Fluxo" Advances in Cryogenic Engineering, 5, pp. 226-234o 8~ Bo Gebhart, "Natural Convection Transients~" ASME Trans., J. Heat Transfer, Series C, 85, ppo 184-1.85o See also Jo Heat Transfer, 85, po 10 and 83, po 61K 9, F. Go Hammitt, "Heat and Mass Transfer in Closed, Vertical, Cylindrical Vessels with Internal Heat Sources for EHomogeneous Reactorso' PhoDo Thesis, Univo of Micho, Dec. 19570 10. J. P. Hartnett, W. E. Welch and Fo W. Larsen, "Free Convection Heat Transfer to Water and Mercuty in an. Enclosed Cylindrical Tube)" Nuclear Engineering and Science Conference. preprint 27, Session XX, Chicago, March 17-21, 1958o 11 J. P..Hartnett and W. E. Welch, "Experimental Studies of Free Convection Heat Transfer in a Vertical Tube with Uniform Heat Fluxo! Trans, ASME, 79 P o 1651, 1957o 42

12 J. o D. Hellums "Finite-Difference Computation of Natural Convection Flowo" PhoDo Thesis, Univo of Mich, Sept. 1960. 13. F. Herzbert, "Effective Density of Boiling Liquid Oxygen." Advances in Cryogenic Engineering, 50 14. Landau and Lifshitz, "Fluid Mechanicso" Adison-Wesly, 1959,.15o F. W. Larsen and J. P. Hartnett, "Effect of Aspect Ratio and Tube Orientation on Free Convection Heat Transfer to Water and-Mercuty in En closed Circular Tubes." ASME, Trans., Jo Heat Transfer, 83, Series C, ppo 87-93. 16. A. F. Leitzke, "Theoretical and Experimental. Investigation of Heat Transfer by Natural Convection Between Parallel Plates." NACA Report 1223, 1955. 17. S. Levy, "'Integral Methods in Natural Convection Flow." ASME paper No. 55g APM-22. 18. M. J. Lighthill, "Theoretical Consideration on Free Convection in Tubes." Quarto J. Mech. and Applo Math., 6, ppo 398-439, 1953. 19. B. W, Martin, "'Free Convection in an Open Thermosyphen. with Special Reference to Turbulent Flowo" Proco Royo Soco, Series A, 230, ppo 502530, 1955. 20. So Ostrach,g "An Analysis of Laminar Free Convection Flow and Heat Transfer About a Flat Plate Parallel to the Direction of the Generating Body Force o NACA Report 11ll9, 1953. 21. So Ostrach, "TLminar Natural Convection Flow and Heat Transfer to Fluids With andWithout Heat Sources With Constant Wall Temperatureo" NACA Tn2863.5 Dec. 1952. 22o So Ostrach,!'Combined Natural and Forced Convection Laminar Flow and Heat Transfer to Fluids With and Without Heat Sources in Channels With Linearly Varying Wall Temperature." iNIACA Tn-3141, April 1954. 235 So Ostrach and P. R. Thronton, "On the Stagnation of Natural Convection Flows in Closed End Tubeso" ASME paper No. 57-SA-2. 24. G. A. Ostromov, "Free Convection Under the Conditions of the Internal Problem,! NACA TMI-1407. 250 Eo Pohlhausen, ZAM 115 (1921)o 43

26. Go Poots,"Heat Transfer by Laminar Free Convection in Enclosed Plane - Gas Layersd" Quart. J. Mech. and Applo Math., XX, Pt. 3, pp. 257-273, 1958. 27. R. D. Richtmyer, "Difference Methods for Initial Value Problem." Interscience Publication, 1957. 28, A. G. Romonov, "Study of Heat Exchange in a Deal End Channel Under Free Convection Conditions," Izvestiya Akademii Nau, USSR. Otdeleni Tekhnisheskikh, No, 6, p'p. 63-76, 1956. 29. R. S. Schechter, "Natural Convection Heat Transfer in. Regions of Maximum Fluid Density." AIChE paper No. 57-HT-25. 30. R. Seigel, "Transient Free Convection from a Vertical Flat Plate." ASME paper No. 57-Sa-8. 31. R. Seigel and R. H. Norris, ":Test of Free Convection in a Partially Enclosed Spaces Between Two Heated Vertical Plates." ASMIE paper No, 56SA-5o 32. RoGoSo Skipper, I.S.C. Holt and 0. A. Saunders, "Natural Convection in Viscous Oilo International Development in Heat Transfer, Pto 5, ppo 1003-1009. 33. E. M. Sparrow and Jo L. Gregg, "The Variable Fluid Property." ASME paper No. 57-A-46. 34. E. M. Sparrow, J. L. Gree, "Similar Solutions for Free Convection from a Nonisothermal Vertical Plate." ASNE paper No. 57-SA-3. 355 E. M. Sparrow and SO J. Kaufmann, "Visual Study of Free Convection in a Narrow Vertical Enclosureo" NACA RM E55 L 14a,9 Febo 1956. 36- J. Todd, "Survey of Numerical Analysis,." McGraw-Hill, Chap. 11, 19620 37. So L. Zeiberg and W. Ko Mueller, "Transient Laminar Combined Free and Forced Convection in a Duct." ASME Trans, Jo. Heat Transfer, 84, Series C, ppo 141-148. 38, "A Compendium of the Properties of Materials at Low Temperature Phase I and II." R. Bo Stewart and VO J0 Johnson. Editors, 44

APPENDIX COMFPUTER PROGRAM The computer program used for the cases of constant wall flux or constant step change in wall temperature and the interfacial temperature specified is given on the following pages. The program is written in MAD language. The symbols UVT,W,To0,sTkA, and B are the same as in the text. The meaning of the principal symbols which are not defined in the program are given belowo Dx Ax Dy = Ay DT = AT M = No of divisions'in the x-direction N 3 0No. of divisions in the y-direction GA Acceleration due to gravity SF Stream function. ST The value of the stream function at the previous time step W The value of W at the previous time step WL= The value of W at the advanced time step TL The value of -temperature at the previous time step NT = Total number of time steps NE Number of interati.ons. TJI = Dimensionless wall temperature Tiw for the case of step change in wall temperature TI = Dimensionless interfacial temperature~ 45

HUSSEIN Z. BARAKAT 5265F 02,0 300 HU-Sfl Z —--------— Al —65F --- Q 3 Q - - - $COMPILE MAD,EXECUTE,DUMP,PRINT OBJECT,PUNCH OBJECT DIMENSION U(2000,'DIM),V(2000,DIM),TL(2000,DIM),T(2000.DIM) 1ST ( 200,DI M),SF( 20 00DIM),WL(2 00 0,DIM) W(2 00,D ),DI(100)........... 2__ ~~_ _L._I _ 1_Q 1 gE Q_,_ E _J_._QO_ L.,___F_ ___________ _.__ VECTOR VALUES DIM=2,1,0.____,,,,,__ ______...1] F -R - -L. M- a * IM — ~ ~ L2 3 L-3- 1 4 A J Cr J_-d- JJ 5[-L Lj_.... 1' N-7,N8 EXECUTE FTRAP. S [ART READ AND PRINT DATA ------------------- __jL_.IE J _C 2___ —-----------— _ -— _ ~ -__ ------------------------- EXECUTE ZE RO.(U (2, 2)...U( M+ 2,N+2),V(2,2 )... V ( M+2,N+2 ),T(,1). 1.. T(M l+2,N+~ 2L TL (1,1)...i L(M + 2,N+2),ST ( 2,2) *..S T ( + 2,N+2), 2SF'(2,2)...SF(M+2,N+2) *W(2,2)..'.W(M+2,N+2),WL(2,2)...WL(M+2, 3N+2),E 1 I (1 I (N+2),Di (1)...D, I ( N+2.) t,F ( 1)... F (N+2 ),EJ (i )... 4 EJ (M+2) ) TS= TS Q S U R F = Q SURF' A=A K=K T I ME=O N —--------------- Q _ _ —-------------- ---------------------------------------------------------- DX= 1./M DY=-1./N _ _ _ _ _ __ __ I 1 =v/ 2 J2=N/2 R=AOVERB DX2= DX*D X DX3=DX2*DX DY3= DY2*DY DXY=DX2/DY2 DYX=DY2/DX2 B I =2.* ( 1+DXYL. BJ=2.*( 1.+DYX) -----------— _ - Ci'=GA*PR-RBEI TA_*QSU RF' (. A P - N E —-------------— /NEw1 THROUGH ZZ,FOR J=3,1,J.G.N+1 EI( 2)= 0 ZZ EI(J)=1./(BJ-EI (J-1)) -Y____YEI________________IHt i_ f__ Y_~E-_13_1_=-_ _li___._+_._ ___ EJ(2)=O._V E. —---- ___-L —- - - --- — Le._-B_ —L_ —-J l —— _________ —— __________ THROUGH TH ETAA,FOR I =2,1, I.G. i+2 TH FTA TI ( T.N+2)=T. J NTC=O a -------— ______-_LLL-L +l_______________________________-_ —-------------------------------------------- THROUGH bE HH,FO.R J=2,i,J.G.N+3 -. —-- -______*-_~ LL Lt+ J__L_ -1______ —-------------- BEHH T(M+2,J)=T1 THROUJGH ALEFF FOR I-=2,1,I G. +l ALEFF T (I,N+2 )=TL( I,N+2 )+2.C 1 D13/ DY+T C 3; 2.'( TL( i,Nti )- LI_-I _+_LLlJD__+K_!_!* T * R*R* ( TL ( 1T L ( I N __ _+;-ii Li i__ 2N+2) )/DX2 T ( I. _ ) ( _ =_ T N3i ji +2 iTH ROUGH AALEF,FOR I=2,1o.G.iv+1 THROUGH AALEF,FOR J=Z,1,J. G.N+1 T( I,J)=T L( i,J)+;DT3;R"T;i R< TL( i j ) -',iTL ( i i +, I)TL( I-i,,J) )/D X 2 46

1+DT3*(L(II,J+1)-2.*TL(I,J)+TL(I,J-1))/DY2 _T ( ____________J___ _L3_JJ_-_ -...... _.... T(1i,1) =T(3 1) AALEF T(I,1) =T (I 3) THROUGH XX,FOR I=1,lI.G.M+l THROUGH XX,FOR J=1,1,J.G.N+2 XX TL( I,J )=T(I,J). —---------------— WhE-NEY-ER LS - -C. L -t TRAjiah LR- TQ A _ _ ___ —----------------------- BACK NT=NT+1 THROUGH SIGMA,FOR I=1,1,I'.G. (M1+2) TiHROUGH SIGMAFOR J= 1,1, J. G. N+2 _____ I MA ____ L ( I, J- L -___I_J _.. —-------------------- WHENEVER NT-.G.N1,DT=DT1 SMAX= _ THROUGH GAMMA,FOR I =,1,I.G. (+1 )'THROUGH GAIMA,.FOR J=3,1,J.G.(N+1) WHENEVER PR.L.1 S'RiR 2./DX2/PR+2/P./ PR/DY2+.ABS.U(IJ)/DX+.AS.V(,J)/DY OTHERWISE S=R*'-R*2_. /D X2+2./ DY2+.AGS. I J / DX+ b ( Y 2 A 1 S L -I,_J_/DY END OF CONDITIONAL WHENEVER SMAX.LE S SMAX=S G A M M A L _-OJI-_-_-__AL-___-_N____ —_I________ SN D TS 5 i AX W H E N E_____RV___K_ -.__1_____HE, Qr, L A-1X tVI X_______ _________1 ___________ TI M TIM E + D T THROUJGH 6FH.FOR.J=1,1.IJ.G.N+? 2 T(M+2,J )=T1..___r __ IILLt 2 lLt- l_________I 1_ —------------------------- ------ TH ROUGH ALEF, FOR I=2,1, I.G. ( M+1) ---------- --------------- 1TL ( I,N+2 ) ) /DY2+TC eDT R'- R'- ( TL ( I +1,iN+2 ) -2.* -TL ( i,N+2 ) +TL ( I-i, 21N+2 ) )/fX 2 THROUGH DELTA,FOR I=2,1, I.G.(i +1 ) ------------—....... Ilt iLl-R IJ fi_ -.- I._,_L- 1 G 1_. _ _ _ __ —---------------------— ____ WHENEVER U ( I,J ). GE.O._______,,,____________ B J -11 1sL;-I — - J-L —------------------------------------------------------------- OTHERWISE RT=TI ( lI +-)-TL (T - I END OF CONDITIONAL -_ELiRYJIALILEQ_W-IE-LE-_V_ G_ __L__ J__-_ E___ C =TL(I, J) - TL (I,J-1) C-T=TL(I,J+1)-TL( 1,J).FNUl) OF CCNDT T ITNAI T (I,J)=T L( I,J),+DT*(T(I+1,J )-2 t*T (I,J)+T(I-1 J ))*R* R/i DX+ U' T* i ______j,_J_____ [_i_ ",._L___+._[jIj-2L-J _.1L_/_~_~ gY L-__(__ l -'Lc__JL/_[iX —- — i- J-J- -I -'_ 2T/DY __ — — ___ _ ____ ___ _IA- ----- — 3 — — _ —----------------------------------------------------------- DELTA T( I,1)=T(I 3) T (1 i,N+2 ) =T ( 3N2 ) T ( 1,1) =T( 3,1) - N'Hw EY_,N EG. P S L N 3 = E- L[N __ ___ THROUGH YY,FOR I=31,I.G.',lM+i T_______________ __H_I__ TFOR J- N I Q z II \ + 1 __ ______ ___ _ ___ _ _ _ W(I,N+2) = TM*U ( I,N+1 ) /DY W ( 2,J ) =- T','I* V ( 3,J ) iDX WHENEVER U(-,J).GE.0 47

~_ " W 9 -J w (I. -it' _F _ —---------— _ —---— i LLESii LSE_-_ ---------------------------------- B=(W( (I+1 J )-W(I,J) )/ DX END OF' CON[DITI ONAL WHENEVER V(IJ).GE.0 OTHERW I SE -------— _ —---------- C _ _-LI- UIJL -1 ~ J i I') __ _END OF CONDITIONAL D=( W(I+1,J )-2,. W(I (IJ) -1 J ))/ DX2. E=(W(I,J+1)-2.*W(I,J)+W(I,J-1))/DY2 1- U I,J)I'-JDT 6-V( I,J)7,-D T'C -— __ —--------— L tI- L+_-2- L -_, L_-L_-_ L2tJ_ ________ _ - YY WL ( 2,J ) W (2,J) NE=CO JJ NE=NE+1.. —-.........-...E.._-_F._N.F._.Li__,__ PRINT RESULTS NE, NE,' N T I,i' IME,S F ( 2 s,2 ). S r ( M S+ 2,N+2 )-,S T ( 2,2) TRANSFER TO EE END OF COND T I ONAL I=2 - _i —---------— J -+ I —_+ —----— _- ---- -- _____________________T I H 6 Q U Grl_ N_t_U --------— _+- _-_-___ __ —-------------------- DI (J)=(S1-(I-1,J)+ST(I+1,J )) *DYX- DY2*WL(I,J) F (2 =0 LM F (J )= (DI (J)+F(J -1 ))iE I (J) _________ _,,,, ___ _I d.R_ E~:_ J __EZ&R__J-_ 1o ] ] t_,_..L- 1__t__..____ -~~~~~~~ —-— 2 —--------- SF( I,H)=EI(H)-'SF(1, H+1)+F(tH) IN.T T(T H =;F(THI WHENEVER I.L.M+1,TRANSFER TO II --— _________ ________ _ ________________________________________________________________-_______________ HH J=J+1 THROUGH ViMFOR I=3,1, I.G.M+1 _I =(TI =(. T 1J-1 )+ST( I, J+1 ) I _DXY-D'X2'WL ( i,j ) F(2)=0 __J__L......... F (-I_ —- ( —— _LD+F-_- J -_L — — JI(I_-___+_ E-_!___________)._ __ ___________ THROUGIH NN,FCOR I=1,Ii.G.M-1 ---— __ —------------ L - — 4____________________- -------------- H=-iM+ 2-I NN SF (HJ)= EI( H) S F ( 1 J ) F ( J F WHENEVER JLN+1,TRANSFER TO HH -J_4AO__; __A _______________________________________ THROUGH AAA,FOR I-3,1,I.G.*M+1 WHENEVER((I-4, )/NS-i1-4)/ Nd).E.O.AND. ((J - 4. )/Ni 8 -iJ-4)/N8) 1.E. 0 R3=.ABS.((SF(I,J)+1.OE-30) /(ST ( I,J)+1.OE-30)) --— w-L-LW.U..EVNE.4_ tV. AX =________________ J 3,E A R = _______________________________ 3_______ AAA END OF CONDITIONAL ___ — WAk,____f G.L/_H__NPEY L N 3-_ —__..-_A_ X _____,_E___ - _ ___ TH ROUGH LL,FOR I=3,.,1.G.,+1 THROUGH LLFOR J=3,1lJ.G.,N+1 LL S.T (I,J)=SF(IJ)

TRANSFER TO JJ CONTINUE END OF CONDITIONAL THROUGH ALPHA,FOR I=1o,1,I.G.(M+2) ALPHA ST(I,J)=SF(I,J),,,,,_,,,,,,,,,,B___ T I-.I_~.g!U __ _ fl FI. AE_ R__ i-aL_L. j_,l.. _L _._G_ CLM_+_ __ _ THROUGH BETA,FOR J=2,1,J.G.N+1'SF (M+3,J) =-SF.( M+ 1,J) U(I,J)=(SF(I,J+1 )-SF (I,J)) /DY......_E_}LA_.....__.__.../_ I __ J_=_ ~ 2E_~I i J_G_-E_(_ I+ __ ____ ______________________________________ WriENEVER NT.L.N2,PRINT RESULTS'NT,DT,TIliE, U(2,2)... -----— _1U( —+2-N+2 -— V(2,2 — V- 2i+2 T E ---- -—.(M~2,N+2J - NE 2 SF-' (2,2).. SF (M+2,N+2),WL(2,2) i..W L(M+2,N+2),RMilAX NlTO= (NT+O. )/N3-NT/N3 WHENEV E R NTO.E.O UR'=- U ( I 1,2)/2. THROUGH DAL,FOR J=3,1,JG. (N+2) DAL UR= UR+U(IiLJ)_ PRINT RESULTS NT,DT,TI'iME, R,U(2,2 )...U (i2,N+2 ), 1V(2,2).* *.V( M +2,N+2), (1,1) r ( -M+2,N+2),SF( 2,2 )...SF(M+ 2 N +2) 2,N WL ( 2,2 ) ~ ~ ~.. WL(+ +2,N+ 2), RAX OTh ER WI SE _ — CONTINUE END OF C OQLD ITQION AL WHENEVER T2.G.2,TRANSFER TO GEEM WHE-'NEVE T ( 1,N+2).G.T3,TRANSF'ER TO LAMDA TRANSFER TO BACK _ -__ ~ zL ------— ____ __._J_7_ ~_ Z~_A..ILI_ T —_. _.__T_- _ — _ —--------------------- THROUGH AAFOR I=3,1,I.G.M+1 WHENtEVER ( (I-3 ) /N8- (1-3)/ N8) EO AND ( (J-3)'/ N8- ( J-3) /N 8) 1 * -E RI=.ABS. ( (SF(I,J)-ST( I,J) )/S1 (J-) ). L____________________-__LA I LJ,___T_ _ I __L:_ L -(_I_IPID_I_/_ L _L-J_L_- _I WHENEVER R1.G.EPSLON.OR.R2.GoEPSLN2;,TRANSFER TO BACK LAMDA 1=3 il/2 THROUGH LAMFOR J=3 1 J G (N+2)..._.S.t___________ __Ut__uL+ 31J_________________________________ __,___-____'_ __ PRINTT RESULTS NT,DTT IME,UR3,U( 2,2).. U(+2,N+2), 2,NE,RMAX -FRANSFF'R TO) TA RT ELE CONTINUE. —-------- -----— _'- -- _ —_(-____ _........... —-------------------------- $ DATA 49

40 30 20 l IFlow Pattern \ 10\ \ s wall flux = 10, Btu/hr/ftt 5\ s \ 1, 1esurface temp. = To T = 6.785 x 10'4 [I\ \ lr I 50 sec 21 x 21 grid Fig. 1 50

wall flux = 10, Btu/hr/ft2 surface temp. = To a= 1.67 x 10-lO 2 min 21 x 21 grid Fig. 2 51

G) W ON wall flux \ 10, Btu/hr/ft2 surface temp. To T=.003 F 3.6 min 21 x 21 grid Fig. 3 52

Flow Pattern wall flux = 1000, Btu/hr/ft2 surface temp. = To T = 7.07 x t051 sec 21 x 21 grid Fig. 4 53

m Flow Pattern wall flux = 10 surface temp. = To T=.003785 % 4.540 min 16 x 16 grid Fig. 5 54

Flow Pattern wall flux = 10 surface temp. - T T - 003880 - 4.545 in 11 x 11 grid /,Fig. 6,

600 l \ It | wall Lux " 4010S i Pte Streamline Pattern wall flux = 1000 Btu/hr/ft" adiabatic upper surface v = 0.00114 \ 3.42 mia 11 x 11 grid Fig. 7 56

24[ 160 ~80 1 r- I Streamline Pattern wall flux = 1000, Btu/hr/ft2 surface temp. = T sat. = 6.6 x l0-4 48 sec 21 x 21 grid Fig. 8

20~F above initial 1.0 temp. 0.9 0.8 4.35 0.7.08 0.6 -.052 0.5 -.0022 0.4.000035 0.3 Isothermals ~~~1.3 x 10-`8 1wall flux = 1000, Btu/hr/ft2 0.2 free surface temp. = Ts T = 6.6 x 100.1 S 48 sec 21 x 21 grid o 2.4.6.8 1.o Fig. 9 58

2000 21 x~l-A IY = y/a 1500 IF 2L a =.25 ft 1000 Cx/b _ 0.8 + b = 1.0 ft 500 ig 141x 16 -500 11 x 11 grid y 1. 1000 x/b = 0.6 500 0 o 1000 x/b 0.1 Effect of Grid Size on dimensionl~Velocity Profile (b) 10 x 10 grid 1000 500 ( 3_____________.6 min dimensionless distance 59

10-5 1041- /0 x / o / x/b=.6, / 16 x 16 0 o / 21 x 21 grid, 10 - / lo- lo- 1 lo-/ dimensionless time T = 2. t a2 Fig. 11 60

0.8 ~~-150 F(a) x -- - -oI. -100 1 _ \ -50 O o 0 Velocity vs. Grid Size, Ax -300 r-4 -250 (b) x.6 -200 wall flux = 10 -150 surface temp. = T.6 - 100 T.00257 5o m ~ 3 min I I l I l il l l I 0.01 0.02 0.04 0.06 0.08 o0.10o grid size Fig. 12 61

water water surfac surface heated wall heated wall (a) (b) water surface heated wall (c) Fig. 13 62

cC)'_ co - ---- C0 -i.ia___ _ C - DD 0)-