T H E U N I V E R S I T Y O F M IC H I G A N COLLEGE OF ENGINEERING Department of Mechanical Engineering Heat Transfer Laboratory Technical Report No. 4 FINITE DIFFERENCE SOLUTION OF STRATIFICATION AND PRESSURE RISE IN CONTAINERS Herman Merte, Jr. John A. Clark Hussein Z. Barakat ORA Project 07461 under contract with: NATIONAL AERONAUTICS AND SPACE ADMINISTRATION GEORGE C. MARSHALL SPACE FLIGHT CENTER CONTRACT NO. Nas-8-20228 HUNTSVILLE, ALABAMA administered t hrough: OFFICE OF RESEARCH ADMINISTRATION ANN ARBOR January 1968

TABLE OF CONTENTS Page LIST OF FIGURES v NOMENCLATURE vi ABSTRACT ix.I.. T.TRODUCTION 1 I o FORMULATION 5 Ao LIQUID REGION 8 Bo VAPOR REGION 10 Co LIQUID-VAPOR COUPLING AT INTERFACE 13 Do SUMMARY OF FORMULATION 14 IIIo TRANSFORMATION OF THE PARTIAL DIFFERENTIAL EQUATIONS 16 Ao LIQUID REGION 16 Bo VAPOR REGION 17 IVo DIMENSIONLESS FORM OF THE EQUATIONS 25 Ao WALL TEMPERATURE 25 B o LIQUID REGION 26 Co VAPOR REGION 28 V. METHOD OF SOLUTION 3 0 Ao FINITE-DIFFERENCE FORMS 30 Bo STABILITY OF FINITE-DIFFERENCE EQUATIONS 32 Co COMPUTATIONAL PROCEDURES 33 Do COMPUTER PROGRAM 39 VIo COMPUTATIONS 41 Ao GENERAL ASSUMPTIONS 41 B. VARIATIONS POSSIBLE IN PROGRAM 41 C. VARIABLES 42 Do RESULTS 43 APPENDIX Ao COMPUTER PROGRAM LISTING 59 APPENDIX Bo FLOW CHART 71 APPENDIX Co PROGRAM NOTATION AND NOMENCLATURE 73 iii

TABLE OF CONTENTS (Concluded) Page APPENDIX D. DATA INPUTS 83 APPENDIX E. TYPICAL OUTPUT 85 REFERENCES 97 iv

LIST OF FIGURES Figure Page 1. Container configuration and coordinate system. 6 2. Container wall section. 7 3. Control volume of vapor region. 21 40 Effect of grid size on computed interfacial velocity. 35 5. Effect of grid size on computed condensed masso 36 6. Effect of grid size on computed temperature distribution, 37 7. Total mass evaporated and pressure rise under high heat flux. 44 8. Total mass evaporated and pressure rise under low heat flux, 45 9. Effect of heat flux on pressure rise and mass evaporated. 46 10. Effect of heat flux on pressure rise. 47 11. Effect of initial ullage on pressure rise and mass evaporated. 49 12. Effect of initial ullage on pressure rise and mass evaporated, 50 135 Effect of acceleration on pressure rise and mass evaporated. 51 14. Axial temperature distribution. 52 v

NOMENCLATURE A -Area, Ft2 a - Cylindrical tank radius, Ft, Acceleration, Ft/Sec2 b - Total tank height;, Ft, Cp,Cv - Specific heats BTU/LBm-~F D- Substantial derivative of function f Dt e - Specific internal energy, BTU/Lbm E - Total internal energy, BTU g - Acceleration due to local gravity, Ft/Sec2 h - Heat transfer coefficient, BTU/HR-Ft2-~F, Enthalpy, BTU/Lbm hfg Latent heat, BTU/Lbm k - Thermal conductivity, BTU/HR-Ft-~F P - Pressure, psia q Heat flux, BTU/HR-Ft2 r - Radial coordinate direction, Ft; Lb f Ft R Gas constant, Lbf Ft Lbm OR t - Time, Hour T Temperature, ~F, OR AT - Tw - Tsat, ~F ATwmax - Specified maximum value of ATw, ~F u - Axial component of velocity, Ft/Sec U - Dimensionless Axial velocity- Eq. (85) v - Raidal component of velocity, Ft/Sec vi

NOMENCLATURE (Continued) V - Volume, Ft3; Dimensionless radial velocity - Eq. (85) w - Mass flow rate, Lbm/Sec x - Axial coordinate direction, Ft. X -Liquid height in container, Ft Ax - Finite-difference grid spacing, Ft Z - Compressibility factor, [1] Greek Letters - Thermal diffusivity, Ft2/Hour 5 - Coefficient of thermal expansion, OF-1 y - Ratio of specific heats Cp/CV, [1] 6 - Tank wall thickness, Ft - Dimensionless radius, [1] - Eq. (85) @ - Dimensionless temperature - Eq. (85) - Dimensionless axial distance, [1] - Eq. (85) - Dimensionless stream function - Eq. (85) f',t" - Liquid and vapor stream functions, Equations (50) and (56) [1] T- Dimensionless time - Eq. (85) p. - Kinematic viscosity, Ft2/Hour v - Dynamic viscosity, Lbm/Ft-Hour p - Density, Lbm/Ft3 W(' - Vorticity, Equation (49) [Ft-Sec]-1 - Dimensionless vorticity - Eqo (85) vii

NOMENCLATURE (Concluded) Subscripts d - discharge f - liquid or fluid g - gas or vapor i - liquid-vapor interface - liquid o - initial state p pressurization r - reduced s - saturation w - wall

ABSTRACT The processes of heat and mass transfer interactions between the gas and liquid phases of a single component in cylindrical containers with axial symmetry are considered. In the general formulation attention is given to the cases of external pressurization with and without liquid discharge as well as to the nonvented condition. The governing equations are cast into finitedifference form and numerical computations are carried out for the case of a non-vented container having an imposed head flux, using ideal gas relations for the vapor properties. Of specific interest is the calculation of the pressure-time history of the container under these conditions. ix

I. INTRODUCTION A number of space missions of current interest will require the storage of liquid propellants for long periods of time, varying from hours to months. The ultimate goal is to maximize the quantity of useful propellant remaining at the end of this period of time. In the limit this maximum corresponds to the non-vented condition in which the original mass of propellant is retained. Whether or not non-venting is practical depends on the maximum internal pressure which results as a consequence of the thermal interaction between the storage container and its ambient. Higher pressures require increased structural weight of the container, and a trade-off may become necessary between tank weight and propellant loss due to venting in order to maximize the mass of propellant remaining. Each case must be evaluated for the specific mission, storage time, environment and propellant. An additional factor is the method of discharge from the tank. Where the propellant discharges to the suction of a pump, some degree of subcooling is necessary to avoid cavitation and, depending on the fluid, may require additional external pressurization just prior to discharge. The pressure at which venting might take place would thus be less than the tank design pressure. The pressure in a tank containing two phase, liquid and vapor, is related directly to the temperature at the liquid-vapor interface for a single component system. For a binary system, the pressure depends not only on the temperature at the interface but the relative liquid-vapor concentration of the two components as well. Pressure changes take place owing to heat and mass transfer interactions between the vapor, liquid and container walls. The processes of heat and mass transfer interactions between the gas and liquid phases of a single component in cylinderical containers with axial symmetry are considered in this study. In the general formulation presented first attention is given to the cases of external pressurization with and without liquid discharge as well as to the non-vented case. Solutions are then presented, utilizing numerical finite-difference procedures for the non-vented case. The initial conditions of the liquid and vapor are considered to be those at equilibrium, with uniform pressure Po and saturation temperature To From these initial conditions, the walls of the container undergo a thermal perturbation, such as a change in temperature or an exposure to an external heat flux, either of which may be an arbitrary function of time and axial location. The perturbations in the boundary conditions lead to a series of non-equilibrium phenomena within the container. Natural convection currents are set up in the liquid and in the vapor spaces. The liquid-vapor system 1

tends to adjust to the new nron-equilibrium conditions within the container by t:rans~ferrzing mass and energy across the i.nterface by either evaporation or _condensation. The cond;tions at the liquid-vapor interface couple simultaneous tr ansport processes in the liq'id and gas phases. In the case of self-pressurn zation of non- vented tanks the rate of pressure rise within the tank is governed by -Thele raite of heat and:mass transfer frorn the walls anid the liqui-d, phase to the gas phaseo The interfacial temperature is essentially that of equilibrium (saturation) conditions corresponding to the system pressure. At t+he same t-ime, the temperature of'the liquid-vapor interface affects both the interfacial mass and heat transfer, as well as the convective processes in both phases. These latte-r processes influence the temperature gradients within botUh phases arnd- JIn turn will have an effect on the rate of pressure rise in the ullage space. Tndeed, all the processes of heat transfer from the ambient -o both phases', the natural convection within the container, the interfacial ph'enormena and'the rate of Lullage pressure rise are all mutually coupled. Such -int-eractnons have been the subject of many experimental investigations. An-y anralytical app-oach that would adequately describe the phenomena tak:ing place in propellant tanks must take into account these interactions. ThLs requires the calculation of the transient velocity and temperature profLles in both. the gas and liquid phases, as well as the corresponding concentration distributi-Tons in the case of multicomponent systems. Owing to the complexity. of these ph.enonerna, ne analytical studies have been presented which consider hei.r simultaneous interaction. The problem nmai be complicated further with -' the presence of turbuLence or boiling of the liquid near the tank walls, as the1ir nLatucre is not yeit fully understood nor adequat;ely described. It is evident t1ha't some assumptions are necessary for the construction of models which reasonably represent the practical situations. A 1-eview -whCtich covered mfuch of the available literature on pressurizaticn, n strati:;fication and interfacial phenomena in propellant tanks is given:e.n Ref. (1)o Therefore tno further survey of the literature will be presented here. However, several studies which are pertinent to the present one will be mentio.nedrie. Thomas and Morse (2) and Knuth (3) have considered the phase change of suddenly pressurized single-component liquid-vapor system. Yang, e'~ alo (4 5,6) have solved the phase change in a suddenly pressurized sringle-component and. binary liqui dvapor system. The model adopted in these cases was one-dimensional, with the origin of the x-axis at -the initial position of the liquid vapor int-erface. Initially, the gas arnd the liquid regions are at the same temperature To. The conditions in the gas space are suddenly changed to a pressure P0 and temperature To. The differential equations describing the temperature and velocity transients in both phases are: aTg aTg a2Tg aT. aT ( at aQ 8x2 (2)

These are coupled at the interface by the condition pa hfg U = k -x = g Tx Thermodynamic eqU-ilibriumrn is assumned at -the interface, i.e",'the irnterfacial t+emperature'is always at the saturation temperature corresponding to the ullage pressure. Experience shows that th:is is a reasonable as-mpti.on (1)>'The details of the sol-;tio-n as well as the differential equat:ions describing t+he oncer:n'raion distribution for the binary case are given in Refl (1). The effect of the changes in the properties of the fluids with temperature on the accuracy of the results was invrestigated by' OVLoughlin and Glenn (7)o The same cne-dimens1orial single compoment mno;el was employed' to study the case of iInt'erfacial phase change with variable ullage pressureo Finite-differences were used'to solve Eqs. (1) through (3). The pressure in the ullage vol.mzie was either constant or a specified funcili.eon of timeo For the case in which the press'ure decreases with time, the interfacial temperature also is reduced with time. In this case'the liquid layers which may exist at a temperature highe t haLn the interfacial temperature are made equal to the saturatncn;.emperature. The excess energy' of this li.quid over that of th e thsaturate. li.quidi'is used to evaporate some liqui.i int.o *the vapor reg-i.o:nr Likewise any supersat.uXrated. vapor that may be present is allowed to condense a fract-ion of t-he vapor, thus provii;i:ng t;he e:.etergy requ.ired- to heat the supersaturated vapor cto the saturation temperatjure. This supersaturation is a result of the expansion of t —he gas during the pressure decrease process, which is assumed t' be isentropic. The resul'ts show {t;hat the variation rin the -thermal condauctivit;y wit'lh temlperature h.as very little effect. Thi-s one-dimensional model does not; acco'un for the effect of wall 1hea`t transfer, Epsitein, e't al., ( 8 ) u se fini eiifferences for the calculation of the praec, oia~tozatoprI pr.ess... t'this model, the axial variation of temppera't;ure is co n. sird G..ble th, he rad-ial variatio n is neglectedo The wall' to fluid heat' jrart'.fer was acco-unted for by iEnt:rodcielng an effectivl.we heat, transfer coeffi.ci ent be~tween t;he fluid. amnd the wall. Al. so, effectiv.e'themal crnlductivities ant mass diffusivities between th.e adjacenut flu.id laye'rs we-re uIsed o sffirmulate thhe effect)sL of randiom fluid motion. Pro'vi s ions for variable tarik cross= sec ti:on as well as variable heat talansfer coeff.Lcients were made i. the progrcm.o'The effec-t orDf boiling near the walls was not considereed, and momientum effects were neglected. The principal advantage of this program s s its simplicity anA the relatively small. amount of machine t-me requi.red to advance thke soluition to a particular time level. However, the determinartion of the effective ther~mal conducrtivities anid heat transfer coefficients and their variation with locatioi require the dete r:mination o. a large n0 bei of erperical constants, whilth TrquirS cosisde&rable experimen-tal experience wi.+th a pa-rt;icular SyStemS 35

Althoughsome analytical models have been used to study liquid stratification in self-pressurized propellant containers, none has considered the simultaneous interactions between the liquid and the vapor phases. Most of the available studies of the pressure and temperature transients in the vapor space are experimental in nature and have served to identify some factors influencing the rate of pressure rise within the containers. The studies of Huntley (9) show that liquid stratification causes the pressure to rise at a higher rate than that calculated using the average or mixed mean liquid temperature. These also showed that stirring the liquid causes a smaller rate of pressure rise, while stirring the vapor resulted in a substantial increase in pressure. Higher rates of pressure rise were obtained with smaller ullage volumes. Liebenberg and Edescuty have shown similar results (10).

II. FOJRMULATION As mentioned earlier, the formulattion of the most general case will be given first. This involves external pressurization, liquid discharge and the consideration of the vented and nonvented container. The heat capacity of the container side wall is also included in the formulation. The heat capacity of the end walls of the cylindrical container is neglected. The assumptions and limitations will be mentioned as appropriate and the resulls and problems encountered in the solution of the unvented, non-discharge case will be given. A cylinderical container of diameter 2a, height b and with wall thickness 6, is partially filled with a liquid as shown in Figure 1. The initial heightc of the liquid is Xo, and that of the vappor is b-Xoo The origin of the coordinate system is taken at the centser of the corntainer base with x-positive in the direction of the liquid~ At any time t, the location of the liquid-vapor interface is given by X = X(t). The fluid in both phases is initially at rest in an equilibrium state at a uniform temperature To and pressure PO. From these initial conditions the tank side walls are subjected to an arbitrary heat flux, qw(x,t). The tank ends are assumed adiabatic but this restriction can be remnved if desired. The differential equation describing the temperature-time histeory of the tank wall obtained by lumping the wall radially, but not axially, is, Figure (2),' w T f + (bk x2T w (pcp6 -(k -) + PPat ( ar ax2 + qw(xt), (4) where the subscripts w and f refer to the wall and the fluid, respect;ively The first term on the right-hand side of Eqo (4) is the rate of heat flow from the wall to the fluid, while the second term accounts for the net axial conduction is usual to calculate the heat flow from the wall to the liquid by introducing the heat transfer coefficient h according to; (k =T) h(Tw - Tf). () However, the determination of the heat transfer coefficient for these cases is difficult and subject to many uncertainties. Furthermore, its a priori use in studies of this kind is not relevant~ The nu merical procedure used in 11S study allows the direct evaluation of the rate of heat flow from the wal1 to the liquid, thus making the use of a heat transfer coefficient unnecessary. 5

9 ~uIsXs aWuTipaooo pui uoTq.qzrjjuoo aauEquoD lT awnrtj, ~1&'j~ ~ I!i_ I~~~~~~~~~~~~~~~ q e'm, I.II

I-I: x 3: I I I" —— ~VZt — qw (x,t) kf(d5 )f 33 _I Figure 2. Container wall section. 7:

A reversed process may be used, however to determine a value for the heat transfer coefficient under these transient conditions. From the rate of heat flow from the wall to the fluid in the container as determined by Eq. (5), the transient heat transfer coefficient h may be computed. This procedure permits the determination of -the variation oJ the heat tra3nsfer coefficient with t;ime arnd space during th:is transient processo The differential equations governing the velocity and temperatu;re distribution in both phases are developed as followso A. LIQUID REGION The following assumptions are made: 1. Constant thermal conductivity and viscosityo 2. Incompressible fluid. Density variations are introduced only in the body force term of the momentum equation. These variations are described by P = po[l + (To - T)] (6) 3o The influence of viscous dissipation is neglected. Governing Equations (i) x-momentum P U+ _ uP Pg ax + 2 1 ( + ~_ ) ( (ii) radial moment'um +t X + + ar t 2 I IA-v +V (8) 8

(iii) continuity +u + av = o (9) ax 6r r (iv) energy 6T + T + T 62T 1 6T 62T - +u- +v + + r2(1 Initial Conditions T(xro) = To, u(xro0) 0, v(x,r,o) = 0 (11) Boundary Conditions Velocity u(ort) = (12) u(X o,r,t) - (13) where wid is the rate of liquid discharge~ u(x,a,t) = O (14) bU'E (,=2 O ( 5) 6r v(Xa,"t) 0 (16) v(x, o,t) O ('7) v(o,r,t) = (18) v(X, rt) = 0 (19) 9

Temperature T(X,rt) = TS(t) (20) T(x,a,t) =Tw(xt) (21) 6T(o,r,t)_ O (22) 6T(x,o,t) 0 ( x,. og (23) br where Ts(t) is the saturation temperature corresponding to the system pressure and Tw is the wall temperature, Eqo (4)o Bo VAPOR REGION The differential equations for the vapor region are obtained assuming negligible changes in the viscosity and thermal conductivity within the vapor phase. The compressibility of the vapor is considered. Governing Equations (i) x-momenturn Pu +'U U + - pg a + ~ 2u +1 a + 2u 3 ax (x (24) (ii) radial momentum + \2 - r2 r ar a-2/ r x r + (25) 10

(iii) continuity Dp + + v + V O (26) (iv) energy PCv + u a + v + P DP + k 2T+ 1 6T 62T'2 Dt (0>7 r ar a-2) *(27) (v) equation of state P = zpRT (28) where z is the compressibility factor which can be evaluated from P. V, T, data or by using an equation of state which accounts for real gas effectso Using the Van der Waal equation of state the following equation relating z to the reduced pressure and temperature is obtainedo 5(3 (Pr+1,2 27 Pr 27 pr2 K8Tr 64 Tr 512 Tr3 where Pr and Tr are the reduced pressure and temperature, respectively~ The value of z obtained from Eqo (29) is quite accurate for engineer:ing applications (11)o Other equations describing z or tabulated values of z can also be employedo However, calculation of z using an equation of the form of (29) offers less difficulties from the programming standpoint. Initial Conditions T(x,r,o) = To p(x,ro) = p u(x,r o) = O v(x,r,o) = O (30) 11

Boundary Conditions Velocity u(b,r,t) 0 (51va) or u(b,r,t) = ugp pgpa2' X (31-b) where Wgp is the rate of mass flow of the pressurant and pgp is its density. Equation (31-b) is written for the case of uniform velocity at the inlet diffuser of the tank as in Fig. 3o u(X+o,r,t) = ugi (32) au(xot) = O(33) 6jr u(x,at ) O (34) v(br,t) = O (35) v(xa9' it) 0 (356) u(xo,t) O (3) ~v(X~rt 80 (38) where ugi is the velocity of the vapor at the interface caused by simultaneous phase change and liquid discharge, Eq. (47). The boundary conditions givenl by Eqs. (19) and (38) assume a zero shear stress in both the liquid anid vapor, at the interface. These are reasonable and are made for the sake of simplifi.cation. Another approximation would be to take the vapor velocity at the interface equal to that of the liquid at that point. The interfacial shear stress would then be determined by the analysis. For flows wit;h large radial gas velocities at the interface2, this latter boundary condi-tilon becomes mof e realistic. 12

Temperature T(X,r,t) = Ts(t) (39) T(x,a,t) = Tw(x,t) (40) aT(x,o,t) = (41) ar 3T(b,r,t) O (42-a) (42- a) or T(b,r,t) = Tgp(r,t) (42-b) Equation (42-a) applies for the case with no external pressurization, and Eq. (42-b) for arbitrary external pressurization. C. LIQUID-VAPOR COUPLING AT INTERFACE The processes in the liquid and vapor regions are coupled by Ts in Eqs. (20) and (39) and by ugi in Eq. (32), which is related to the rate of mass transfer at the interface. The rate of mass transfer by evaporation or condensation across the liquid-vapor interface depends on the relative rates of heat transfer by diffusion from each phase at the interphase. Conservation of energy at the interPface determines the rate of phase change as well as the interfacial displacement, according to; hWi = oa [ r[k (XT(X r,t)] [k IT(' )]v 2Ardr (43) hfgw. = f {[k ] [x ~x r where w. is the rate of interfacial phase change. According to Eq. (43), w. will be positive if condensation takes place. The interfacial displacement with no liquid discharge is then given by: p ta2 dXi =w (44) 13

Should it be desirable to include the process of liquid discharge, Eqo (44) would be written Pi ~a2 dXd wd (45) dt U d dt = tldd (46) where dXi --- = rate of interfacial displacement due to phase change, dt dXd d-d — rate of interfacial displacement due to liquid discharge, dt dX dt; = combined rate of interfacial displacement due to phase change and liquid discharge. Equation (46) is obtained assuming that the liquid surface rema.ins flat during the discharge processo This is an approximation which neglects the influence of viscosl.ty near the walls and surzface tension effects. Such effects are negligible except for very low gravity levels and small Bond:-numberso The vapor velocity ugi at the i:nterface is related to the rate of interfacial displacement due to phase change and discharge by, U PPgs- d.Xi'' d.Xd ugi = d i +Pgs dXd (47) Do SUTVMRY OF FO:RMUTATION Taking a broad perspective of the foregoing formrulation'the following observations might be made0 For convenience the liquid and vapor regions are considered separately, recognizing that they are coupled at the liquid-vapor interface o 1l Liquid Region Four urnknowns exist: u,v, T, P, each of which are functions of x, r, t The four Eqso (7)-(1)), together with the boundary conditions, are available for the solution of these unknownso The solution of this system of equations 14

by finite-difference procedures has been attempted with no success. The reason for this is believed due to the extremely small variations in pressure within the domain, P(x,r,t), compared to the level of pressure. With the assumption of an incompressible fluid (except for changes in density due to temperatures in the body force term of the momentum equation) the pressure terms occur only in the x-and r-momentum equationso The pressure distribution is not of particular interest except that it is a factor in the fluid motion. This motion is influenced by the system pressure level only indirectly in that it establishes the temperature boundary condition at the liquid-vapor interface. It is, however, possible to eliminate the pressure term by cross-differentiating and combining the two momentum equations. Introducing the definitions of vorticity and stream function, the system of 4 equations with 4 unknowns is reduced to a system of 3 equations in 3 unknowns, temperature T(x,r,t), vorticity CD(x,r,t), and stream function I'(x,r,t). This procedure will be outlined below, and has been successfully applied to a single uncoupled domain (12,16). 20 Vapor Region Since the assumption of an incompressible fluid is not a valid approximation in the vapor region, five unknowns exist: u, v, T, PF p, each of which are functions of x,r,t. The five Eqs. (24-28), together with the boundary conditions, are available for the solution of these unknowns. As in the case with the liquid region, the pressure distribution is of interest only insofar as it gives rise to the fluid motion. The influence of the spatial pressure variation within the vapor space on other properties can reasonably be neglected. Thus, the pressure term can be eliminated from the two momentum equations by combining them, and introducing the vorticity and stream function. The temporal variation of pressure, however, still remains in the energy equation. As will be demonstrated, the system is now described by 4 equations with 5 unknowns, i.e., T(x,r,t), D?'(x,r,t), l"(x,r,t), p(x,rt), P(t). The additional needed relation is furnished by the energy equation written for the entire vapor space as a control volume. These procedures will now be described. 15

II o TRANSFORMATION OF THE PARTIAL DIFFERENTIAL EQUATIONS The momentum and the continuity ectations are combined to obtain the vorticity equation and an elliptic equation relating the vorticity and the stream function. This is accomplished in the same manner as described in Refo (12)o The x-momentum is di.ferentiated with respect to r, the r-momentum equation is differentiated with respect to x and the two combined to eliminate the pressure termso~ Ao LIQUID REGION Equations ('7) and (8) reduce to: 1 ____ - 6 + V N ____ (48) 6at u 6x v r pr r L[aX2 + r ar r2 j where (' is given by: 1 (49) Introducing the stream function fr', defined by u = r; v - (50) Equations (10), (48), and (9) are transformed respectively to: aT +1 V' a6T 1 a* 6T (2T 1;T _ 62T> +at r 2 r r /x - r r r (51 ) at ar a6 x r ar ar r ar r = r r r 6r a2 ax2 -r r + Wr2 (53) The initial and boundary conditions zmust also be transformed appropriatselyo Equations (51), (52), and (53) are solved numerically by finite differences to obta~in the temperature, vorticity and strean functior distributions0 From 16

Eq. (50) the velocity distributions can then be computed. B. VAPOR. REGION Combining the momentum Eqso (24) and (25) in the manner described above) the following is obtained: ) 3' t o'')' Dp 1 Du 6p 1 Dv ap t +u +U V7 - P Dt pr Dt; T + p-r DT Pr r + 7 (54) where w(, is given by Eqo (49). Equation (54) is the vorticity equation for this case. The first terms within the brackets on the right-hand side of this equation account for the compressibility effects due to density variation. Should the compressibility effects be neglected, Eq. (54) will reduce to Eq. (52). It should be noted that the presence of these terms does not introduce any additional difficulties as far as the solution of the vorticity equation is concernedo The difficulty will be in combining the continuity Eq. (26) and Eqo (49), in order to ob`tain an equation similar to Eq. (53). However, this difficulty can'be overcome if the term ap/1t is neglected in combining the continuity Eqo (26) with the definition of vorticityo The rigorous jus tification for this approximation has not yet been established. In this case EqO (26) will be rewritten as: 3(pu) _ 1 (prv) 0 (55 ax ar The continuity Eqo (55) is combined with Eq. (49) by introducing the stream function *K, which satisfies (55). This stream function K" is defined by: u p-; v' - p - (56) Su;bstituting (56) and (49) into (55), the following equation which relates i" to 0w', is obtained, ___"_ 1 It 2\ 2 Lv (57+2) 17

The last terms within the brackets in Eqo (57) account for the density changes in the vapor region. The terms ap/)r and ap/)x are calculated in the finite-difference procedure using the equation of state and the temperature distribution from the previous time stepo The coefficients u and v in u ap/ar and v ap/ax will be taken equal to their values at the previous time step, which is the same procedure used in calculating'the non-linear -terms u aT/)x, v aT/)r, u 1)'/)x, and v -I~'/)r, thus in effect linearizing these terms. The spatial variation of pressure in the vapor region can be neglectedo Hence, DP dP dP- = — ~~~ (58) Dt dt The density terms in Eqso (54) and (57) are evaluated using the equation of state (28), and the derivatives given by: ar RZ2T2 )T - r (59) x RZ2T2 T )X (6) p Substitution in EqO (54) gives: De w' dP WY' (Z T) DT Dt; P dt ZT T Dt g (ZT> aT 1 aT (zT rZT aT p )r rZT ar Tp Du 1 (zTa) T Dv ODt ZT a x Dt + v Lx2 + r ar ar2a (61) Substitution of Eqso (28), (59), and (60) into EqO (57) gives: 3r 2 -r x2 RZT — l) r2 R 2TP2 (T.) T v — 18

The terms u and v in Eqso (61) and (62) could also be expressed in terms of t" by Eq. (56). Substitution of the equation of state (28) into the energy Eqo (27) gives: LCv+ (KT pi DT dd + kV2T (63) Taking the thermodynamic relation Cp - Cv T ( (64) and the ratio of specifie heats Cp (65) and substituting into Eq. (63) and rearranging, another form of the energy equation is: DT 1 0 RZT dP + T (66 Dt (7I): z cVP dt +_-1)z /6( - -v (ya) V It might be noted that the body force term in Eq. (61) could also be expressed as: rZ TP (6 7) Since the assumption of ideal gas behavior is adequate for the vapor, * Eqs. (61), (62), and (66) reduceg respectively, to: De' ( J' dP w)' DT g aT 1 6T Du D- - P dt;+ T Dt r-T ar rT Dt 1 i? Dv + (6(w? 3 6 (68)6w? rT ax Dt r r 68 *For saturated 02 vapor at atmospheric pressure, Z = 0o97, and [a(ZT)/~Tp1o015o For an ideal gas, Z = Lo00, and a(ZT)/aT = 1o00o 19

2~,,, br1 ~" a2,, P ar2 - r Or aX2 = RT r2 Pr I T PRT- - (69) DT = _1 dP + 2T (70) Dt P'dt + Equations (61), (62) and (66) (or (68)-(70) for ideal gases) together with the transformed boundary conditions constitute the system of 3 equations with 4 unknowns for the vapor space, w'(x,r,t), r"(x,r,t), T(x,r,t), and P(t) referred to earlier. The finite-difference solution of Eqs. (61), (62), and (66) for w'I, 4" and T requires that a function for dP/dt be available. This is obtained from the First Law of Thermodynamics written as an instantaneous rate equation, taking the vapor space as the control volume. Referring to Figure 3 this formulation is written: d fVg(t)(Pgeg)dV - Ai hgsPgs(ugi - ui)dA + hgp gpugpdA = qdA - P V (71) Ap gp Pgp~gp Ac.s. dt This formulation includes the following generalized variation: eg = e(x,r,t) Pg = p(x,r,t) Ugi = u(r,t) hgsPgs are uniform over the liquid-vapor interface ui is uniform for a flat interface hgp,pgpugp are specified by the particular pressurization process. q is the local heat flux on the control surface, and its integral includes that from the tank wall and through the liquid-vapor interface. Given property information on eg = fl(P,T), Pg = f2(P,T),. hgs = 3(p) Pgs = f4(P), the desired term dP/dt can be extracted from the first term of Eq. (71). The necessary properties can come either from specific heat data plus an equation of state, or from tabulated values. It may be anticipated that different computational procedures would be required for each of these. 20

~ wgp P gp A CONTROL IV SURFACE VAOR qwg q Ugi1 qwQt.. dX w.d +wi Wd. Figure 3. Control volume of vapor region. 21 —,

The procedure will be demonstrated for the case of self-pressurization, that is with no liquid discharge, no external pressurization, and a flat interface. Thus, wd = O, wp = O, and ugp = O0 Then,. u dXi (72) dt dt and from Eq. (47), assuming a uniform vapor velocity across the interface Pi ugi - ui P ugi (73) P I-'Pgs For the coordinate frame selected dV dX dt i dt) For brevity qtOt = fA q d A (75) Substituting Eqs. (73)-(75) into Eq. (71) Vg(t) (Pgeg)dV q + hgsPgs pig -p Ugi Ai + P Aiui (76) Using Leibnitz's rule for differentiation of integrals having variable limits the first term of Eq. (76) can be written as d ~ Igt(ggd P Jeg dt fV (t)(Pgeg)dV = fVg Pg t dV 6Pg + fVg eg atg dV - PgsegsUiAi (77) Substituting Eq. (77) in Eqo (76) and rearranging with the use of Eq. (47) gives, 22

V9'dV V eg ttdV tt + hgs Aiui (pgs P.) (78) For the case in which the use of an equation of state (28) is desired, the entire left-hand side of Eqo (78) can be expressed as *the following by means of general thermodynamic property relations: Cv dP dV PCv 1 R dt Jv 6 (ZTR p + TT(JTL ZT dV p h9 F(zTA R vg z2T2 T v (-;9) The enthalpy terms in Eqs (78) and (y79) can be expressed by: h-hTR Cp dT + /PR F-(ZT) dp hhR JTR PR d + Z L T p (80) where the subscript R refers to a low pressure reference conditiono For an ideal gas with constant specific heat expression ('79) reduces'to: R- )vg (81) Substituting expression (81) for the left-hand side of Eq, (78), writing the enthalpy by a specific heat term, and rearranging, Eqo (78)'reduces to: dP _ R RTs dt CvVg qtot + T Aiui (Pgs p,) (82) d~t Vg For the case far from the region of the critical state, such that pi >> Pgs) Eqo (82) can be written as dP R wiRyT s dt, CvVg oQot Vg (8 2j5r

where wi is the rate of interfacial mass transfer (having a negative value for evaporation) and is given by PRPgs Wj (PPgs) Aiugi (84) 24

IV. DIMENSIONLESS FORM OF THE EQUATIONS The governing equations for the wall, liquid and vapor regions (4), (51), (52), (53), (61), (62), and (66) are next made nondimensional, along with the boundary conditions. The substitutions necessary to nondimensionalize the differential equations are: c~b 5~ v~G~b u = U; v = V; T-To = 4 a2 t a; x = b; r a =T c;' r'br (85) The resulting dimensionless equations are given below. Ao WALL TEMPERATURE The dimensionless equations describing the wall temperature, from Eqo (4), are given by: x(t) (i) 0< ~5t b dQw a (PCp)~ 9 a2 (PC0 kw _2_w pw ( 2 aT - W (PCp)w (P) pWw + a a (P, Pr Gr*(,t) (86) b 6w (PCp)wW( (i) bx(t) < < 1 b aGw a. (p.Cp) k g o - _ a2 (pc) kw a26 a s (pcp)W k- k' -'n'+1 b2 (pcp)~ k- a. 2 p apw a a (pc ~ Pre Gr(~,t) (87) 25

where Gr*(1,t) is the modified Grashof number and is given by: w G:r*(,t) -- gha (88) W k.~ I, %, (iii) Boundary and Initial Conditions T7) "0 (89) o(3,o) 0 (90) Bo LIQUID REGION (i) The energy equation, Eqo (51) a0 1 ~a ao 1 ai aG a2 2 1 ao a2g —,s), -- *-(91) (ii) The vorticity equation, Eqo (52) a) 1 a o 1 a Ca c a 2 a23- a~ aO21 + —~ ~ P- 2 (9rl2Q + a )2u2 at rl a) r a e an ) Pr5 b2 a2 rl a 2j (92) (iii) The vorticity-stream function equation., Eq. (55) ~ a - = a +. c (93) b2 )2 2;3q U -; V (94) (iv) The Boundary Conditions (1) Stream function boundary conditions,(0 7 1 T)= — =O o- 0 (95) (? (= 2 0O (96) 26

t(,o, 0-T) - ll ( i 2fl. (97) T( n,lT) = Av(,l, T) 0 (98) (2) Vorticity boundary conditions )(1,I, T) = O (99) W(,O, T) = O (100o) Two additional vorticity boundary conditions are required at the tank wall and bottom. An explicit expression for the vorticity at any of these locations is difficult to obtain. The method of solution used in this work, which follows that of Ref. (16), overcomes this difficulty. (3) Thermal boundary conditions e(,i, ) T= G,r(:,T) (101) = GS(T) (102) 8G(0fT) = 0 (103) wae(eoe) i e (4) Initial Conditions = G(,T,o) -= T(,,0o ) = 0 (106) 27

C. VAPOR REGION (i) The energy equation, Eq. (66) La 1 a ag 1 a__ ag 1 RZ (+gOo) aP 6T -P1 7 -- 097 16(zT) Cv P T 1++((1)Z( a(zT Tv _+ __Y_ _g _a2 a2_ 2 s2g (107) ((T)) ~ Lb 2 n arl n2 where 0 = g a4 T0 (108) (ii) The vorticity equation, Eq. (61) a 6 1 a a 1 aPr a vg 2 3 aC. aL 1 2 aD 4 j(ZT' a aTPri a a 7 P q v 22 P rl ai C PbZ1g o 1 (a(zTh ag Du 1 ( (ZT Y DV WdPo'qz(o+@) kT a DT-rZ(G+Q0) ~o T 7 -DT \dT) (109) where P is given by: 1 ((Z) (110) (iii) The vorticity-stream function equation, Eq. (62) a2 62* 1 6* _2 P p __ 2 T2 F rP aq + F6 ZR(- ( (+0 v2 ab~ZaR(G (T@ ) aT -L a7 V. p _~T)a2 F 0a'6 (iv) Boundary Conditions (1) Stream function boundary conditions ( ni T) -= =(1.9T) = o (112) t ( ti,71,T)- (Pt -Pg.S) vi 2L PgsUgi _l2( 28

( T) = 1 24(,O, r) -0 (114) =50 T) ~Pgj 2r1 (,:.1,-T) = 2(1r, T) = o (115) (2) Vorticity boundary conditions (t go07) = O (116) ~(1)(,.-T) = 0 (117) (3) Temperature boundary conditions @(Siyrl7T) = S(T) (118) - T) = o (119) (, _,.-T) = oW(, -,T) (120) (.k0, T)2. = o (121) (4) Initial Conditions (,.,o) U(=,n,O) = V(S,,o) -= w(~,,O) = O (122) From the above results, it is clear that the temperature, pressure and velocity within the container are functions of the parameters Gr*, Pr, (pCp)/(pCp)w, (kw/k~), (g/ (g/vj) (Pg/) /), (a/b) and (a/65) Gr* is defined in Eq. (88)o 29

V. METHOD OF SOLUTION A. FINITE-DIFFERENCE FORMS The finite-difference method of solution used in Refo (16) is adopted here. A complete discussion of the application of difference methods for the solution of the energy and vorticity equations is given in this reference, along with the problem of stability of the difference equations. A brief discussion of the method of solution will be made here. The basic concept in the application of finite-difference methods for the solution of partial differential equations is the use of Taylor Series Expansion to approximate the derivatives at a point in terms of the value of the function at that point and/or at its neighboring points. This may be demonstrated as follows: o1 The time derivative is represented by af f(T+LT) - f(T) + O(AT) (123) dT AT 2. The first order derivatives a and can be approximated by 3f f f(x+Ax,RT) - f(XR,T) + O(Ax) (forward differences) ~ax nAx (124) or ( f =f(xR T) - f(x-AxR -T) + O(Ax) (backward differences) x(i ) = ~Ax (125) or af = f(x+Ax,R,T) - f(x-AX, R,T) + 0(Ax)2 (central differences) (i) x 2Ax (126) 35 The second order derivatives are replaced by finite-differences according to the formula. 30

a f =f(x+Ax,R, T) - 2f(x,R, 9T) + f(X AX,R, T) + (Ax)2 (127 6X2 (ax)2 + O( x) ((27) The function f represents either C or wn, and Ax is the spatial increment in the x-direction~ The lastt;erm on the righthanrd side of Eqso (123) through (127) indicates the order of t;he truncation error involved in replacing the derivatives by finite-differenceso It is clear that the central differences offer a better representation of the first order derivative f than the forward or backward differences. However, the form used to appriximate the first order derivatives is usually determined by stability considerationso Similar formulations can be written for the derivati;ves f and 20 The substitution of the above formulae in the energy and vorticity equations produces a set of explicit difference equations~ However, if the values of the function f in Eqso (124) through (127) are taken at the time level T+AT instead of being taken at time level Tr the resulting finite-difference methods may require the use of small time increments and consequently large machine time. Certain implicit formulations may permit the use of large increments0 The application of both explicit and implicit methods to the present problem has been extensively investigated in Refo (16)o It was concluded that the lack of explicit boundary conditions for the vorticity at the solid boundaries prevent the use of large time increments, ioeo, implicit methods0 Therefore, it was decided to employ explicit methodso It is clear from Eqso (124) through (126) that more than one explicit finite-difference formulation can be constructed for each of the vorticity and energy equationsO The finite-difference formulation chosen for the solution of the present problem is dictated by stability, as well as practical considerations, which will be shown below in studying the stability of the finite-difference equationso The method of solution used in the present problem can be summarized as follows: 1o The time derivatives - and are approximated by Eqo (1235) 20 The nonlinear terms U — V —, UN- and V - are linearized by consider= ing the velocity components U and V to be kno and are taken equal to their values at time level Tr The order of the error introduced by this linearization. can be obtained from Taylor Series Expansion0 If UO and U are the values of the axial velocity component at time levels' and,T0+ATo, respectively, then 30 3o 0u 30 ua = U +aT( a) + a (128) 31

where 0 < K A T The last term on the right side of Eqo (128) represents the linearization error and is of order O(ATr) 3o The nonlinear terms U ~ and U - are approximated by backward differences, Eqo (125), if the coefficient velocity U is positive and by forward differences, EqO (124), if U is negative The same procedure is followed for approximating the terms V and V acceording to the sign of the velocity componrent VO 4[ Central differences are used to approximate the first order terms 1 a6 35 __ and 1 6ao No stability problems will be encountered in this gad' RA~Rthe centerline where both R and are zero, 1 is replaced by its limit according to Limit 2 RO F R2 R O0 (129) 5, The second order terms 2 2 a a re represented by Eq. (127)o xx2 -R2 x- - 6o Although the first order derivative -I in Eqso (86) and (87) can be represented by any of Eqso (124), (127 =1 or (126), the following formula, which has a higher truncation error is used, al al= 6A-9 (130) S imilarly 8@2 - llQ( gl yg-2AtgrP>2G( r t-5A8grj) (151) (131) B0 STABILITY OF FINITE-DIFFERENCE EQUATIONS The stability of the finite-difference equations is an important consideration in establishing the size of the grid and the time steps, and the type of differencing used, The details are given in Ref0 (16)o As might be anticipated, the size of the grid and time steps will also be governed by the storage capacity of the machine and limitations of time and cost; The necessary requirements for stability are given by the following: 52

2a2 2.. AT a + 2 2 <+1 (132) 2a2Fr 2 Pr U..1 IV AT + + + At < 1 (133) For Prandtl number less than unity, inequality (132) is more restrictive and therefore should be used. For Prandtl numbers greater than unity inequality (133) must be used. C. COMPUTATIONAL PROCEDURES The sequence of steps in establishing the numerical calculation is as follows: (1) A convenient grid size is selected, dictated by the machine storage capacity. (2) A suitable time increment is chosen. This may be altered during the course of computation as necessary to maintain numerical stability. (3) The temperature distribution is computed using velocities, temperatures and pressure from the previous time step. (4) These temperatures are used to compute the vorticity at the interior nodal points for the current time step. (5) The stream function is computed at the interior nodal points. (6) The vorticities at the solid boundaries are calculated. (7) The velocity components are calculated. (8) The rate of phase change at the interface is determined using the computed temperature distribution, with Eqo (43), (9) The pressure rise is computed with an equation such as (83) using the parameters from the previous time step. (10) The above procedures are repeated successively. A problem inherent in the use of numerical methods is the accurate determination of temperature gradients. In the present application temperature gradients are computed in the liquid and vapor at the liquid-vapor interface to determine the rate of phase change, from Eq. (43), and in the liquid and vapor at the wall to determine the heat transfer to the bulk liquid and vapor, respectivelyo The general problem of the effect of grid size on the interfacial heat and mass transfer is presented in some detail here. The formulation with Eq. (43) was used by a number of investigators (2 through 8) to determine the rate of interfacial phase change for the case of a suddenly pressurized one-dimensional modelo However, the determination of the temperature gradients by numerical differentiation of the calculated temperature distribution in the case of self-pressurized containers may be 33

difficult for two-dimensional problemso This would occur in cases in which the temperature gradients near the interface, that cause the phase change, are large. These large temperature gradients exist in a very thin layer near the interface, as has been shown by experimental measurements (13914). This requires the use of a very small grid size in order to obtain an acceptable approximation for the temperature gradients near the interface~ Computations have shown that the spatial-grid size has a considerable effect on the calculated rates of heat transfer at the interface. These calculation have been performed for the one-dimensional problem described by Eqs.(l), (2), (3), and the results are shown in Figures 4, 5, 6. Figure 4 gives the interfacial velocity and indicates that the effect of grid size is large at early times. As time progresses the interface velocity becomes less influenced by choice of grid size. This is also seen in Figure 5, which shows that for suitably small grid size the rate of condensation approaches a value independent of grid size. However, the influence on total mass condensed (or evaporated) is a cumulative one, as also noted in Figure 5, and can give rise to serious error in cases where the ratio of interfacial area to vapor volume is large or for small times. Figure 6 shows the influence of grid size on the csomputed temperatures for two different time periods of 10 and 100 seconds. As the grid size decreases the finite-difference calculations of the rate of interfacial phase change may be expected to approach. the exact solution. The example presented above dealt with an externally pressurized system, and the severity of the temperatures gradients at the interface will depend on the magnitudes of the pressure change and the superheat of the pressurant. In the case of self-pressurized containers the liquid is considered initially saturated. Any disturbances from these conditions such as imposed pheat transfer from the ambient will cause the liquid to become superheated. This is an unstable nonequilbrium state, and the liquid will attempt to adjust to a new equilibrium condition by convection and evaporation. The evaporation of the liquid will cause the pressure in the ullage space to rise. The net effect will be to decrease the degree of superheat of the liquid, and increase the mass of the vapor in the ullage volume. If the rate of heating is low9 the conditions in the liquid will be very close to equilibrium conditions. In this case the degree of superheat will be small. This condition may be expected to exist in high altitude flights with a well-insulutted container. Under these circumstances evaporation takes place such that the pressure in the ullage volume rises at a rate which keeps the liquid near an equilibrium state. It thus might be anticipated that the temperature gradients in the liquid and vapor at the liquid-vapor interface wilL be relatively small, and the errors associated with reasonably sized grid spacing will be small. The only true test available is to vary the grid size for given conditions and compare the computed results. In the case of the temperature gradients in the fluid adjacent to the wall two aspects must be considered when either the heat flux or the war,' temperatures are imposed. One is similar to the liquid-vapor interface problem in that the grid spacing should again be small enough to permit sufficiently accurate representation of the temperature distribution in the

LIQUID OXYGEN INITIALLY SATURATED AT 15 PSIA to 7- PRESSURIZED AT TIME = 0 TO O 30 PSIA WITH GASEOUS OXYGEN -_ AT 200 R 6 I 1 i - -xt —-'lAx, v.0167 FT., At ~ 2.5SEC x -2X I1/2 Ax x.00835 FT 5 /t =.625 SEC E \*xi t 1I/2 Ax2 w.004175 FT Jo A \t..15625 SEC. 4- - 3.. —Ax3 Z -- SI 01. 1 1 I I o 20 40 60 80 100 TIME- SECONDS Figure 4. Effect of grid size on computed interfacial velocity.

x- ax.0167 FT At a 2.5 SEC. Ax2 ~ 1/2 A, ~.00835 FT At t.625 SEC. Ax3 a 1/2 Ax 2.004175 FT At,.15625 FT. 0 - 14 u. 12 /7 E 10 1/! 0 Z< / NOTE: CONDITIONS SAME o 2 AS IN FIGURE 4. 0 20 40 60 80 100 TI ME -SECONDS Figure 5. Effect of grid size on computed condensed mass.

-VAPOR 200 TEMPERATURE AT 8 ^ a~ ZERO TIME 195 - o0 at 10 SEC. 0 t90 t 100 SEC A J > O cr185 O VAPOR - - LIOQUID A I 180 / INTERFACE L- |o Axl:.0167 FT. _. At 2 2.5 SEC. F 175 F o ax2 1/2 Ax,.00835 FT At =.625 SEC. Ax3,. 1/2 Ax2.004175FT.0167FT t 100 SEC. At =.15625 SEC 170 10 LIQUID TEMPERATURE AT ZERO TIME t 10 NOTE: CONDITIONS SAME AS IN FIGURE,4. 160 Figure 6. Effect of grid on computed temperature distribution. 37

fluid, either liquid or vapor, and in general dictates that the grid spacing be as small as possible. The degree of being sufficiently small will be dictated in part by the magnitude of the imposed heat flux or the imposed temperature, and the response of the transient convective process to these disturbances under the prevailing effective gravity levelo Again, the only true test available is to vary the grid size for given conditions and compare the coupled resultsO The other aspect to be considered is a physical one involving only the liquid, but is also related to the problem of grid spacing. If the temperature of the solid wall in contact with the liquid exceeds the saturation temperature by some amount, dependent upon various parameters including liquid and solid properties and the configuration, nucleate boiling will be initiated. This particular heating surface superheat might be called the incipient boiling point. If information on the Incipient boiling point is available for the prevailing conditions, an imposed wall temperature below this point then represents no additional problem beyond that of having sufficiently small grid sizes, as discussed above. For the case of a heat flux imposed on the outer surface of the container, however, the resulting wall temperature is a variable dependent upon a number of parameters such as wall thickness and heat capacity, fluid properties, acceleration level and container geometry. Whether the wall temperature will exceed the incipient boiling point will not be known a priori. since the wall temperature is computed during the course of the computations. The procedure by which the possibility of nucleate boiling is taken into account is based on the following physical assumptions: (1) should nucleate boiling begin, further increases in heat flux generally result in relatively small increases in surface temperature as compared to nonboiling convection. This has been observed widely (e.go, Ref. (15)). (2) The vapor bubbles formed are transported by bouyant forces to the ullage volume quite rapidlyo The extent to which this may occur is as yet uncertain, but may be anticipated to depend upon the degree of subcooling present, the pressure and the effective gravity level. These are implemented in the computational procedure by considering that should the tank wall temperature, and hence the liquid adjacent to the wall, exceed the existing saturation temperature by some arbitrary amount, this excess is eliminated by the evaporation of the appropriate amount of liquid directly into the ullage space. In effect, then, a portion of the vapor by passes the liquid-vapor interface. The arbitrary amount referred to above, the symbol for which is given as ATwmax might be considered as the incipient boiling point, and no longer is arbitrary when sufficient information on its behavior is available. They physical phenonmena described above is simulated in the computer program as follows: (1) The container wall temperature is calculated using Eq. (4). (2) The liquid temperature is obtained using Eqo (51). 38

(3) The calculated wall temperature in the liquid region is examined. If it exceeds the saturation temperature-byfmore than the prescribed temperature difference AT max, it then is reduced such that it equals the saturation temperature plus the prescribed temperature difference. (4) Part of the heat added to the liquid region appears as enthalpy in the wall and liquid and the rest is used for evaporating some of the liquid. The portion of the heat transferred to the liquid and resulting in evaporation is determined by setting an energy balance according to fX2Aa qw(x,t)dxdt + 2itrq idr ~a XdT I X 6Tw - f 2dxdr+ 2p taCpr dxr + Wihfg (134) where qil is the rate of heat flow from the interface to the liquid and is given by, qi = k/i a X (135) Equations (134) and (135) are used to determine the rate of evaporation from the interface, wi. If the difference between the wall temperature and the saturation temperature is less than the specified maximum, ATmax, then the procedure above is bypassed and computations proceed as described earlier~ An implicit assumption in the use of this procedure is that the laminar flow conditions described by the momentum and energy equations are not affected. This may be reasonable if the container is relatively large compared to the "bubble boundary layers" region next to the wall. In other words, if the vapor bubbles remain in the vicinity of the wall, the major bulk laminar motion of the liquid will not be influenced by their movement to the ullage space. An accurate physical description of this behavior requires additional analytical and experimental investigation in incipient boiling and the departure and motion of vapor bubbles under low gravity fields with various patterns of subcooling. Do COMPUTER PROGRAM The computer program listing with the assIumption of ideal gas behavior is given in Appendix A and a simplified flow sheet is shown in Appendix B, The MAD (Michigan Algorithmic Decoder) language is used here~ In Appendix C are listed the meanings of the symbols employed, and Appendix D lists the 39

required program inputs. Appendix E presents the results for a single time step during the course of computations of a typical run, and includes the input conditions used. 40

VI. COMPUTATIONS Computations were carried out using the program listing of Appendix A, varying several different parameters so as to demonstrate general trends. As a review, the general assumptions incorporated into this particular program are listed below. A. GENERAL ASSUMPTIONS 1. Cylindrical tank having flat ends~ 2. Acceleration or body forces act along the axis. 3. Two dimensional conditions prevail, with variations only along the axis and radially. 4. Laminar flow conditions within the entire container~ 5. The vapor behaves as an ideal gas. 6. The liquid has constant properties except in the body force terms. 7. The tank side-wall is uniform in thickness and has constant properties. The wall is lumped in the radial direction but axial conduction is taken into consideration. 8. The imposed heat flux on the outside of the tank wall is uniform, but may differ in those portions in contact with the liquid and vapor. 9o The ends of the tank are adiabatic, and the heat capacity of the ends is neglected 10. The grid size varies in the liquid and vapor region as liquid fraction changes, in order that a nodal point always exist at the liquid-vapor interfaceo 11. Initial conditions of uniform temperature and zero velocity exists within the container. 12. The Bond number is sufficiently large that a reasonable approximation to a flat liquid-vapor interface exists~ B. VARIATIONS POSSIBLE IN PROGRAM By relatively minor modifications to the program, some of the above listed general assumptions can be relaxed, providing additional flexibility. In any case, however, axial symmetry must be maintained~ 1o The imposed heat flux can be varied axially and with time. 2. Specified initial conditions of temperature and velocity can be utilized. 3. Specified heat flux to the tank ends, including variations with radius and time, and including heat capacity can be incorporatedo 4. If the description of the physical process warrants, imposed temperatures of the container walls can be /utilized, with variation axially and with time. 41

5. Axial variation of tank side wall thickness and variation of specific heat with temperature can be accounted for. 6. Radial variations in the tank wall temperature can be taken into consideration. This may be particularly dc~'~r-able if the wall is of composite construction. Although the influence of variations in liquid properties with temperature and pressure can be incorporated with minor changes, the use of real gas properties in the vapor space will require major modifications to the program. Major modifications also are necessary to handle the spacewise variation of grid size in either the liquid or vapor domains. C. VARIABLES The variables listed below were maintained constant for the computations presented here: 1. Grid size: 21 Radial x 31 Axial Nodel Points 20 Tank Diameter: 5 feet 3. Tank Height: 10 feet 4. Tank Wall: Aluminum —OoO1 feet thickness 5. Fluid: Liquid and Gaseous Oxygen, initially saturated at 15 psia, with zero velocity. The grid size was made as small as possible for reasonable computational times on the IBM 7090 computer. It might be noted that decreasing the grid size by one-half would result in an increase in computational time by a factor of approximately 16. The relation between real physical system time and computational time depends on a number of factors, bl.t primarily on a/g level and heat flux. This arises from the stability requirement on the computational time intervalthe higher the fluid velocity the smaller the time steps, hence the greater computer time required to cover a given amount of real time. For example, at a/g = 10-5 and q/A = 1 BTU/hr-ft2, 30 minutes of computational time results in 200 minutes of real time, for a ratio of approximately 7 to 1o On the other hand, at a/g = 10-1 and q/A = 1 BTU/hr- ft2 this ratio was 1 to 35 The parameters listed below were varied in order to demonstrate their influence on the pressure rise and the mass of liquid evaporated: 1 AT wmax: 0 to 4~F. This is the maximum permissible wall superheat. The value of zero corresponds to neglecting the the heat capacity of the container walls. 2. q/A: 1 to 72 BTU/hr-ft2 imposed on wall. The heat flux to the liquid and vapor portions of the tank was made the same~ 3. Ullage Fraction: 33% to 67yo 4. a/g: 101 to 10-5. 42

D. RESULTS 1. Influence of ATwmax To demonstrate the influence of the specified maximum wall superheat, ATwmax on the computed pressure, computations were carried on for several values of this parameter at a relatively high and low heat flux. Figure 7 shows the pressure rise and total mass evaporated for the heat flux of A = 72 BTU/hr ft2. A continuing effect of ATwmax up to 4[F is noted The computer output listed finite nonzero values of the quantity designated DMB, indicating that evaporation was by-passing the liquid-vapor inter:ace. It is ant1icipated that a sufficientl.y large value of AT.wrax exists such that a further change in.it would produce no change in the pressure-time behaviior. This resulting wall temperature would probably be considerably above the incipient boiling point, so that the results would have no physical significance. Figure 8 shows the corresponding results with a relatively low heat flux q/A - 1 BTU/hr-ft2. No differences accrue when ATwmax is changed from 1~F to 2~F, indicating that the "natural" wall superheat is less than 1~F for the given ccnditionso This is also demonstrated by the fact that the quantity DMB is zero for ATwmax of 1~F and 2~F. The results for ATwwmax O0 correspond -to the case where the heat capaicty of the wall is neglected,, After 120 minutes, it is noted that the rate of pressure rise is approximately the same for both cases where the wall heatJ capacity is considered and neglected, although the pressure level is somewhat differentr It is also noted in the lower part; of Figure 8 that the t;otal mass evaporated levels off for ATwmax = O. This also occurs with the high heat flux shown in F:i.gure 77, and it is believed will occur with -the low heat flux case for ATL max = ~F, in PFigure 8, at longer periods of time. Even though the rate of evaporation levels off, the rate of pressure rise continueso This is attributed, to th e increasing effect on pressure rise of the heat t5ransfer tjhe ullage space as compared to the heat-, t ransfer tI.-o the liquid. 20 INFLUENCE OF HEAT FLUi3X Figuire 9 shows the influencee of heat flux on pressuroe and lass evaporatred for the case with ATw max= OF,'which. corresponds to neglecting the wall heat; capacity. An additional heat flux of q/A:- 1 BTU/hr-ft2 was -used, but the effects were too small to be shown on Figure 9 because of -the scale. A crossplot of Figure 9 is shown on a log-log scale in Figure 10 for various time periods, with this additional data showing the effect of heat flux. It is noted that the pressure rise is approximately a linea:r~ funct-vion of heat flux. Whet her this would hold for Iother values of ini tial ullage cannoIt be stated at this time~ It might be anticipated that these results would be modified at the early periods of time with values of ATwiax other than zerom

q/A = 72 BTU/HR-FT2 a/g =10-5 33% Initial Ullage max 30 -0 F 40F 25 20 15 0 20 40 60 80 100 120 TIME - MINUTES 4 ATwmax 0 ~F 3 2 5 2 \4 F < 1 0 20 40 60 80 100 120 TIME - MINUTES Figure 7. Total mass evaporated and pressure rise under high heat flux. 44

15. 4p q/A = 1 BTU/HR-FT2 AT a/g = 10-5 wmiax 33% Initial Ullage 0oF V2 15. 3 3 1OF and 2~F 5 15. 2 CO q/A = BTU/HR-FT2 15. 1 15.0 0 0 40 80 120 160 200 TIME - MINUTES 0. 3 m X-~ WTlax ~ 0.2 CO. / 1~F and 2~F 0 0 40 80 120 160 200 TIME - MINUTES Figure 8. Total mass evaporated and pressure rise under low heat flux.

a/g - 10-5 ATvwriax = 0~F 30 33% Initial Ullage 25P |q/A = 72 BTU/HR-FT2 20 q/A 10 BTU/HR-FT2 15 0 20 40 60 80 100 120 TIME - MINUTES;4 4 q/A 72 BTU/HR-FT2 3 O0 > 2 q/A - 10 BTU/HR-FT2 0 0 120 40 60 80 100 120 TIME - MINUTES Figure 9. Effect of heat flux on pressure rise and mass evaporated. 46

a/g = 10-5 Twvmax 0~F 33% Initial Ullage 10 n CC q/A BTU/HR-FT 2 Figure 10. Effect of heat flux on pressure rise. 47

3. INFLUrENCE OF ULLAGE FRACTION F'Kigure 11 shows t5hat the irnfl-uence of initifal ullage vol-ume on pressure rise and mass evaporated for the case where q/A - 1 BTU/hr-ft2, a/g = 105 5amd ATwmax = 1~F. In this case th.e wall superheat never reaches the specified maximum value of ATwmax = 10F, and hence all of the evaporation takes place at; te llqudvapr L.nterface.. In terms of the compfut-er program the quantity DMP is always zero. Figure 12 is a crossplot of Figure 11 for three levels of time. It is noted that the rate of pressure rise in.creases as the initial ullage fraction increaseso This is contrary to the trend observed in an experimental study (9), and may indicate that the level of heat flux is an additional factor governing the effects of ullage volnume. Added to'the upper part of Figure 11 are -the cases for 100% ullage, as computed from Eqo (83)., and 0% ullage consLsten-tr with the assumption of an inec:mpressible liquid, The behavior of the intermediate values of ullage is consistent with these limitso It is further noted from Figure 11 that after an initial start:ing transient of approximately 100 minutes that both the pressure and mass evaporated become essentially linear functions of time. 4/o INFLTJENCE OF ACCELERATION LEVEL The effect of accelerati;on or body force level on the pressure rise and mass evaporated is shown in Figure 13o Again, the wall superheat never reaches the specified maximrum value of4 AT wmax= 1iF. As might be anticipated, increasing the acceleration level:increases bot;h the pressure rise rate and, -the interface evaporation rate, owing to the increasing convecti.ve velocities n=uaced in both th.e liquid and vapor Figure 13 also demonstrates, as discussed earlier, the relation between ccmpu`t~ational t.:ime and real timeo In each of'-he three cases shown the comutoat.ional time was 30 minuteso The increased fluid velocities associated with -the larger body forces result in smaller incremental time steps from the st+abilitys equirements. If the body forces present in a particular physical appl-ication are considerably smaller -than a//g = 10-5 as in deep space probes, it, can be anti.cipated that the behavior for l.ong periods of real time canl be described with reasonably small amounts of total computer time. 50 REPRESENTATIVE TEMPERATURE PROFILES In FiFgure 14 a e are plot~tfed t;he axfal temperature distribut;ion for several levels of time, for five different radial locations from the centerline to the wall. All temperat-ures are expressed as the temperature increase above the initial temperature.

q/A = 1 BTU/HR-FT2 a/g = 10-5 AT.wmax = 10F 15. 3 Initial Ullage;2 167%Jo 10 33 w 15.2 100%o 5 4 ~~Initial; 15,1 Ullage P- 15.1 0% Initial Ullage 15.0 0 40 80 120 160 200 TIME - MINUTES 404 /0.04 0.03 0 a~C 0.02 0 010 40 80 120 160 200 TIME - MINUTES Figure 11. Effect of initial ullage on pressure rise and mass evaporated. 49

q/A = 1 BTU/HR-FT2 a/g = 10-5 ATvmax =1 OF 15. 3 Cn 15.2 DI uD: |.120 Minutes;.4 15.1 80 Minutes 15.0 15.0 L I 1 40 Minutes I 1.0 0.8 0.6 0.4 0.2 0 INITIAL ULLAGE FRACTION 0.04 H 0.03 0 > 0.02 r.2 t~Q 0.01 120 Minutes O iH |40 Minules 80 Minutes 1.0 0.8 0.6 0.4 0.2 0 Figure 12. Effect of initial ullage on pressure rise and mass evaporated.

q/A = 1 BTU/HR-FT2 15.3 33% Initial Ullage AT-wmax 1 OF C 15.2 X 15.1 15.0 I I I I 0 20 40 60 80 100 120 140 160 180 200 TIME - MINUTES 0.03 0 0.02 P 0.01 0 0 20 40 60 80 100 120 140 160 180 200 TIME - MINUTES Figure 13. Effect of acceleration on pressure rise and mass evaporated. 51

0 0.2 0.4 0.6 0.8 1.0 3L,-Container Top 31 29 - -,-..-.. Vapor 27 ~E - Region.C LO 25 - c\ w co 23W1C, J, Liquid Vapor //~ /Interface c21 1 9 9 f j\oes n I r/R =O (Centerline) 7 L33/O Initial Ullage q/A= IBtu/Hr-Ft 2 ~5 a ~a/g=10-5 3 ATwmax = IqF rI.Container Bottom 0 0.1 0.2 0.3 TEMPERATURE RISE (T-To)-~F Figure 14. Axial temperature distribution. 32

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 I I I I I I I I I 3 1 ___,,"_ Container Top 29 (U27 vVapor o G F27 ~' /- Region 25, 0 = l 23'- rf) 05 Liquid-Vapor W Interface m 21 2t- 32.5 Minutes 0 13 7 q/A =IBtu/Hr-Ft2 3 (b) 19 M Container Bottom 0 0.1 0.2 0.3 0.4 0.5 R0.6 FiguCAre ]IDCD Initr ined.5

0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 31 31 Container Top 29 27? C2 3 C L Vapor 25 i 5 Region 23I 23, liquid-Vapor r- 21 m~ ~ 32.SMinutes D19 /'-96 Minutes z z 15 | 17 ~es l e, i~~ Liquid i3o~ - - ~~~Region 3,_d /\Tma =r/R =.5 > 9 21x31 Grid 7 33% Initial Ullage q/A= I Btu/Hr-Ft a/g = 10-5 3 ATwmax = IF (I) I"-Container Bottom 0 0.1 0.2 0.3 0.4 0.5 0.6 TEMPERATURE RISE (T-To)- ~ F Figure 14. (Continued).

0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 ~3 ~~1 L ~~,,Container Top 29 27 so 25 to Vapor Region a z 153...-'~9'~"" /Liquid Vapor X 17 z 15 \e' Liquid Region r/R =.75 30,,' 21 x31Grid m> 9 33% Initial Ullage q/A= I Btu/Hr-Ft2 7 aa/g = 10-5 5 ATwmax = I ~F 3 (d) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 TEMPERATURE RISE (T-To)-~F Figure 14. (Continued). 55

0 I 2 3 4 5 6 I ~~~I III ~~31 2~,.Container Top 2931 ~29- U Vapor 0 Q 23 o = C ~ Region 27 Va p 25 -~ o~ c cn -jz ~ ~ c rOI 23 \9I,\~3 o 7/R-(Sde all 9) Liquid Vapor > 9 33O/1 Initinterface a0~0C~ I~ d I 7 r() ~~~~~q/A~l Btu/Hr-Ft2 a /g=10-5 Container Bottom o 0.2 0.4 0.6 0.8 1.0 1-.2 1.4 1. 6 1.8 2.0 2.2 2.4 TEMPERATURE RISE (T-T0)-OF gre 21 (CnlddI LI..5 m x 7 -o. _. ~ iqi ~~~~~ ~~~~~~ _~~ = Regio TEMPE~~~~r/TR- RI(Sie Wall)k —~~~~Fgr I1_. SCnlue

The relatively rapid increase in temperature of the ullage space as compared to the liquid region is noted. In fact, about one half of the liquid undergoes no change in temperature at all during this period of time. The temperature gradients near the centerline are such that condensation may be taking place there simultaneously with evaporation toward the walls. Of particular interest is the superheated layer of liquid existing just below the liquid-vapor interface at the later time periods in Figures 14 b-d. The influence of the stagnation region in the tank is observed in the lower part of Figure 14e. In Appendix E the computer input and output are reproduced for the case represented in Figures 14 a-e, for computational step number 154, which corresponds to a total elapsed real time of 192 minutes.

APPENDIX A COMPUTER PROGRAM LISTING

$ CC&PILE MAD,PRINT OBJECT,EXECUTE,DUMF,I/C DUMP,PUNCH OBJECT 01S1el 07/C8/67:2. PM M"D (17 MAY 1967 VERSION) PRCGRAM LISTING...DIMENSION U( 980,'OIM),~{(;EC',DIM)',1( 9g80,DIM),SF( 980,01DM), *CC1 1 W( 980,DIM),TO( 980,tDIM),WC( _80,DIM),SFI(31), *CO1 2 EJ'31),F(31),DI(31 ),DR(31 ),R2(3 1 ),R3( 1,CR4 (1),R1(31 ),R *CC1 3 2(31),R3(31),R4(31),R5(31),A1(31},A2(31),A3(31),C2(31 *C01 4 ),C1 80iM),C3( 980,DIM),EJ1(31),C(l31),D2(31),R6(31),CRP *CC1 5 (31),UL(31), DR5(31),DR6(31),CR7(31),CR8(31},R6(31),A4(31) *001 6,VG(31),RC'( 980,DIM), OTCPh(46),CTCXG(31),UO( 9EC,DIMI,VC( S *CC1 7 80,DIM),DTDXL(31),HL(31),HG(31},C4( S8C,CIM) *C01 VECTOR VALUES DIM=2,1,0 *C02 INTEGER I,J,MN,N'TNE,NE.1,NE2, L,NR,NCH,P,J1,J2,I1,14,NP, *C03 i N1,N2, N3, TAU1,NMAX *C03 EXECUTE FTRAP. *CC4 START REAC AND PRINT DATA *CC5 D IM( 2)=N+1 *CC6 EXECUTE- ZERO.(U {I,1- )... U(P+I,+',V"i,l-'.'1) [ P+N+T(l....-*IN+1(11 C07.... 1...T(P+1,N+I), SF(1,1)...SF(P+1,N+1)}, W(1,1)...dP+I1,N+1), *CO7 2 TO (1,1 )....TC P+.,N+I), NT,TIE,NEl,NE2,N2,U,TAU, *CC7 3 UL(0){..UL(N+I),T0XG(0)... ETXG(N+1),CTORW(O)...DTDRW(N+I), *C07 4 RO( 1,1'j..)... RO(P+'1,'N+I),P1, UO(1,1)... UG(P+1,N-+1),VC(1,1 )...VC *07 5 (P+IN+1),D MTDMBTDM,A6,P6,C4(1,1)...C4(M+1,N+1),NPI,MB ) *C07......i i="M/ 2*+1 *........................................ Y=AOVERB *.c.s ON PR=NEW/ALPHA *010 I- PRV=NEWV/ALPHAV *011 RNEW=NEWV/NEW *C 12 RALPHA=ALPIAV/ALPHA *C 13 DR- 1.0/Nq C14 CR2=DR*DR *C15 GAM=CPV/CV....-*C16 DYY=2. O*DR2 *C 17 0R6=DR2 /2.0 *018 PRESS=PC *C19 TSAT=C5 *'IPRESS.P.C6) *C20 TINIT=C5*(PO.P.C6} *C21 CONST = A*A/ALPHA *022 CONSTI=Y*PR*G*BEITA*(A.P.3 ) /NEW/NEW *023 CONST2=1. O/CONST *024 CONST3= 1.0/CONST1 *C25 PRINT RE-iES"LTS COL"TS-"'S T "c'Jt-, C"NSTi, CONST3,P,PRV *C26 M1=AOVERB*QSURFL* (A.P.5)*BEITA*G/ALPFA/ALPhIA/NEW/ROW/CPW *C27 1 /DELTA *027 M2=N<L*A/ALPHA/ROW/CP D/EL TA *028 M3=N'G*A/ALPHA/RCW/CPW/EELTA.*C29 M7=AOVERB*QSURFG*(A.P.5 )*BEITA*G/ALPF-A/ALPFA/NEW/ROW/CPW *C30 1 /OELTA *030 *C30 MS=CONST1 *T IN IT *021 02R=DR/2.0 *032 AOVRB2=AOVERB*AOVERB *033 THROUGH CC,FOR J=2,1,J.G.N+1 *034 AI( J )=DR2*DR2*(J-1.)*(J-1.) *C35 01 C2(J)=1.+1./(2.*(J-1.)) *C36 01 A4(J)=1.0/(iJ-1.0 }D/OR/4.0 *037 O1 CC RI(J)=12*(J-1.0)*CR2 *C38 01

THROUGH GAST9FOR I=M+*1,1,I G.Pl *C39 THROUGH GAST,FOR J=1tl,J.G.IN *40 01 CAST T(IJ)=0.*C41 02 WHENEVER CODE E.. 1 *C42 REAC FORMAT INP'UTAPRESS, PC, X, NT, TIME, U, *043 Cl 1 T(ll)...T(P+1,N+i1),SF(1,1)...SF(P*1,N+l) *C43 THROUGH HUSSFOR I=1,1,I.C.F+1 *C44 01 THROUGH HUSS, FOR J=1,1,J.G. N1 *045 01 01 FUSS TO(IJ)=T(IJ) *C46 01 02 END CF CONDITIONAL * C C47 01 N3=1 *C48 PV=PRESS *C49 BACK NT=NT+1 *C50 HFG=.00745*(-1.1383* (TSAT-18O)*( TSAT-198 )*(TSAT-216)*3.238* *051 1 (TSAT-162)*(TSAT-198)*(TSAT-216)-3.025*( TSAT-162)*(TSAT-180) 2 *( TSAT-216 )*.925*( TSAT-162) * (TSAT-18C) *( TSAT-198)) *051 PP=PRESS/PO *052 I=(P+M)/2+1 *053 J=N/2+1 *C54 BEITAG = 1./(TINIT+T(IJ)/CCNST1) *055 M4=(GAM-1)/GAM *056 CT=1. *057 DX=X/M *058 EX2=DX*DX *05G DZ= (1-X)/(P-M) *060 Z12=0L*0Z *061 ORZ=CR2/0Z2 *062 A3=ORZ*Y*Y *063 8J1=2.*( 1*A3) *064 \XR=0X2/DR2 *C65 r\) CRX=0CR2/0X2 *066 A1=,RX*Y*Y *067 A2=0XR/Y/ Y *68 D2Z=.5 *069 ZOVR2=Z/2.0 *070 BJ=2. *( 1.+A1 *071 80=2.O*PR*Y*Y/CR2+2. C'*PR /2 *0C72 EJ(1)=O *073 EJ1(1)=O *074 THROUGH JE, FOR J = 2,1, J.0. N+1 *075 A2(J)=CX2*DR2*(J-1.)*(J-1.)/Y/Y *0C7 0 A3(J)=DX2/Y/Y/(J-1.)/CR2/2. *077 01 D1( J)=A J-C2( J)*EJ ( J-1 *078 01 R2(J)=12*(J-1.)*DR*DX *079 Cl D2 (J)=8J1-C2 (J )*EJI( J-I) *80 01 EJ1(J)=(1-1/(2.C*(J-1.C)))/02(J) *081 01 R4(J)=Y*Y/(J-1.)/(J-1.)/CR2/CX2/2. *082 01 R5(J)=YsY/(J-1.)/(J-11.)/ /R2/02/2. *083 01 R6(J)=12*(J-1 )*OR*DZ *084 01 JE EJ(J)=(1-1/(22.*( J —1. )) /D1(J) *085 01 WHENEVER COCE.8. 1 *C86 THROUGH FOUR, FOR J=211,J.0.N *0 FOUR SF(M+1,J)=O *088 01 01 THROUGH MPM, FOR 1=2,1,I.G.M *0CE 01 THROUGH MMM, FPR J=2,1,J.G.N *090 C] 01 YMM W(IJ)=((SF(I+1IJ)-2 *SF(IJ)*SF(I-1,J))*Y*Y/DY2-(SF([,J+1)- *091 O1 02 I SF(.1,J-1))/DR2/(J-1.0)/2+(SF(IJ 1)-2.C*SF(IJ)+*SF(IJ-1) )I/ *91 2 DR2)/DR2/(J-1.O) /(J-1.C) *C91 THROUGH QQC, FOR I=2,1,I.G.M *S92 Cl CQ W(IN+ 1J= (.*SF(IN)-5.F(I1N-1)-7.*SF( IN+1))/DYY *093 CL 01

THROUGH. QQQQQ,r FOR I=M+l,1,I.G.P+l *094 0 THROUGH QQQQQ, FOR J=191,pJ.G.N+l'a'- C5 1 0 cc(;Cc RO(I,J)=CCNST1*PRESS*144./RCAS/(T(-I,J)*M5)_____ *C96 01 0 THROUAGH QQQcO 1+,;.G.P *C97 01 Qc~c( W(IPN-.Il)=(8*SF([,fN)-SF(1,N-1)-7*SF(I,N*1))/OYY/RC(I,N+lI *C98 Cl 0 UG=-U*(ROL/RO(m+11) 1,.O) *099 Cl THROUGH FIVE, FCR J=2,11,J.G.N+l *100 Cl FIVE SF (M+1, J)=RO'(M*1,j )*UG*DR2*(J-1. 0) *(J-1.C) /2.0 *10 Cl 0 THROUGH WWW, FOR I=M+1,1,I.G.P *102 __C U(i,N).=2.*(3*SF(1,N,)-6*SF(I,N-1)+SF(I,N-2))/Rl(N)/RC(1,N) *103 01 0 U (I,l) = SF (I,2) /DR 6/R0( I Il1) *104 01 0 WWW U(1I,2)=2.*(6.*SF(1I,3)-SFII,4)-3-.*SF( 1,2) )/Rl(2)/RO(1,2) *105- Cl 0 THROUGH WWWW, FOR I=t'+1,l,1.G.P *106 01 THROUGH WWWW, FCR J=3,1,J.G.N-1.. *107 Cl 0 WWWW U(I,,J)=(SF(I,J-2)-8*SF(IJ-1)48*SF(_I,J*1)-SF(I,J-2))/Rl(J) *10801 2 1 /RO(I,J) *108 THROUGH EFF, FOR J=2,1,.J.G.th *109 Cl V(M+2,J)=2.*(3.*SF(M+2,tJ)*SF(M*4,J)-6.*SF(M+3,J).2.*SF(M+1,J) *110 1 )/R6(J)/RO(M+2,j) *110 VG(J)=2.*(SF(M*3,J)-8.*SF(I'*2,J)*7.*SF(M+1,J))/R6(J)/RC(M*1,J *111 ci 0 EFF v(tP J=.*4(6.'4* FP1, 3." SPFtPJ-SF(F-,J ) )/Ff'(J )/Rt P, J *12C 0 THROUGH EFFF, FOR J=-2,1,J.G.N *113 01 THROUGH EFFF, FOR I=tM+3,1,I.C.P *11-4 Cl 0 EFFF V( I,J)=(SF(1+2,J)-8*SF(1il,J).8*SF(I-1,J)-SF( I-2,J))/Rl(J)/RC *115 01 0 1 (Ili) *115 THROUGH SSS, FOR I=M*2,1,1.G.P *116 0 THROUGH SSSv FOR J=2,1,J.G.. *1701 0 sss W( [,J)=((SF(1+1,J)-2.*SF(T,J)*SF(I-1,J))*Y*Y/0Z2-(SF(I,J+1)- *118 Cl 0 1 SF(1I J-1))/DR2/(J-1. )/2+(SF(I,J*1)'-2*SF( I,J)tSF( I,J-1) )/DR2)/ *118 2 0R2/(J-1.)/(J-1.)/RO(I,J) —(U( I,J)*(RC(I,Je-1)- *1i8 3 RO(I-,J-1) )/02R-V(1,J)*(RC(I-tl,J)-RO(I —1,J) )/C2Z)*A4(J)/RO(I,J *118 4 ) *118 THROUGH PPPP FOR j=2,1,J.G.N *119 0 W(P+1.,J)=(-SF(P —1,J)+8.*SF(P,j) —7.*SF(P*1,J)) *12001 1 1 *R5(J)/RC(P+l9J) *120 Fpp W(1,J)=(-sFI 3,pJ)8.*SF(2,J) )4R4(J) *121 Cl 0 COOE=2 *122 01 TRANSFER TO M~L *123 Cl END OF CON.C-ITICNAL*140 SMAX=0 *125 THROUGH COOOFOR 1=2,1l,I.C-.M41 *126 THROUGH COO,FOR J=1,1,J.G.N *127 0 S=BO+,.ABS.U(1,J)/DX4.ABS.V(1IJ)/OR *128 0 co0 WHENEVER SMAX.L.S, SMAX=S *129 0 WHENEVER PRV'OG'1. I "" - - --- -*130d BOO=4.*PPV/CR2+2.*PRV*Y*Y/0Z2 *131 0 BO=2.*PRV*( Y*Y/0Z2*1./ER2) *132 01 CTHERWI SE *133 01 BOO=4./DR2+2.*Y*Y/0Z2 *134 01 BO 2.*Y *Y / OZ2 +1./ D R2) __ *135 Cl END OF' CONbITIO'NAL *136 Cl THROUGH OCC, FOR I=M*2,1,I.C.P - *137 THROUGH OCC, FOR J=2!y1,j.G. *138 0 S=BO+.ABS.U( I,J)/DZ+.AES.V (I,J)/[CR *139 0 ECCC WHENEVER SM4AX.1.. 5,SM= — * 140 THROUGH MAXS,FOR I~=M+2,1l,I.G.P *141 ___ S=BOO+.ABS.U( 1,1)/DZ *142 - 0 WHENEVER (SMAX*OT).G. C.8,DT=0.8/SMAX *143 C MAXS WHENEVER sMAx.Lo S,-lSMAX=S~ *144 0

TIME=TIME+OT *145 TAU=TII4E *OS X=X*U*DT __ _ ____*146 _ OX 1=DT/DX *148 071=01/07. * 149 OZ3=CT*DT*Y*Y*RALPH-A/DZ2*15 OZ5=2.*DZ3 *151 07? = 0? *Y Y FR * RNEh 072 *5 C=T *Y *Y*A L PHA W /AL PH A / 0 2 ___*153 GY3=CT*DT*RALPI-A/0R2. —----— __DY1=2.*0Y3 *155 0Y7 = DT * RNEW PR / 02 *5 DY5=4.*DY3 *156 OX 3=DT*Y*Y/DX2 -*158 0R3=N*N*DT *159 CR4=N*DT....*160 0DR5=4.*0R3 *161 OX5=2.*DX3 *6 DR1=2*DR3 *163 OX7=DX3*PR *164 CX8=DT*Y*Y*ALPHAW/ AL PHA/CX2.. ~*165 DR7=0R3*PR *166 C1=1.O-0X5-0R5 *167 C2=1.0C-X5-CORI *168 C3=1.0-2. C*CX7-2.O*0R7 *169 THROUGH CD,FOR J=211, J.G.K*1 *17C FI=0.5/(J-1.O) __*171__0 DI=1.51( J-1.0) *172 0 DR1( J)=0R3*( 1.0-Fl *173 0 DR2(J)=CR3*(1.O0-FI) *174 0 DR3(J)=0R3*(1.0-CI) *175 0 DR4(J)=0R3*( 1.0*01) *176 0 R3( J)=PR*0P3/(J-1. C) /2.0 *177 0 DR5fJ)=CT*RALPHA*0RL(J) *178 0 R 6 (J )= C T*R A LP hA *0R 2 (J) *179 0 DR?7U) = CR3(J) * PR * RrNEW * 180 0 cc DP8(J) = CR4(J) * PR * RNEW *181 0 THROUGH BBB,FOR 1=24,1I.C.t *182 THROUGH 888,FCR J=l1,1J.G. N *183 0 Cl( I,J)=U(I1,J)*DX1 *184 0 888 C3C1,J)=V(1,J)*DR4 *185 0 T(1,N*1)=TC(1,N*1)+DT*tMl-CT*pJ2*(11*T(1,N*1)-l8*T(1,,N)*9* *186 1 T( 1,N-13 —2*T(1,N-2))/6/CP*2*CX8*(TO(2,IA*1)-TO( 1,N*1)) *186 THROUGH CBAFCP 1=2,1,I.G.M+1 *187 CEA T 1,Ni-1)=T(1,N+1)4-CT*M1-OT*tM2*(11*T(1,N*1)-18*T( 1,Nd*9*T(l,t\-..*188 0 1 1)-2*T(1,N-2))/6/OP*DXe*(TO(1*1,N*1)-2*TO,( 1,N+l)*TC(I-1,N*1)) *188 THROUGH BAAA.,FOR 1=P,*2,1l,1.G.P * 189 EAAA T(1,N-I1)=TOC I,N*1)*OT*M7-O1*t'3*(11*T(I,N*1)-18*T(1,N)*q*T(l, *190 0 1 N-1)-2.*T(I,N-2) )/6/CR*0Z8*(TC( I*1,N-+1)-2.*TO(1,N*1)*TO(I-1,N *190 2 *1)) *190 T(P+1,N*1)=TO(P*1,N*1)*OT*M7-ET*M3*(.11*T(P*1,N*1)-.18*T(P.1,N) *191 1 *9*T(P+1,N-1)-2*T( P*1,N2 )//R2C8(TO(P,N t1)-TOtP+1l,N+l) *191 2) *191 THROUGH BA,FOR 1=1,1,I.G.F*1 *192 THROUGH BA,FOR J=l,1,J.G. Nil *193 C RO(I,J)=CCNST1*PRESS'*144.0/RC-AS /(T(1,J)+M5) *194 0 WO(1I,JA)=W (IJ) - -----— ~*195 -0 EA T0(IJ)=T(IJ) *196 0 THROUGH GRAD, FOR J= 1,1, J.G. N *197 0TDXL(J)=(11*T(M1,PJ )-18*T(NtJ)-.9*T(MJ-l,J)-2*-T(M,-2,J))/6/OX - *198 -0

DTOXG(J)=(-11*T(M*1Jj)*18*T(td*2,J)-9*T(M*+3,J)+2*T(M,+4,PJ))/CZ/ *199 0 1 6 *199 GRAD UL(J)=(CTCXL( J)-OTDXG(J)*KG/KL)*CP*NEW*NEW*ACVE~e/PR/ *200 0 1 B EI TA /G /H F GA.P. *200 U0O *201 THROUGH AVER,FOR J 1,11, J.G. N *202 AVER U=U*+(UL(J*1)+UL(J) )*(2.C*J-1.C)*DR2/2.C *203 0 THROUGH ACAA,FOR 1=1,1,I.C.P*1 *204 ACAA DTDRWf(I)=(11*T(IpN*1)-18* 1(d**( N1-2TI~-)//R*2C5 0 Qw=0 *206 THRCUGH CAAA, FOR I=iv*1,1,I.CG.P _*207 CAAA OW=QW+(OTDRW(I)*DTDRW(1*1))*CZ/2. *208 0 01=0 *209 DQ1L = 0 *210 THRCUGH DIAM, FCR J=1,I,J.C-.t *211 001L = OQI.L+(DTO-XL(J)*DTODXL(J41) )*(2.0*J-1.C)*0)R2/2.0 *212 0 CAAA QI=QlI(TCTXC(J )*OTDXC-(J+1) )*(2.C*J-1.C)*0R2/2.C *213 0 OQWDT=QWif*KG*6.28*NEW*ALPHA/EEITA/G/A/A/AOVERB/AOVERB *214 RQIN=DQWCT*AOVER8/6.28/A/A/CSURFG./(l1-X) *215 DQ1LOT.= 3.14*OQIL*KL*NEW*ALPFA/BEIIA/O/A/A *216 DQIDT=3.14*C1*KG*N\EW*ALPHA/EEI'hAlG/A/A *217 OME = DT*U*ROL*(A.P.3) * 3.14/AOVERE *218 DMT=DMT*OME *219 EM=-EMB*OME *220 PRESS=PRESS*RGAS*ACVERB/(3.lA*144.*A.P.3.* (1.-M)*((COWCT- *221 1 DOIDT)*DT*A.P. 2.I( CV*ALPIHA )-COM*GAMw*TSAT) *221 A6=( PRESS-P1l)/PRESS *222 P6=( PRESS-Pi) /P1 *223 P1=PRESS * 22,4 TSAT=C5*(PRESS.P.C6) *225 A5=2.*DT*AOVERB*G*(A.P.3. )/ALPHA/ALPH-A *226 THROUGH INTERFFOR J=1,I1,J.G.N *227 INTERF T1M*1vJ)=(TSAT-TINIT.)*CONST1 *228 0 T(11,1)=TC(1,1)*CI*TC(1,2)*EXS+TC(2,1)*CR5 *229 THROUGH OO,FCR 1=2,1,1.G.M *230 6HENEVER U(I,1).G. 0 *231 0 C=(TO(1-1,1)-TO(I,Mf*Cl(1,1) *232 Cl 0 OTHERWISE *233 Cl 0 C=(TCO(1,1)-TO(I*1,1)*C1(1,1) *234 01 0 END CF CCNOITIONAL *235 01 0 EC T(1,1)=TO(1,1)*CL*TOU1,2)*ORS*+(TC(1+1,1)4TC(I-1,1))*0X3- +C *236 0 THROUGH EEFOR J=2,1l,J.G.N *237 EE T(1,J)=TO( I,J)*C2*TC(2,J)*CXS*+TC(1,J-1)*CR1(J) *238 0 1 *TO(1,J+l)*0R2(J) *2 3 8 THROUGH FFFOR J=2,1,J.G.N *239 THROUGH FF,FOR 1=2,1,1.G.P *240 0 WHENEVER U(,J) C.. C *241 0 C=(TC(I-1,J)-TO(1,J))*CI(I,J) *242 01 0 OTHERW ISE *243 01 0 C= (TO IIvJ )-TO ( +1*1,J) )*C 1 (I,tJ *244 Cl 0 END OF CONDITIONAL. *245 01I 0'WHENEVER V(I,J).0. 0 *~24 6 0 D=(TO(1,J-1)-TC(I,J))*C3(1,yJ) *247 01 0 CTHERW ISE *248 01 0 D=(TO(I,J)-TO(1,,J+1)*C3(IJ) *249 el 0 END OF CONDIT ICNAL *250 Cl 0 FF T( I,J)= TO(I,J)*C2*(TO( 1*1,J)*TC(I-l,J) )*0X3* *251 0 1 TO'(1,J-1)*DR1(J)*TO( 1,J41I)*CR2(J)1-C*O C THROUGH GGG,FOR 1=21,1I.G.M, *252 THROUGH GGG,FOR J=2,1,J.G.N *253 0

FI=R3(J)*(T(I,J*1)-T(IJ-1)) *254 02 WHENEVER U(IJ).G. 0 *255 02 C=(WO(I-1,J)-Wc(IJ))*CI(IJ) *256 Cl 02 OTHERWISE *257 01 02 C=( WOI I,J)-WC( I +1, J) )*C1 (I,J) *258 01 02 END OF CONDITIONAL *259 0 1 WHENEVER V(IJ).G. 0 *260 02 D=(WO( I J- 1)-WO(I,J))*C3(I,J) *261 CL 02 OTHERWISE *262 - C 02 O=(WC(I,J)-WO(IJ+.l))*C3(I,J) *263 01 02 END CF CONDITIONAL *264 01 02 COG W(IJ)=WO(IJ)*C3*(WC(I+lJ)*kC(I-1,J))*CX7.C+D* *265 02 1 WO(IiJ-1) *DR 3(JJI+WC (I,1J R4 (JFI *265 THROUGH MMMM, FCR J=2,1,J.G.N+L *266 MMMN SF(M+1,J)=O *267 01 NEI=O *268 JJ NE1=NE1+1 *269 WHENEVER NEL.G. NE T TRANSFER TO PRINT *270 THROUGH LL,FCR I=1,1,.G.M*1 *271 THROUGH LL,FOR J=i,1vJ.G.N+i *272 01 IL. C3 (1,J)=SF(I,J) *273 02 1=1 *274 II 1=1*1 *275 F(1)=0 *276 THROUGH L1,FOR J=2,1,J.G.N *277 OI(.J)=(C3(I+1, J)+C3(I-iJ)* Al-W(IJ)*AI(J) *278 01 LP F(J)=(CI(J)*F(J-1)*C2(Jfl/C1(J) *28790 THROUGH LNFOR J=1,1,J.G.N-1 *280 H=N+1-J *281 01 IN SF1 IH)=EJ(H)*SF(I,H+1)*F(H) *282 01 O\ WHENEVER I.L.MTRANSFER TO II *283 RMAX=O *284 THROUGH AAAFCR I=2,NRI.G.M *285 THROUGH AAAFOR J=2,NRJ.G.N *286 01 R3=.ABS. ((SF(IJ)-C3(1,J) )/(C3(IJ)+1.CE-20)) *287 02 AAA WHENEVER RMAX L1. R3,RMAX=R3 *288 02 WHENEVER RMAX.0. EPSICN, TRANSFER TO JJ *289 THROUGH GFH, FOR 1=1,1,I.G.M *290 WHENEVER T(I,N*1).G. (T(tP*1,N)*DTW*CCNST1),T(IN*1)=T(M.,N)*TW *291 01 1 *CONST1 *291 GFH WHENEVER T(IN).G.GT( 1,N)4CTS*CONST1),T(IN)=T(M1,N)+CTS*CCNST1 *292 01 DQS = 0 *293 001 = 0 *294 DMB=O *295 CT2=6. 28*A*A*CX*RC*CPW*CELTA/ACVER/2.C/CCNSTI *296 THROUGH SQC5 FOR I=1, 1, I G.M *297 C=TU,N+1)+T( I +1,N*) *298 01 0QS=DQ54(C-C4( I,N+1) )*CT2 *2990 SQO C4(INN1)=C *300 01 THROUGH GHF,FCR J=1,1,J.G.N+l *301 GHF CRI(J)=3.14*ROI*CP*(A.P.3)*C *CR2*(2.0*J-1.0)/ AOVER8 *302 01 THROUGH FHG, FOR I=1,1,i.G.M - 303 THROUGH FHG, FOR J=l,1,J.G.N *304 01 C = (TiIJ)*T(IJ*1)+T(I I1,J)*T J++1, *305 02 DQL = DCL +(C-C4(I,Jfl*0R1(J) *306 02 FIC C4(IJ) = C *307 02 CQI=(6.28*A*A*X*QSURFL/ACVERE*CQILCT)*CT*CONST _*308 DM8 = ICCI-OCL-DQS)/EFG *309 WHENEVER 0MB.1. 0.0, 0MB=C.O *310 DMBT=DMBT+DMB. *311

DMTT=OMT-DMBT *312 OWDT=-.581*PRESS*ALPHA*A*U/ACVER8 *314 LG=-U*(ROL/RO(M*1,1)-.0) *31 5 Cl= 1-DZ5-DY5 *316 C2= 1-Dt5-Dt~l *317 C3=1-2*OZ7-2*DY7 *3218 T HROCU GH BBBbR!=M*-2,.l,I. C.P *1,9 THROUGH BBBB,FCR J=l,1,J.G.N *320 0 ClIIJ)=DZ1*U(ItJ) *321 0 EBBB C3(IJ)=V(I,J)*DR4 *322 0 T(P+l,1)i=TOIP1 0+!,) ~*Ci*+TC(Pt *CZ5+-TC(P 41,2)*DYS *323 1 4-t'4*M6*(M54-TO(P*1,1)) *323 THROUGH BIB,' FUR J=2,1,J.G.N *324 818 T(P*1,J)=TO(P+1,J)*C2*TC(P,J)*DZS*+TO(P*1,J-1)*CR5(J)+TC(P+1,J *325 —0 I1 *1*Ri)N4M*.CP+ I j).) *325 THROUGH TWO, FOR I = I'*2vl,I.G.P *32,6 WtEN-EV-ER —-U1 1)i'-.G. 0 *3 27 0 C=(TO( 1-1,1)-TO(l, ) *C1( 1,1) *32801 1 OTHERWISE *322901 0 C=(TO(Ivl)-TO(1+1,1)*Cl(I,ll *23301 0 END OF CONDITIONAL -— *3 l 0 TWO T(1,1)=TO(I,1)*Cl+TO(I.,2)*0V5*(TO(I*1,1)*TO(1-1,91))*CZ3 *C *332 0 1 *M4*M6*(M5*TO([1,1)) *3232 THROUGH ONE,FCR I=M*2,11,.C.P * 333 THROUGH ONE,FOR J=2,1,J.G.N *334 0 WHENEVER U(IvJ).G. C *335 0 C=(TO( I -1tJ )-T O( I,J) )*CI( I,J ) *33,6Cl 0 OTHERWISE *33701 2 C= (TO (I,J )-TO ( I+I,J))*CI( I,J ) *323801 2 END OF CONDITIONAL *33901 2 WHENEVER V(,J).G. C *340 0 D=(TM(,J-1) —T0(I,J))*C3(I,J) *34101 2 OT HERW I S'E'- -.l.l-__. —.._ * 34201 2 D= (TO( IvJ )-TO( IpJ~lfl)*C3 (ItJ) *34301 2 END OF'CONDITIONAL *344Cl 0 CNE IT I,J)= TO(I,J)*C2*(TO(1*1,J)4TC(I-1,J))*0Z3* *345 0 1 TO(IvJ-1)*DR5(J)+TCCI,J*1)*ER6(J)*C*O *M4*t'6*(IM5-TC(I,J)) *345 THROUGH THREEFCR I=tM+2,1,I.C.P *346 THROuGHT'HREE,9-FCR J=2,1 1,j.'G.N- *'347 0 FI=A4(J)*(-(2*(U(I,J)-UC(I,J))*C1(I,J)*(U(I*1,j)-U(I-I,Jfl4 *348 0 1 C3(IJ)*'(U(IJ.-1)-U(I,J-1)))*(RO(1,J+1)-PO(I,J-1))*tACVRB2*(2* *348 2 (V( IJ)-VC,(1,J) )+C1( I,J)*(V(I*1,J)-VCI-1,jfl*C3(1I,J)*(V(lI,j41 *348 3 )-V(1IJ-1)) )*(RO(1*1,J)-RO(I-1,J) ))/RC(IJ)-W(I,J)*(2*(T(1,J) *348 4 -TO(I,J) )*Cl(I,J)*(T(I+1,J)-T(I-l,J))*C3(IJ)*(T(I,J*1)-T( IJ *348 5-1) ))/ZCVR2/( T(tI,'J)*P5) +A6*W(I J) * WHENEVER U(I,J).G. 0*490 C= (WO II- 1,J )-W ( I, J) )*C1( I,J) *5 l 0 OTHERWISE *35101 2 C=(WO(I,J)~-wO(I*1,J))*CI(IJ) *5 1 0 END OF CONDITIONAL __ *35301 2 W HE NEV ER V (I J )."G. C *354 0 D=ktWC ( I,J-1I)-WO I tJ) ) *C 3 It J *35501 2 OIH ER W I SE *5 1 0 E= (WO ( I,J)-WO( I,J+1 ) ) *C 3 I, - *35701 2 END OF CONDITIONAL *315801 0 ThREE W(IJ)=WC(I,J)*C3+(WO(I*1,J)*WC(I-1,J))*CZ'7_*C+D_+ *359 0 I WO (I iJ-.1)*DR7(J) *WCII + 1) *CR8 IJ) *FI *359 2 eA5*A4(J)*(T(1,J+1)-T( I,J-1))/(T(IJ)*P5)__*5 THROUGH-NNNv F-OR -J2,JG *360

SF(M+l1,J)=RO(M*1,J)*UG*0R2*(J-1.0)*(J-1.C)/2. *361 01 NNN SFI(J)=SF(M*1,J) *362 THROUGH DDAA, FOR I=M42,19,I.G.P *363 CDAA SF(IN+1)=SF(M+1,N+1) *364 01 NE 2=0 _365 JJJ NE2=NE2+1 *366 WHENEVER NE2.G.NETRANSFER 10 PRINT *367 THROUGH LNI, FOR I = M1411, I.G. P +1 *368 THROUGH INI, FOR J=1,1,J.. N1_ _ *369 01 LNL C3(IJ)=SF(IJ) *370 02 I=M+1 *371 IIi 1=1*1 ~.. 8*372 THROUGH 1114 FOR J = 2,1, J G.. N - *373 OI(J)=(C3(1*1,J)+C3(I-1,J))*A3-(W(IJ)*RC( IJ)*+A4 (J)*(U(IJ)* *374 01 1 (RO(1,J+1)-RO(IJ-1))/D2R-V41,J)*(RO(*1,1J)-RC(i-1,J))/02Z))*, *374 2 A1(J) *374 LLM F(J)=(OI(J )F(J- 1)*C2(J) )/L2(J) *3765 THROUGH LIN, FOR J = N,-I, J.1. 2 LLN SF(IJ)=EJ1(J)*SF(IJ)+1)+F(J) *37 01 WHENEVER I.1. P, TRANSFER TO III *378 WHENEVER NT.E.1.ANC. NE2 E. It TRANSFER TO JJJ *379 RMAXX= 0380 THROUGH AAAA, FCR 1=M429NRtI.0. P *381 THROUGH AAAA, FOR J = 2,NRj J G. N *382 01 R3=.ABS.(SF(IJ)-C3(IJ))/(C3(IJ)1.E-2C)) *383 02 AAAA WHENEVER RMAXX.1. R3, RNAXX= R3 *384 02 WHENEVER RMAXX.G. EPSLON,_TRANSFER TO JJJ *385 THROUGH BEEFOR I=2,1,1.G.M 1*386 BEE W(IN+l)=(8.*SF(I,)-SF(9IN-1)-7.*SF(IN*1))/PYY *387 01 \ THROUGH UJUUU, FOR 1= M+2,1,I.G.P *388 CIUUU W( IN*1)=(8*SF(IN )-SFtIN-1)-7*SF( I,N4l))/ODYY./RU1I,N1) *389 01 THROUGH CEE FOR J=2,1,J C-. N *390 W(P+1,J)=(-SF(P- 1,J)*8.*SF(PJ)-7.*SF(P+1, J)-)*R5(J)/RC(P*+1,J _ *391_ 01 CEE W(1( J)= (-SF(3, J )+8. O*SFF( 2, J) )4*R4( J) *392 01 FL THROUGH UUU, FCR I=M+l, 1, I.0.P*1 *393 THROUGH UUU, FOR J=1,1,J.G.N*1 *394 01 VOC( IJ)=V( I-,J) *395 02 LUU UO(IvJ)=U(IJ) *396 02 THROUGH BE, FOR 1=2,1,1.0G.M *397 MA S, N -6. *SF,N-11SFtN-22/Rl(N) *398 01 t(I,1)=SF (I 2) /DR6 *399 01 BE UIl,2)=2 O*(400*SF(1,3)-3.C*F(I,2) -FC / - *0 THROUGH VVV, FOR I=10*1,1,1.G.P * 401 U(ItN)=2.*(3*SF(I,'vN)-6*SF(IN1)+SF(1N-2))/Rl(K)/RO(IN) *402 01 U(I,1)=SF(I,2)/OR6IRO(I___ *403 01 %ivv U(II2)=2.*(6.*SF I,3)-SF(1,')-3.*SFII,2))/R1f2)/RO(I,2) *404 Oi THROUGH B'BETA, FOR 1=29,11.0.M4 *405 THROUGH BBETA,FCR J=391,J.G.N-1 *406 01 BEETA U(1,J)=(SF(1,J-2)-8*SF(IIJ-1)*8*SF(IJ*1)-SFIIJ2))/R1(J) *407 02 U( +l1N)=U(M+1,N-1) *40B THROUGH VVVV, FOR I=M1,1,I.G.P *'s09 THROUGH VVVV, FOR J=3,1J.G.N-1 *'410 01 ~YVV U(1,J)=(SF(IJ-2)-8*SF(1,J-1)*8*SF(IJ+1)-SF(IJ-2))/R1(J) *411 02 1 /RO(IJ) *411 THROUGH BFF, FOR J=2,1,J.G.N *412 V(M+2,J)=2.*(3.*SF(M,2, J)+SF(M+4,J)-6.*SF(M*3,J)4t2*SF(f*1,rJ)) *413 01 1 /R6(J)/RO(M*2,J) *413 V(P,J)=2.*(6.*SF(P-1,J)-3.*SF(PJ)SF(P-2,J))/6(J)/R(PJ)1 01 8FF VG(J )=2.*(SF(M+3,J)-8.*SF(1M*2,J)4-7.*SF(*1, J))/R6(J)/RC(M+e1,J _ *415 01 1) *415

THROUGH BFFF, FOR I=M+3,1,I.C.P-1 *416 THROUGH BFFF, FOR J=2,1,J.G.N - *417 01 BFFF VII,J)=(SF( + 2,J)-8*SF(l+ 1,J)+ 8*SF(I-1J)-SF(I-2,J))/R6(J)/RC.*418 02 1 (IJ). *418 THROUGH FBF, FCR J=2,1J.G.{*+1'419 FBF SFI,M+J)=O 20 O!A0 THROUGH BF,F]R J=2,1,J.G.N *421 V(2, J)=2.- -*(SF(4,J)-6.0o*SF ( 3,J ) +3.0* OSF(2,J))/R2 ( *422 01 V(M,J)=2.0*(6.O*SF(M-1,J)-3.C*SF(M,J)-SF(M-2,J))/R2(J) *423 01 BF V(M+1, J)=2.0*(8*SF (MJ )-SF(M- 1,J) )/R2( J) *424 01 THROUGH BETAFOR I=3,1,I.G.{-l *425 THROUGH BETA tFOR J=2,1,J.G...426 01 EETA V( I,J)=(SF( I+2,J )-8*SF ( I+ 14lJ )+8*SF( I-1 J )-SF( I-2,J ) )/R2(J) *427 02 WHENEVER'NT.L.N1; TRANSFER TO PR I NT 428 WHENEVER NT E.NMAX. *429 PUNCH FORMAT ODATAA,PRESS, PC, X, NT, TIME, U, *430 01 1 T(1,1)...T(P+1,N+I1),SF(,1)...SF(P+lN+1) *430 T R ANSF ER TO PRI NT 43 01 END OF CONDITIONAL *432 01 NP=NPi'433 NP1=TAU/(TAU1*N3) *434 WHENEVER NP1i'-NP-1. G.O,'''..'TRA'~S'FER "TO PRIN T *435 TRANSFER TO BACK. *436 PRINT WHENEVER NT;L. NCTRANSFER 1C NEXT *437 R1=0 *438 THROUGH YZFOR I=2,1, I.G.M.439 THROUGH YZ,FOR J=2,1,J.GNN+ I *440 01 R2=' BAS.-(i(T(I"J)- TO( I, IJ )'/TC ( I, )) * 441 02 YZ WHENEVER R1.1. R2,RI=R2 *442 02 ON I'EXT NUG=O......'-.'*443 NUL=O *444 RAL = 0 *445 RAG = 0 *446 THROUGH CLH, FCR I=l=],[,,G.'*+1 *447 HL(I)=.ABS.(KL*DTDRWI )/A/(T(I,N+ 1)-T(I I,)) *448 01 RAL = RAL+DX*T(I,N+1) -*449 01 CLF NUL = NUL+D'X*HLI)*A*X/ACVERB/KL *450 01 RAL = RAL* ( X..P.3) / iAOVERE.P.4)...*451 THROUGH GLH, FOR I=M+ltl,,I.C.P+1 *452 J="M................................. *453Y 01 HG(J)=.ABS.(KG*DTDRW(I)/A/(T(ItN+ 1)-T(I, ) )) *454 01 RAG = RAG+DZ*T(IN+1).455 01 CLI NUG = NUG+DZ*HG(J)*A*(I-X)//CVERB/KG *456 01 RAG = RAG * ({'-X).P.3) * BEITAG/RALPHA/RNEW/BEITA/ *457 1 (AOVERB.P.4) *457 CL = N'UL/(RA-,P.25)........ *458 CG = NUG/(RAG.P.O.25) *459 UR=U(II,1)*RI(l) /2.0 *460 THROUGH DAL, FCR J=2,1,J.G.N+1 *461 CAL UR=(UR+*U( I 1 J)*R-*iJ) )/2. 462 01 PRINT RESULTS NTTIME,TAUCTURU,PRESSZ,XDX, *463 1 DMEdMTDB', DMBTbCMMMT TT, *463 2 RQIN,,DQIOi,OWOT, QWDT-TFG, U(1,1)...U(P+I1Nl+)V(1,1)...V( *463 3 P+1,N+1) T(1,1)". T(P+ N+ 1),f( )..SF' (+iN+i),VG(1)... *463 4 VG(N+*I),UL(1)...UL(N+1),HL(1)...HL(M+I),HG(1)...HG(P+I-M), *463 5 NUL,NUG,RMAX,R'I,NEE2;'"'-RH - L, W R L G, C''-CL-,GR- PR 1, *463 6 SFI(1)...SFI (N+1(N+ ) CQILCT *463 WHENEVER NE. NE, TRANSFER TO ENC *464 WHENEVER NT.E. NMAX, TRANSFER TO END *465 WHENEV ER' N-2.E'1-*" *466

N3=2.*N3 *467 01 N2=0 *468 01 OTHERWISE _469 01 N2=1.*470 01 END OF CONDITIONAL- *471 Cl TRANSFER TO BACK *472 END CONTINUE *473'TRANSFER TO START *474 VECTOR VALUES INPUT=$3F12.4,F12.8,I 4,Fl2.8,E16.8/(5E16.8)*$ *475 VECTOR VALUES OOATA=$3F12.A,F12 2,I 4,F12.8,E16.E/(5El6.8)*$ *476 END OF PROGRAM *477 THlE FCLLCWING NAPES' HAVE OCCURRED.CNLY CNCE IN IHIS PROCRAM. CCMPILATION WILL CONTINUE. CPV *016 Ni *428 NC *437 TAU1 *434 -1

APPENDIX B FLOW CHART 71

T READ FORMAT INPUT 3 PRINTRESULTS FI t J...UI N APRSS T(IJ. SF(I.J) - ETC. V = IMI I II-. READ AND CONST,CONST 1 O = STRO(TI,J) ="-'- 1PRINT-DATA 3NS 2 CONST 3, I., PRV, _ TDX(J = LJ — _UU= —_-. DTCR =nt *- C-V -. =.. _ QL= -. _ l=..V._ DWT=*-_ RIN =- T0 T L. Wall L. CODE= Ss ~:~~~.+V.. Wal SMA ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~~~SA= Sq S~~~-(D=SA SA<_ BO='DTDxL(J)....)TI (+1Nl DTX1 J H ~~ -~ TO(I'J)=T(I'J) H- WO(I,J)=W(I,J) H RO(I,J):':'H~- T(P+I,N+I)....'~t (I'N+I)='{'H~T(;i'N+I)-'-'' T(i'N+I)....'~~ IME=TIME+DT " Entire Region Entire Region Entire Region'Top Surface V. L.'DTDXG(J).... UL(J)....''U'=U'=+'.' "'-DTDRW=~... aW....'' OQ'L~~L...' 0=":'- DawDT'-'".'~ R'i~.... L.+V. V. Wall II J'1) 4 ol M_ = NE1=NEI~ T+1,)=)T(M+1,N)+DTW*tONST 1 TLO RM X<R R3 NE1>F(,J -t-DTW, H..H -I,.- HT(IN)T(4AN) DS=O B< ~ l+. Df~. ~ F ON*CONS / F s0iM+1,'J)=f(RO(M+f'j)) HWIJ=~H ilj=~HTI1=~ T(P +i,J)= ~ HT(P- Ll) _L HUU~ D__~ H UU DT=~ DMBT'=... SF(P iM,, _,, NE2= H ~ I- A —.. T=H H _ AN d_ NE2!I PLO U(IN)=~~~,, (1,4 =~~~ I 1U(I,1) ~ ~'~-1 U(IN E ~~ I I UO 4 7=U (I, I- V(J)(,) WJz P,) wo -+)=~ V. L.'- L V. V.. V.+L " @~~~~~~~~~~~~~~~~w ) - T=m+f=N I FH) I V. V. L. r- t, V.int. V. PUNCH FORM.T LV(M - (M_41j)=b A APRgSSPO' NT=NMAXY) L BACK ~~~~~~~~~X,NT,TIMEUU TI J,J),SF(a J) Ilt. T<NC~~ ~r R~~~~~~~~~~~~~RI<R2 PRINT' UR = ~~ H CG= ~~H CL= ~~ ~ -RAG= ~~ ~ H UG= ~~ H RAG ~~~t~ HG(J)= ~ L= =..._H-9-G HL= = |PRINT RESULTS'' /H TU H.?:0 < T(I,J),SF(IJ) >NE T1 7=0 /u(l,J~,J) I ~' I X 72

APPENDIX C PROGRAM NOTATION AND NOMENCLATURE 73

Symbols preceeding the expression indicate the following: * Input data ** - Printout Symbols following the expression indicate the following: [0] - Integer [1] - Dimensionless or units (150) - eog., - This expression appears in statement number 150 in the program of Appendix A. A *A - [ft] - Tank radius *AOVERB [1] - A/B AOVERB2 -[1] (A/B)2 Al - [1] - (67) A2 - [1] (68) A3 - [1 - (63) A5 - [1] - (226,359) A6 - [1] - (222,348) A1(J) - [1] - (35) A2(J) - [1] - (76) A3(J) - [1] - (77) A4(J) - [1] - (37) ALPHA - [ft2/sec] - Thermal diffusivity of liquid ALPJAW - [ft2/sec] - Thermal diffusivity of wall ALPHAV - [ft2/sec] - Thermal diffusivity of vapor B B - [ft] - Total tank height BJ - [1] - (71) BJ1 - [1] - (64) BO [1] - (72,132,135) BOO [1] - (131,134) BEITA - [oF-l] - Volumetric coeff. of expansion-liquid BEITAG - [OF-l] - (55) - Volumetric coeff. of expansion-vapor C C - [1] - Temporary storage C1 - [1] - (167,316) C2 - [1] - (317) c3- [1] -(318)

C *C5 - [~F] - (20,225) - Coefficient for computing TSAT as function of pressure *C6 - [1] - (20,225) - Exponent for computing TSAT as function of pressure C1(I,J) - [1] - (184,242,321) - Temporary storage C2(J) - [1] - (36) C3(I,J) - [1] - Temporary storage 4(I,J) - [1] - Temporary storage **CL - [1] - (458) - Coefficient of Rayleigh number correlation-liquid **CG - [1] - (459) - Coefficient of Rayleigh number correlation-vapor *CP - [BTU/lbm-~F] - specific heat-liquid *CPV - [BTU/lbm-~F] - specific heat-vapor-constant pressure *CV - [BTU/lbm-OF] - specific heat-vapor-constant volume *CPW - [BTU/lbm-~F] - specific heat-wall CT - [1] - (57) - Property term for vapor (=1 for ideal gas) CT2 - [BTU] - (296) - Property term for wall *CODE - [0] - (86) - = 2 for a new run = 1 for continuation of prior run vs READ FORMAT INPUTJT **CONST - [SECONDS] - (22) - conversion constant for dimensionless time **CONST1 - [OF-1] - (23) - conversion constant for dimensionless temperature **CONST2 - [SECONDS-1] - =.1/CONST **CONST3 - [~F] - = 1/CONST1 D D -[1] -- Temporary storage Dl(J) - [1] - (78) D2(J) - [1] - (80) DI(J) - [1] - (278,374) DRI.(J)- [1] - (173,302) DR2(J)- [1] - (174) DR3(J)- [1] - (175) DR4(J)- [1] - (176) DR5(J)- [1] - (178) DR6(J)- [1] - (179) DR7(J)- [1] - (180) DR8(J)- [1] - (181) I - [1] - (172) 75

D DR - [1] - (14) - Radial grid space D2R - [1] - (32) - DR/2 DR1 - [1] - (63)DR2 - [1] -(15) - DR2 DR3 - [1] - (159) DR4 - [1] - (160) DR5 - [1] - (161) DR6 - [1] - (18) DR7 - [1] - (166) DRZ - [1] - (62) - (DR/DZ)2 **DX - [1] - (58) - Axial grid space-liquid Dxl - [1] - (148) - DT/DX DX2 - [I] - (59) - DX2 DX3 - [1] - (158) DX5 - [1] - (162) DX7 [1 - (164) DX8 - [1] - (165) DXR - [1 - (65) - (DX/DR)2 DRX - [1] - (66) - 1/DXR DYY - [1] - (17) - 2(DR)2 DY1 - [1] - (155) DY3 - [1] - (154) DY5 - [1] - (157) DY7 - [1] - (156) DZ - [1] - (60) - Axial grid space-vapor DZl - Il[] - (149) - DT/DZ DZ2 - [1] - (61) - DZ2 D2Z - [1] - (69) - DZ/2 DZ3 - [El - (150) DZ5 - [1] - (151) DZ7 - [1 - (152) DZ8 - [1] - (153) **DM - [lbm] - (220) - Net phase change in time DT. (minus sign for evaporation) **DMB - [lbm] - (309) - Phase change in time DT from energy balance on tank wall and liquid system (plus sign for evaporation) **DME - [lbm] - (218) - Phase change at L-V interface in time DT (minus sign for evaporation) **DMT - [lbm] - (219) - Cumulative total of DME (minus sign for evaporation) **DMBT - [lbm] - (311) - Cumulative total of DMB **DMTT - [lbm] - (312) - Cumulative total of DM (minus sign for evaporation) 76

D ***DT - [1] Time step interval *DTS - [ F] - (292) Max. permissible liquid superheat at node adjacent to wall *DTW - [ ~F] - (291) Max. permissible superheat of wall in contact with liquid DTDRW(I) - [1] - (205,448,v454) Temperature gradient in fluid at wall DTDXG(J) - [1] - (199) Temperature gradient in vapor at Liquid-vapor interface DTDXL(J) - [1] - (198) Temperature gradient in liquid at liquid-vapor interface DQIL - [1] - (212) Mean temperature gradient in liquid at L-V interface DQILDT - [BTU/sec] - (216) Heat transfer rate in liquid at L-V {ate~race (plus for heat transfer to i.quid **DQIDT - [BTU/sec] - (217) Heat transfer rate in vapor at L-V interface (plus for heat transfer out of vapor) **DQWDT - [BTU/sec] - (214) Heat transfer rate to vapor from wall (plus for heat transfer to vapor) **DWDT - [BTU/sec] - (314) Work rate by vapor control volume DQI - [BTU] - (308) Sum of heat transfer to wall in contact with liquid and to liquid at L-V interface in DT (plus for heat transfer in) DQL - [BTU] - (306) Enthalpy rise of liquid in DT DQS - [BTU] - (299) Enthalpy rise of -wall in contact with liquid, in DT> *DELTA - [feet] Wall thickness E EJ(J) - [1] - (85,282) EJl(J) - [1] - (81) *EPSLON - [1] - (289,385) Maximum fractional change in stream function between iterations before iteration is terminated F F(J) - [1] - (2799,375) Temporary storage FI - [1] - (171,254,3 48) Temporary storage 77

G *G - [ft/sec2] - Acceleration corresponding to effective body force acting on container GAM - [1] - Ratio of specific heats-vapor H H - [0] - (281) - Index **HFG - [BTU/lbm] - (51) - Latent heat of vaporization corresponding to TSAT **HG(I) - [BTU/sec-ft,2-~F] - (454) - Local heat transfer coefficient of vapor based on AT between wall and centerline **HIL(I) - [BTU/sec-ft2- F] - (448) - Local heat transfer coefficient of liquid based on AT between wall and centerline I I - [0] - Axial nodel index number I1 - [0] - Mid vertical height of liquid J J - [0] - Radial nodel index number K *KG - [BTU/sec-ft-~F] - Thermal conductivity-vapor *KL - [BTU/sec-ft-OF] - Thermal conductivity-liquid M *M - [0] - Number of vertical divisions in liquid (M+l = L-V interface) M1 - [1] - (27) - Liquid-wall property M2 - [1] - (28) - Liquid-wall property M3 - [1] - (29) - Vapor-wall property M4 - [1] - (56) - Vapor property M5 - [1] - (31) - Initial temperature M6 - [1] - (223) - Pressure rise ratio in DT M7 - [1] - (30) - Vapor-wall property - [0] - Number of radial divisions (N+1 = wall) *N1 - [1] - (428) - Minimum number of time steps before punch is permitted 78

N N2 - [0] (466,468) Control number for establishing time steps between printouts. N3 - [0 (48, 467) Control number for establishing time steps between printouts *NC - [0] (437) Minimum number of time steps NT before fractional temperature changes between time steps is computed *NE - [0] (270,367) Maximum number of iterations permitted in computations of stream function **NE1 - [0] (269,270) counter on iterations on stream functionliquid **NE2 - [0] (366,367) counter on iterations on stream functionvapor NP - [0] (433, 43 5) control for printout instruction NP1 - [1] (433,434,435) control for printout instruction *NR - [0] (285,381) Index steps of grid spaces for which fractional change in stream function with iterations is computed **NT - [0] (50) Time step number *NMAX - (465) Maximum number of time steps NT permitted before program is ter;minated *NEW - [ft2/sec] kinematic viscosity-liquid *NEWV - [ft2/sec] kinematic viscosity-vapor **NUG - [1] Mean Nusselt number-vapor **NUL - [1] Mean Nusselt number-liquid P *P - [0] Total number of vertical divisions *PO - [lbf/in2] Initial system pressure P1 - [lbf/in ] (49, 224) Pressure previous time step **PRESS - [lbf/in2] (221) Current system pressure **PR - [1] (10) Prandtl number-liquid **PRV - [1] (11) Prandtl number-vapor Q QI - [1] (213) Mean temperature gradient in vapor at L-V interface QW - [1] (208) Mean temperature gradient in vapor at wall *QSURFG - [BTIJU/sec-ft2] Imposed heat flux on container exterior on vapor portion *QSURFL - [IBTU/sec-ft2] - Imposed heat flux on container exterior on liquid portion 79

R **R1 - [1] - (442) - Maximum value of R2 R2 - [1] - (441) - Fractional change in temperature between time steps R3 - [1] - (287,383) Fractional change in stream function with iteration R1(J) - [1] (38) R2(J) - [1] - (79) R3(J) - [1] -(177) R4(J) - [1] - (82) R5(J) - [1] (83) R6(J) - [1] - (84) RO(I,J) - [lbm/ft3]. - (96) - Local vapor density *ROL - [lbm/ft3] Liquid density *ROW - [lbm/ft3 ] Wall density **RP - [1] - (52) Ratio of current pressure to initial **RMAX - [1] - (289) - Maximum value of change in stream function with'iteration-liquid RMAXX - [1] - (385) - Maximum value of change in stream function with iteration-vapor **RAG - [1] - (455) - Mean Rayleigh number-vapor **RAL - [1] - (451) - Mean Rayleigh number-liquid *RGAS - [ft-lbf/lbm-~R] - Gas constant-vapor **RQIN - [1] - (215) - Fraction of total heat transfer rate to vapor part of tank going into vapor RALPHA - [1] - (13) - Ratio-thermal diffusivity of vapor to liquid RNEW - [1] - (12) - Ratio-kinematic viscosity of vapor to liquid S S - [1] - (139,142) - stability criteria SMAX - [1] - (143,144) - Maximum value of S **SF(I,J) - [1] - stream function **SFI(J) - [1] - vapor stream function at L-V inter face T **T (, J) - [1] Current temperature TO(IJ) - [1] - Temperature of prior time step **TAU - [seconds] - (146) - Current time *TACJ1 - [seconds] - (434) - Basic multiple of time for which printout occurs **TIME - [1] - (149) - Current time TINIT - [ ~F] - (21) - Inl.tial saturation temperat'ure corresponding to initial pressure 80

T TSAT - [OF] -(20,225) - Current saturation temperature corresponding to current pressure U **U(I,J) - [1] - Axial component of velocity **U - [1] - (147,203,313) - Mean velocity of liquid-vapor interface UG - [1] - (99,3 15) - Mean velocity of vapor at liquid-vapor interface **UL(J) - [1] - (200) - Local L-V interface velocity due to interface phase change UO(I,J) - [11 - (396) - Temporary storarge of U(I,J) **UR - [1] - (460,462) - Mean axial velocity at mid-point of liquid (serves as check on continuity) V **V(I,J) - [1] - Radial component of velocity VO(I,J) - [1] - (395) - Temporary storage of V(I,J) **VG(J) - [1] - (111,415) - Radial velocity of vapor at L-V at interface W W(I,J) - [1] - Vorticity WO(I,J) - [1] - Vorticity of prior time step X ***X - [1] - (147) - Liquid fraction in container Y Y - [1] - (9) - A/B z Z - [1] - Compressibility factor ZOVRZ - [1] - (70) - = Z/2, = 0.5 for ideal gas 81

APPENDIX D COMPUTER PROGRAM DATA INPUTS 83

P - Total number of vertical grid spaces M - Number of vertical grid spaces in liquid N - Number of radial grid spaces CODE - = 2 for a new run = 1 for continuation of prior run via READ FORMAT INPUT Nl - Minimum number of time steps NT before PUNCH FORMAT DDATA is permitted NC - Minimum number of time steps NT before fractional temperature changes between time steps is computed NMAX - Maximum number of time steps NT permitted before program is terminated NE - Maximum number of iterations permitted in computations of stream function EPSLON - Maximum fractional change in stream function between iterations before iteration is terminated NR - Index steps of grid spaces for which fractional change in stream function with iterations is computed DT - Initial estimate of dimensionless time step TAU1 - seconds - Basic multiple of real time for which printout occurs A - feet - Tank radius APVERB - Ratio of tank radius to height X - Initial fraction of liquid in tank QSURFG - BTU/sec-ft2-Imposed heat flux on container exterior on vapor portion QSURFL - BTU/sec-ft2 - Imposed heat flux on container exterior on liquid portion G -ft/sec2 - Acceleration corresponding to effective body force acting on container PO - lbf/in2 - Initial system pressure DTW - 0F - Imposed maximum permissible temperature difference between container wall and saturation DTS - OF - Imposed maximum permissible temperature difference between liquid nodes adjacent to container wall and saturation CPW - BTU/lbm-~F - specific heat of container wall ROW - lbm/ft3 - Density of container wall ALPHAW - ft2/sec - Thermal diffusivity of container wall DELTA - ft - container wall thickness ALPHA - ft2/sec - Thermal diffusivity of liquid ROL - lbm/ft3 - Liquid density CP - BTU/lbm - ~F - Specific heat of liquid KL - BTU/sec ft - OF - Thermal conductivity of liquid BEITA - ~F -1 Volumetric coefficient of expansion of liquid C5 - Coefficient for computing TSAT C6 - Exponent for computing TSAT NEW - ft2/sec - kinematic viscosity of liquid ALPHAV - ft2/sec - Thermal diffusivity of vapor CPV - BTU/lbm -F - Specific heat of vapor, constant pressure CV - BTU/lbm - ~F - Specific heat of vapor, constant volume KG - BTU/sec - ft - ~F- Thermal conductivity of vapor RGAS - ft-lbf/lbm - ~R Gas constant of vapor Z - Compressibility factor 84

APPENDIX E TYPICAL OUTPUT

S DATA MAP ACCEPT 0000C* DPUNCH 00000* SPRINT 00000* ERROR 00000* SKIP6 00000* SPEEK 00000* SCARDS 00000C* FTRAP 00000* SYSTEM 00000* (PGCM) 10000 (MAIN) IC000.JR4 64412*.IOH 64412* VER.87 64412*.LINK 64412*.IR1 64412*.IR2 64412* (SUBT) 67622.SCOM. 71135* DFMIOH 71466* DFDIOi 71466*.ERR 71541* ~RTN 71541*.ERPRT 71541* ATLOC 71666*.03311 71714*.ICSET 71733*.IOFMT 71751*.RST4 71751*.PRINT 72014*.WR.- 72025*.READ 72041* ~ RD, 72056* ~RDSW 72056* GTCTLW 72140*,GTCLW 72140*.EOF7 72153*.PUNCH 72211*.PC. 72222*.LBUFF 72236*.PRSLT 72367* ~RPDTA 73207* USQRT 74423* SCRT 74473*.EXIT 74515* ELOG 74572*.01300 4.705*.C1301 74773* EXP 75071* ZERO 75204* (PROG) 75240 (ERAS). 77741 EXECUTION BEGINS WITH ROUTINE (MAIN) 02501 LOCS. CAN BE SAFELY USED IN EXPANDING FROG. IOCTAL) FULL CCRE LOADING PROCEDURE USED LADnER- OVERLAPPED BY 05416 LOCS. (OCTAL), FULL CCRE LCADING PROCEDURE USEC A=2.5,AOVERB=.25,P=30,M =2C N=20,X=.67,CC'DE=2,DT=10.E-6, NNMAX=500,Pe=15.,G=32. 2E-5QSURFG=2.78E-4,CSURFL=2.78E-4t DTW=1.,DTS=1.,C5=116. C6=.123,ALPHA=8.35E-7,NEW=1.63E-6, KL=2.33E-5,BEITA=.4SE-2,RCL=7C. CP=.4,ALPHAV=2.14E-5, NEWV=2. 11E-5,KG= 1 4E-6,CPV=. 216,CV=. 153,RGAS=48.3, ALPHAW= 4.16E-5,CPW=.lI,RW=4EE.,DELT=. C1,TAU1=30,Nl=6, NE=50,EPSLON=.01,NC=5 NR=4,Z= 1.0 CONST = 7.485030E+06, CONST2 = 1.3360C00E-07, CONST1 = 4.52832EE+06, CoNST3 = 2.208320E-07 PR = 1.952096, PRV =.985981 NT = 154, lIME = 1.539999E-03, TAU = 1..152654E+C4, CT, i.oeooE-e5 UR = 2.329724, U(1,0) = -2.198956E-C 3, PRESS = 15.2356e8, Z - 1.000000 X =.66S99E, OX =.033500, OME = -3.020815E-C4, D'T -=.024537 DMB =.000000, DMBT =.000000C, CM = -'.020815E-C4, CMTT = -.024537 ROIN =.058085, DQIDT = 6.577250E-05, WUT = 1.625330E-C7,:wT = " 8'T.366078E-04 HFG = 296. 621334 U!1,1)...U(31,21).CCCCOCE-00.0000CE+00.E+.000000E+00.000000E+00.OOOOOE+0.OOOOCOE+CC. OCOOOOE+O.CCOOOCE+O'.000000 +00.OOOOOOE+00.OOOOOE+00.OOOOOE+00 OOOOOOE+00.0000E+00.CCOOOOE+00.CO(0OE+00.CCOCOOE+00.OO OOOE+00. COCOOE+CC.COOO OE+00.OOOOOCE+CO

-2.211137E+01 -2.211628E+01 -2.225188E+01 -2.244571E+01 -2.272659E+C1 -2.310834E+C1 -2.361043E+01 -2.425C56E+01 -2.509219E+01 -2.615838E+01 -2.752779E+01 -2.929880E+01 -3.1 i61301E+01 -3.467842E+C1 -3, 880801E+O 1 -4.44s614E+01 -5.269644E+01 -6.635'45E+01 -5.981325E+01 i.179201E+C2.OOOOOOE+OC.... -4.683610E+01 -4.685800E+01 -4. 7085C5E+C1 -4.T740OC8E+01 -4.785116E+01 -4.845818E+C1 -4.524635E+C1 -5.C24.81E+C1 -5.150318E+01 -5.306329E+01 -5.499125E+01 -5.736472E+O1 -6.027796E+01 -6.384294E+01 -6.818893E+O01 -7.346853E+O0 -7.52CgE+01 -8.827012E+01 -7.0334fCE+C1 2.C37808E+02 OCOOOOE+00 -7.C13079E+01 -7.015323E+01 -7.037743E+01 -7.C68790E+01 -7.113194E+O1 -7.172829E+CI -7.249926E+01 -7. 34713E+G1 -7.467518E+01 -7.614510E+01 -7.791836E+01 -8.003371E+01 -8.252880E+01 -8.'o543651E+Cl -8. 818148E+01 -9.25826E+O -9.682913E+01 -1.00642E+C2 -6.765599E+Ol 2.721730E+02.CGCOOCE+00 -9. 198101E+O1 -9.199214E+01 -9.211.884E+C1 -9.229820E+01 -9.256296E+01 -. 253C65E+C1 -9.342255E+Cl -S.4C6132E+O1 -9.487201E+01 -9. 587981E+01 -9. 710866E+01 -9.857919E+01 -1.003064E+02 -1.022977E+C2 -1.045527E+02 -1.070612E+02 -1.O961e E+C2 -1.092829E+02 -6.147666E+01 3.2S7831E+C2.000C OOE+400 -1.1321C6E+02 -1.131973E+C2 -1.131234E+02 -1.130357E+02 -1.129417E4C2 -1.128634E+C2 -1.128279E+C2 -1. 128651E+02 -1.130063E+02 -1.132824E+02 -1.137217E+02 -1.143478E+02 -1.151781E+02 -1.162229E+02 -1.174879E+02 -1.189573E+02 -1.202147E+C2 -1.165049E+02 -5.393336E+C1 3.O804950E+C2.C0OOOOCE+CO -1.355505E+02 -1. 354936E+02 -1.350659E+C2 -1.345101E+C2 -1.337946E+02 -1.52568E+C2 -1.320439E+C2 -I. l1c5EC2 -1. 302104E+02 - 1. 294025E+C02 -1.287380E+02 -1.282618E+02 -1.280097E+02 -1.280079E+C2 -1.282771E+C2 -1.287s85E+02 -1.289337E+02 -1.222121E+C2 -4.553632E+Cl1 4.26f123E+02.0000C0CE+400 -1.619037E+02 -1.617711E+02 -1. 607155E+C2 -1.593285E+02 -1.575151E+02 -1.553476E+C2 -1.529197E+C2 -1.5C34C2E+; 2 -1.477259E+02 -1.451929E+02 -1. 428477E+C2 -1.407794E+02 -1. 390547E+02 -1. 377175E+C2 -1.367927E+02 -1.362332E+02 -1.351375E+C2 -1.253352E+C2 -3.555222E+C1 4.69C483E+02.OOOOOOE+400 -1.971351E+02 -1.968796E+02 -1.946884E+02 -1.S17840E+02 -1.879560E+02 -1.e3347]E+O2 -1. 78143E+02 -I.725956E+02 (o -1. 669441E+02 - 1.614516E+02 -1.563457E+02 -1.518035E+02 -1.479428E+02 -1.448252E+C2 -1.424738E+02 -1.408340E+02 CC) -1.388662E+02 -1.271841E+02 -2.547719E+01 5.101930E+02.COOOOOE+CC -2.483276E+02 -2.479085E+02 -2.437584E+C2 -2.381751E+02 -2.307176E+02 -2.216422E+C2 -2.113332E+02 -2.CC3C16E+02 -1.891394E+02 -1.784355E+02 -1.686878E+02 -1.602406E+02 -1.532615E+02 -1.477534E+C2 -1.435848E+02 -1.404458E+02 -1.367071E+02 -1.2ClqC5E+02 -8.311547E+CC 5.455704E+02.CCOOOCE400 -3. 232453E+02 -3. 226866E+02 -3.155704E+02 -3.057719E+02 -2.923645E+02 -2.757CC lE402 -2. 564844E+C2 -2.358041E+02 -2.149971E+02 -1.954091E+02 - 1.781428E+02 -1.638970E+02 -1.529398E+02 -i.451S73E+C2 -1.4C4146E+02 -1.382898E+02 -1.374C97E+02 -1.230529E+02 -6. 175CC6E-01 5.894478E+02.OOOOOCE+CC.......... -4.241774E+02 -4.235865E+02 -4.128358E+C2 -3.q75921E+02 -3.759794E+02 -3.48225CE+C2 -3.154311E+02 -2.1 E27E+C2 -2.436424E+02 -2.101698E+02 -1.815030E+02 -1. ~89480E+02 -1.427801E+02 -1.324584E+02 -1.269832E+02 -1.252167E+02 -1.248605E+02 -1.0500C2E+02 5.428528E+01 6.517C09E+02.CCCOOCE+CO -5.4C5803E+02 -5.400C068E+02 -5.257460E+C2 -5.049056E+02 -4.741055E+C2 -4.33C652E+C2 -3.832464E+02 -3.280C22E+02 -2.720399E+02 -2.203770E+02 -1.771503E+02 -1.447378E+02 -1.235249E+C2 -1.123189E+C2 -1.091181E+02 -1.118555E+02 -1.175'97E+C2 -9.864233E+01 1.4620S5E+02 6.995412E+02.OOOOOOE+OO -6.497908E+02 -6.491102E+02 -6. 320017E+02 -6.063197E+02 -5.668796E402 -5.127491E+C2 -4.45756CE+02 -3.70579E02 -2.940208E+02 -2.235700E+02 -1.656213E+02 -1.239646E+02 -9.921260E+O1 -8.937681E+C1 -9.114832E+01 -1.0C11170E+02 -1.14S081E+02 -9.488580E+01 2.2278C2E+C2 7.536559E+02.CCOOCCE+CC __ _ -7.280661E+02 -7.271375E+02 -7.07S985E+02 -6.786406E+02 -6.323588E+C2 -5.679241E+C2 -4. E1621CE+02 -3.7158'E+O2 -3.048780E+02 -2.202081E+02 -1.515213E+02 -1.040198E+02 -7.861821E+01 -7.256481E+Cl1 -8.132342E+01 -1.C02'01E+02 -1.232050E+02 -9.S17574E+01 2.S66782E+02 8.C65073E+C2.COOOOOE+OC -7.595894E+02 -7.584945E+02 -7.385789E+C2 -7.075731E+02 -6.582059E+02 -5.ES5087E+C2 -5.041374E+C2 -4.CE1C42E+02 -3. 100772E+02 -2.200262E+02 -1.471928E+02 -9.783018E+'O -7.337138E-01 -7.C58148E+01 -8.382255E+0i -i.C73184E+02 -1._167S8E+C2 -7.827946E+01l 3.9045C5E+C-2 8.410642E+02.OOOOOOE+OO

-7.365917E+02 -7.357071E+02 -7.171696E+02 -6. 880736E+C2 -6.418908E.02 -5. 1EIE'7E+C2 -4.996360E+02 -4.116660E+02 -3.218854E+02 -2.389192E+02 -1.708284E+02 -1.235208E+02 -9.871720E+01 -9.338325E+01 -1.014688E+02 -1.169637E+02 -1.201232E+02 1.846268E+O1 5.125SC4E+C2 8.235279E+02.CCCOCOE+00 -:.53652CE+02 -6. 531456E+02 -6.389111E+02 -6.163503E+02 -5.805661E+02 -5.1 14536E+C2 -4.711419E+02 -4.C377C3E7 C2 -3.350413E+02 -2.712689E+02 -2.180214E+02 -1. 787402E+02 -1.537808E+02 -1.4C1121E+C2 -1.308327E+02 -1.261681E+02 -3.E78747E+01 2.1414EEE+C2 6.C78245E+02 7.137864E+02.000000E+0C -5.025515E+02 -5.024320E+02 -4.944684E+C2 -4. E14649E+02 -4.604321E+02 -4.3122C5E+C2 -3.c4S327E+C2 -3. 38452E+02 -3. 111886E+02 -2.705911E+02 -2.351936E+02 -2.064881E+02 -1.827778E+02 -1.568788E+C2 -1.131325E+02 -2.789796E+01 1.13E81C2E+02 3.117109F+02 5.186524E+C2 4.901771E+02.OOOOOCE+CC -2.788324E +02 -2.789231E+02 -2.764835E+02 -2. 721145E+02 -2.645342E+02 -2.534676EE+C2 -2.3C0216E+C2 -2.21]443E+C2 -2.020164E+02 -1.808089E+02 -1.583320E+02 -1.341289E+02 -1. 065836E+02 -7.26514EE+C1 -2.794357E+01 3.272106E+01 1.126C64E+02 2.041811E+02 2.7351S7E+02 2.315450E+02 CCOOOCE4CO 5.473984E-01 5.473984E-01 7.298645E-C1 7.298645E-01 7.298645E-C1 7.2C8645E-C1 7.2C8645E-Cl1 7.298645E-C1 7.298645E-01 7.298645E-01 7.298645E-01 7.298645E-01 7.298645E-01 7.298645E-01Cl 7.2S8645E-C1 7.298645E-01 7.2 SE45E-01 7.298645E-01 7. 298645E-01 7.305762E-C1.OOOOOCE+CO -5.388525E+01 -5.38782CE+01 -7.230521E+01 -7.278958E+C1 -7.333054E+C1 -7.373435E+C1 -7.373837E+01 -7.3CCS18E+01 -7.114455E+01 -6. 768447E+01 -6.21330OE+Cl -5.399423E+01 -4.280847E+01 -2.819807E+C1 -9.8S8295E+OC 1. 219124E+01 3.776487E+01 6.496241E+01 e.625,C1lE+C1 6.400463E+01.COOOCCE+CC -1.032714E+02 -1.033163E+02 -1.388776E+02 -1.39S729E+02 -1.411530E+02 -1.42C267E+C2 -1. 4206CeE+C2 - 1.405781E+02 -1.367669E+02 -1.297086E+02 -1.184231E+02 -1. 019304E+02 -7.932722E+01 -4,988296E+C1 -1.318201E+01 3.05665CE+01 7.976737E+01 1.2899% eE+C2 1.61S013E+02 1.132115E+C2.OOOOOOE+00 -1.435387E+02 -1.436568E+02 -1.934623E+C2 -1.953036E+02 -1.S7283CE-t+02 -I.987956E+C2 -1.S90207E+02 -1.S6S252E+02 CD -1 -912928E+02 -1. 807804E+02 -1.639973E+C2 -1.395962E+-02 -1.063793E+02 -6.345126E+01 -1.051333E+01 5.146846E+01 \,0 11.189829E+02 1.826529E+02 2.18S738E+C2 1.4641C4E+02.OOOOOCE4CC -1. 706188E+02 -1. 708263E+02 -2.306679E+02 -2. 334531E+02 -2.365116E+02 -2.3SC278E+C2 -2.3S8944E+02 -2.7728-EE+G2 -2.309281E+02 -2.177695E+02 -1.965354E+C2 -1.656456E+02 -1.237893E+02 -7.011673E+C1 -4.667783E+00 7.C60270E+0I 1.501719E+02 2.212577E+C2 2.557183E+02 1.650289E+C2.COOOOOE+OC -1.803302E+02 -1.8065C8E+02 -2.449.155E+C2 -2.488407E+02 -2.532972E+02 -2.73056E+C2 -2.5S5247E+O 2 -2.E2754E+02 -2. 516187E+02 -2.374921E+02 -2.138921E+02 -1. 790706E+02 -1.317248E+02 -7.122623E+C1 1.876205E+C 8.452577E+01 1.693S5EE+C2 2.413482E+02 2.703626E+C2 1.692162E+02.OOOOOOE+OC -1.6S0696E+02 -1.695218E+02 -2.311775E+02 -2.362296E+02 -2.421800E+02 -2o4800C8E+C2 -2.22555E+C2 -2.31150E+02 -2.484290E+02 -2.35865CE+02 -2.131226E+02 -1.782149E+02 -1.297996E+02 -6.757025E+C1 7.124181E+00 9.002655E+01 1.72523gE+02 2.3868E5E+02 2.6CC340E+C2 1.584337E+02.CCOOCCE+CC -1.355081E+02 -1.360790E+02 - 1.871253E+02 -1. 927695E+02 -1.996694E+02 -2.C69481E+C2 -2.133281E+02 -2.171215E+02 -2. 162533E+02 -2.083408E+02 -1. 908680E+02 -1. 615086E+-02 -1. 186343E+02 -6.2C0152E+01 6.426422E+0C 8.152475E+01 1.543793E+02 2.1028C5E+02 2.25C899E+02 1.348011E+02.COOOOOE+CO -8. 3E7775E+01 -8.452203E+01 -1.176191E+C2 -1.225290E+02 -1. 287539E+02 -1. 3579C08E+C2 -1.4283c1E+C2 -1.487704E+02 -1.520871E+02 -1.508894E+02 -1.428992E+02 -1.256437E+02 -9. 69359 1E+01 -5.575205E+Cl -3.438022E+00 5,5CS378E+01 1.11753SE+C2 1.548514E+02 1.665842E+C2 1.C00982E+C2 oOOOOOCE+00 -2.E7388CE+01 -2. S364C3E+01 -4.163849E+C1 -4.404742E+01 -4.720672E+01 -5.102794E+C1 -5.53G838E+01 -5.17261CE1 -6.379774E+01 -6. 679627E+01 -6. 763685E+01 -6.478414E+01 -5,630510E+01 -4.02742 lE+01 -1.578596E+01 1. 528582E+01 4.768112E+01 7.417348E+01 E.5C6414E+01 5.398312E+01 COOOCCE+CC.CCOOOCOE+00.OOOOOOE+O00.OOOOCOE+C.COOOOOE+00.OOOOOOE+CC.CCCCCCE+CC CCOOCCE+00.CCCCCCE+0C.OCOOCOE+00.OOOOCCE+.OOOOOOO00.OOOE+00 COOOO0E+C0.CC0000E+00.COCCCE+C.oo ooE.ooce.CCOOCCE+00 0OCCCCE+CC.OOOOOOE+oCC.CCCOCOE+00.000000E+00

V(1,1)...V(31,211.CCOCOE+OQ.OOOOOOE+OO0.OO.OOOE+OO.OOOOOE+OO.OOOOOCE+OO.COOOOOE+CC.OCOOOCE+OC.CCOOOCOEOO.COOCOCE+00.000000E+00 OOOOOOEO COOOOOE+ OOOOOOE+00.0 00+0. O+ 0.CCOOOOE+CO.000000E00 C.CCOOE+CE.CCOOCCE+OO.COOOCOE+00.OOOOOE+OO.COOOOOE+OO.OOOOO.E+.O.CCOOOOE+OO 1.797913E+01 3.605490E+01 5.427706E+01 7.272249E+C1 9.14782IE+C1 1.1C6439E+02 1.303354EC2 1.506878E+02 1.718600E+02 1.940366E+02 2.174298E+02 2.422754E+02 2.688196E+02 2.972836E+02 3.277821E+02 3.600955E+02 3.926962E+02 4.165739E+C2 3.752857E+02.COCOOOCE+00.COOCCOE+00 1.817011E+01 3.63 e3 8C+O -01 5.465256E+01 7. 300314E+01 9. 146446E-+C1 1.100661E+02 1. 288359EO02 1. 477965E+02 1.66c581E+02 1.863070E+02 2.057866E+02 2.252661E+02 2.444858E+02 2.629933E+02 2.7596'E+02 2.o38736E+02 3.0137696+02 2.922639E+02 2. 389425E+02.OOOOOOE400. CCOOOOE+OO 1.67S4E83E+C01 3.355145E+C1 5.024855E+01 6.68571CE+01 8.234622E+C,1 9..68163E+O0 1.158229E+2 1.317190E+02 1.473027E+02 1.624817E+02 1.771288E+02 1.910681E+02 2.040611E+02 2.157917E+02 2.258431E+02 2.335592E+02 2.370666E+02 2.280949E+C2 1.760413E+02.OOOOOCE+CC.0OOOOCE+00 1.591571E+01 3.171957E+01 4.7351C8E+01 6.273177E+01 7.177846EE+C1 9.243477E+C1 1.CE6C82E+02 1.202308E+02 1.332259E+02 1.455114E+02 1.569969E+02 1.675813E+02 1.771523E+02 1. E55864E+02 1.927354E+02 1.982667E+02 2.006018E+02 1.918786E+C2 1.435951E+02 *CCOOOOE+OO.CCOOOOE+OO 1.596c15E+01. 173353E+C1 4.718525E+01 6.218840E+01 7.661750CE+C1 9.036083E+01 1.C33226E+02 1.154239E+02 1.266028E+02 1.368135E+02 1.460248E+02 1.542193E+02 1.613920E+C2 1.675497E+02 1.726900E+02 1.766260E+02 1.,77S213E+02 1.691638E+02 1.237354E+C2.0000OOE4O0.COOOOOE+OO 1.7686C6E+01 3.501280E+01 5.1792C8E+01 6.779265E+01 8.281035E+01 9.667588E+Cl 1.CS2609E+02 1.204823E+02 1. 303025E+02 1.387284E+C2 1.458051E+02 1.516087E+02 1.562375E+02 1.5S80OCE+02 1.623805E+02 1.638483E+02 1.628217E+02 1.523781E+C2 1.C88075E+C2.CCCOCCE+CC o.OOOCCE+OO 2.223636E+01 4.382790E+C1 6.442394E+01 8.360399E+O1 1.C10071E+C2 1.163475E+02 1.2S4315E+C2 1.401694E+02 1.485796E+02 1.547798E+02 1. 589686E+02 1.614004E+02 1. f235S8E+C2 1.6214C8E+02 1.610212E+02 1.591273E+02 1.5556q5E+02 1.4451C0E+C2 1.021438E+02.CCOOOCE+00.CCOOOOE+00 3.132454E+01 6.151852E+01 8.989685E+0C1 1.15657CE+C2 1.381066E+C2 1.567C24E+C2 1.7110C E+02 l.812234E+02 1.872278E+02 1. 895190E+02 1. 886605E+02 1.852918E+02 1.ECC387E+C2 1.734282E+02 1.57c46E+02 1.57C658E+02 1.458218E+02 1. 25S425E+02 8.33:03CE+C1. OOOOOOE+CC.COOOOOE+OO 4.643CSSE+01 9.1088(7E+01 1.326891E+02 1.697371E4C2 2.C08586E+C2 2.249270E+02 2.412485E+02 2.497205E+02 2.50880gE+02 2.458202E+02 2.359900E+02 2.229793E+C2 2.C83248E+02 1.934004E+02 1. 793S43E+02 1.672834E+02 1.56S816E+02 1.4236e1E+C2 1.106457E+02. CCOOCE+o00. CCOOOOE+OO 6.612913E+01 1.299616E+02 1.894636E+02 2.421613E+02 2.854941E+*02 3.11701E+C2 3.356802E+02 3.407357E+02 3.334292E+02 3. 160430E+02 2.915811E+02 2.632056E+02 2.337640E+C2 2.055158E+02 1.600576E+02 1.58255E+0C2 1. 386515E+02 1.045679E+C2 -7. 327448E+00.000OOOE+OO.00000E+OO E.316C54E+01 1.63c1C2E+C2 2.398027E+02 3.074459E402 3.628215E+C2 4.020263E+C2 4.2227C2E+02 4. 227358E+02 4. 049808E+02 3.727080E+C2 3. 309638E+02 2.850687E+02 2.396650E+02 1.981469E+02 1.624870E+02 1.331311E+02 1.C69411E+02 t.023417E+Cl -8.506C4E+01.OOOOOCE+CC.CCOOCCE+OO 8.707033E+01 1.71S168E+02 2.524018E+02 3.246708E+02 3.E372ClE+C2 4.245281E+02 4.433533E+02 4.388508E+02 4.127171E+02 3.695947E+02 3.161538E+02 2.596420E+02 2.-064842E+C2 1.614406E+02 1.274284E+02 1.055203E+02 9.265724E+C1..2905C6E+01 -5.734726E+01.OOOOOOE+OO.CCOOOOE+00 7.241778E+O1 1.4286S0E+C2 2.099681E+02 2.699878E+C2 3.1814866+C2 3.498324E+02 3.616535E+02 3.523636E+02 3.235230E+02 2. 76525E+02 2.275508E+02 1.748180E+02 1.282046E+C2 9.261472E+01 7.0S0565E6O1 6.:6023EE+C1 (. 5C3634E+01 4.023e41E+C1 -1.089186E+02 000.00E+CC.COCCCOE+00 4.243656E+01 8.348672E+01 1.225113E+02 1.567647E+02 1.832560E+C2 l.9c4633E+02 2.C36689E+02 1.952523E+02 1.750988E+02 1.459855E+02 1.125707E+02 8.059971E+01 5.529744E401 4.020561E+01 3.692125E+O1

4.44741E+01 5.365844E+01 1.58c166E+01 -1.28'977E+C2.000000E+CO.COOOOOE+OO 4.001483E+00 8.0874'SE+CC 1.229838E+01..610O74E+C1 1. S4337E+C1 2.419991E+01 2.'58566E+1 3.629777E+01 4.393424E+01 5.168775E+01 5.900409E+01 6.635274E+01 7,.447852E+C1 8.341156E+01 9.161104E+O1 9.408867E+01 7.3177C8E+01 -6.125195E+01 -1.81671OE+02.CCOOOCE+00.CCOOCCE+00 -3.868078E+01 -7.5089C5E+01 -1.079223E+C2 -1.344901E+C2 -1.512252E+C2 -1. 4549EE+C2 -1.4181C1E+02 -1.120631E+02 -6.667315E+01 -9.568711E+00 5.308745E+01 1.139443E+02 1.661068E+C2 2.037071E+02 2.202396E+02 2.C7S497E+02 5.676240E+01 -2. CS94771E+C2 -2.392325E+C2.OOOOOCE+CC.COOOOOE+00 -8. 654C35E+01 -1.689059E+02 -2.441496E+02 -3.071017E+02 -3.520374E+C2 -3.738814E+02 -3.692818E+02 -3.374649E+02 -2. 807693E+02 -2.047398E+02 -1.178885E+02 -3.147351E+01 4.C73254E+C1 7.S51935E+01 5.433086E+01 -7.572167E+01 -3.C089484E+C2 -4.758461E+02 -2.976958E+02. CCOOOOE+00.CCOCCCE+00 -1.4067COE+02 -2. 767917E+C2 -4.047851E+02 -5.188198E+02 -6.1285S7E+C2 -6.818697E+02 -71.22S282E+02 -7.362C85E+02 -7. 255982E+02 -6.986833E+02 -6.,66 C067E+02 -6. 398136E+02 -6. 327979E+02 -6.567641E+02 -7.201281E+02 -8. 2518B56E+02 -8.501942E+02 -6.8618C3E+02 -2. S17352 E+C2.000000E+00.000000E+00 -1.896971E+02 -3.763765E+02 -5.572607E+02 -7.276717E+C2 -8.826922E+C2 -1.017989E+C:3 -1.130626E+03 -1. 219845E+03 -1.287602E+03 -1. 338571E+03 -1.379340E+03 -1.416343E+03 -1.451701E+C3 -1.475672E+03 -1.456473E+03 -1. 337459E+03 -1.091213E+03 -7. CllOCIE+C2 -2.404403E+02.0 COOOCE+CC.CCOOOOE+00 -2.149394E+02 -4.286857E+C2 -6.396252E+02 -8.445816E+02 -1.C3SE77E+C3 -1.221690E+03 -1.3,E37SE+C3 -1.530739E+03 -1.652030E+03 -1.747620E+03 -1. 814185E+03 -1.846498E+03 -1.836060E+0C -1.77C14CE+03 -1.632076E+03 -1.402658E+03 -1.0601SlE+C3 -6.395190OE+C2 -2.16C580E+C2.000000CE+00.CCOOOOE+00 3.980614E+01 7.986476E+C1 1.203395E+02 1.613177E+C2 2.027352E+C2 2.4431C4E+02 2.E54742E+02 3.253002E+02 3.624434E+02 3.950967E+02 4.209737E+02 4.373268E+02 4.41C090OE+02 4.285964E+02 3.966111E+02 3.420188eE+02 2.633540E+02 1.6f39575E+C2 6.C80399E+01.00000OE+00 O.COOCCE+00 3.444746E+01 6.9254C9E+01 1.045782E+C2 1.405059E+02 1.769641E+C2 2.1364S5E+02 2.4SS566E+02 2.849144E+02 3.171420E+02 3.448324E+02 3.657638E+02 3.773385E+02 3.766593E+02 3.6C6926E+C02 3.266378E+02 2.727886E+02 2.000551E+02 1.151421E+02 3.637832 E+01 000000E+CC.CCOOOOE+00 2.602231E+01 5. 248262E+C1 7.957383E+01 1.074341E+02 1.36C316eE+C2 1.650905E+02 1.S4C211E+02 2.218672E+02 2.472875E+02 2.685802E+02 2. 837467E+02 2.905806E+02 2.E67782E+C2 2.7C1239E+02 2.389229E+02 1.929744E+02 1.349641E+02 7.265779E+01 2.102769E+01.CCOOOCE,00.CCOOOE+00 1.442435E+01 2.938584E+01 4.511662E+C1 6.184262E+C1 7.963665E+01 9. E3456CE+Cl 1.175219E+02 1.363748E+02 1.537600E+02 1.682225E+02 1.780979E+02 1.816580E+02 1.772745E+C2 1.636082E+02 1.39S774E+02 1.C73195E-+02 6.909727E+01 3.235465E+Cl 7.129430CE+00.000000CE+CC.COOOOE+00 -2. 775363E-01 5.090634E-02 1.212859E+00 3.522747E+CC 7.;04759E+CC 1.23358EE+Cl 1. E7632E+01 2.616520E+01. 3.378378E+01 4.065572E+01 4.556467E+01 4.722914E+01 4.456714E+C1 3.7C2926E+01 2.502154E+01 1.049882E+01 -2.815433E+00 -9.726333E+CC -6.659998E+00. COOOOOE+00. C COCOOE+00 -1.716295E+01 -3.381292E+C1 -4.983201E+01 -6.489109E+C1 -7.E62297E+C1 -9.063675+C1 -1.CC5472E+02 -1.080282E+02 -1.128973E+02 -1.152216E+02 -1.153902E+02 -1.140525E+02 -1.118095E+C2 -1.086485E+02 -1.C33523E+02 -9.35325E+Cl01 -7.633154E+01 -5.007C75E+C1 -1. 883813E+C1.000000E+00.CCOOCOE+00 -3.345687E+01 -6.679347E+01 -1.002028E+02 -1.335400E+02 -1.664423E+C2 -1.982659E+02 -2.280391E+02 -2. 544612E+02 -2.759855E+02 -2.910132E+02 -2.981810E+02 -2.966315E+02 -2.E60587E+02 -2.663245E+C2 -2.367891E+02 -1.961341E+02 -1.4414CgE+02 -8.426626E+C1 -2.803569E+01.000000E+C0.0000CCCE+00 -4.325045E+01 -8.705460E+01 -1. 321944E+02 -1.790024E+02 -2.2749GCCE+C2 -2.771912E+02 -3.26SS73E+02 -3.750058E+02 -4.184328E+02 -4. 536587E+02 -4.764983E+02 -4. 827701 E+02 -4.691332E+C2 -4.339582E+02 -3.177586E+02 -3.027076E+02 -2.127670CE+02 -1.176122E+02 -3.644418E+01.000000E+00.000000CE+00 -3.554796E+01 -7.242999E+01 -1.112670E+02 -1.527822E+02 -1. S74746E+C2 -2.454674E+C2 -2.93295E+02 -3. 488698E+02 -4.009089E+02 -4.49C613E+02 -4. 886143 E+02 -5.136775E+02 -5. 178313E+02 -4.954398E+02 -4.434598E+02

-3.63C791E+02 -2.597532E+02 -1.4635312E+C2 -4.68S714E+0 1.OOOOOCE+CC.CCCOCCE+OC.OOOOCCE+CO.OOOOOOE+O.COCOOOE+00.OOOOOE+OO.CCCCCCE+CO.CCOCCCE4OC.CCOCCCE+aC.CCOOOOE+OO 00 OOOE+ OOOOOOF+00.OOOOOOE+00.OOOOOOE+00.COOOOE+Co.CCOOCCE+OC.OOOOOOE+00 OO.OCOCCCCE+00.COOOOOE+OO.CCCCCOE+CC.OOOOOOE+00.OOOOOoE+OC T(1,1)...T(31,21) 4.30060EE-C5 3.495397E-06 3.406996E-06 3.707045E-06 4.161366E-06 4.E2CISCE-C6 5.768783E-06 7.252350E-06 1.21740gE-05 7.695154E-05 1.271062E-C3 2.216198E-02 3.537634E-01 5.C631386E+CC 6.4C864OE+01 7.C5S813E+02 6.633168E+03 5.177224E+04 3.23S763E+05 1.546724E+06 5.185355F+06 7.674411E-04 7.783252E-04 8.193104E-04 8.853089E-04 9. e33987E-C4 1.123774E-C3 1.321699E-c3 1.6CC3C3E-03 1.9S5269E-03 2.564005E-03 3.420809E-03 5.094375E-03 1.329515E-02 1.C03278E-C1 1.13866EE+OC 1.3C9884E+O1 1.5C663CE+02 1.S12554E+03 3.C2518SE+C4 6. 18E326E+05 4.690227E+06 2. 67184CE-02 2.908460E-02 3.030168E-02 3.225529E-C2 3.50S364E-02 3.S034CgE-C2 4.4383S6E-02 5.157609E-02 6.122528E-02 7.421677E-02 9.186434E-02 1.166033E-01 1.592966E-O1 3._06734E-O1 1.336052E+00 1.985418E+Ol 2.20227CE+02 2.5454C.E+03 3.248546E+C4 4.775746E+05 4.453294E+C6 5.248538E-01 5.306C73E-01 5.471323E-C1 5.738845E-O1 6.119075E-01 6.f 313CEE-CI 7.301706E-01 E.1644S3E-CI 9. 264603E-Oi 1.066175E+00 1.243558E+C0 1.470293E+00 1.777847E+00 2.411qC7E+CC 6.075274E+00 4.155294E+O1 3.955206E+02 3.E6S63CE+03 2.966219E+C4 4.576776E+C5 4.385143E4C6 6.422473E+00 6.481852E+C0 6.617509E+00 6.843161E+00 7.156796E+C0 7.566233E+CC E.CE2343E+CC e.718C4CE+O0 S.488847E+00 1.041438E+01 1.151933E+01 1.283687E+01 1.443739E+01 1.67614CE+C1 2.434161E+01 8.310969E+01 6.177296E+02 5.335430E+C3 4.75gC61E+04 4.700313E+05 4. 388604E4C6 \~0 5.9z7374E+01 5.'93C47E+O 1 6.0616?SE+Cl 6.184987E+Cl1 6.352008E+01 6. 614e7 E+ C 1 6.813218E+01 7.1C6361E+01 R) 7.439723E+01 7. 813045E+01 8.227758E+01 8.68R036E+01 9.205595E+01 9.849T64E+01 1.129650E+02 2.010087E+02 9.353266E+02 6.'3C2E5E+C3 5.570224E+C4 4.q36915E+05 4.420132E406 4.4C6950E+02 4.434953E+02 4.448417F+02 4.484357E+02 4.53C37EE6C2 4.582313E+C2 4.626433E+C2 4.688275E+02 4.733346E+02 4.7686C6E+02 4.793369E+02 4.809604E+0O2 4.822134E+02 4.84465SE+C2 4.973095E+02 6.C75C71E+02 1. 53693E+C3 8.768130E+03 6.392758E+C4 5.223C42E+05 4.462839E+06 2.664961E+03 2.67938CE+03 2.669586E+03 2.~64593E+03 2.6555;5E+03 2.3898e6E+C3 2.611385E+03 2.56S306E+03 2.510136E+03 2.433354E+03 2.341024E+03 2.237381E+03 2.127870E+03 2.C18793E+02 1.92568CE+0C3 1.569596+03 2.S67320E+03 1.117732E+04 7.17C4ECE+C4 5.50S611E+05 4.51C352E+C6 1.306714E+C4 1.313308E+04 1.303303E+C4 1.292740E+04 1.276979E+04 1.25:65CE+C4 1.220414E+04 1.17522CE+C4 1.117195E+04 1.047359E+04 9.686345E+03 8.851478E+03 8.012171E+03 7.205527E+03 6.468265E+03 5.S314SE+03 6.544982E+03 1. 115F2E+C4 8.101297E+C4 5.852914E+05 4.564790E*06 5.054096E+04 5.081729E+04 5.043146E+C4 4.S98930E+04 4.92c865E+04 4.E22071E+C4 4.66C716E+04 4.433E67E+04 4.'138460E+04 3.783737E+04 3.389757E+04 2.981866E+04 2.584144E+04 2.214681E+C4 1.884830E+04 1.6C4$87E+04 1.485732E+C4 2.1996S4E+04 8.63S72sE+C4 6.006047E+05 4.631771E40( 1.4e4G46E+05 1.494C66E+05 1.490015E+05 1.486647E+05 1.477733E+05 1.456681E+C5 1.4157CSE+05 1.348253F+C5 1.252014E+05 1.130617E+05 S.928182E+04 8.497587E+04 7.117965E+04 5.863324E+04 4.772458E+04 3.864c11E+04 3.228273E+04 -. 7685CE+C4 1.C32260E+05 6.882082E+C5 4.78657CE+C6 3.227506E+05 3.255672E+05 3.276763E+C5 3.313884E+05 3.352458E+C5 3.7C738E+C5 3.339841E+05 3.232280E+05 3.033167E+05 2.747695E+05 2.400400E+05 2. 026592E+05 1.660939E+05 1.328873E+C5 1.043622E+05 8.CSCC60E+04 6.113S6E+C4 5.832360E+04 1.433713E+C5 1.082424E+06 5.167368E+06 5.3C7542E+05 5.366651E+C5 5.462996E+05 5.624808E+05 5.829203E+C5 6.C28320E+C5 6.1501CEE+C5 6.117041E+05 5.E73684E+05 5.410229E+05 4.768951E+C5 4.028932E+05 3.277559E+05 2.584516E4C5 1.989384E+C5 1.504593E+05 1.1313C5E+05 9.221266E+04 1.8449C9E+C5 1.34590C9E+06 5.558392E+C6

7.090756E+05 7.188928E+05 7.402508E+05 7.76SC53E+05 8.263311E+4C5 8.EO76S9E+C5 9.272671E+C5 $.5C5715E+05 9.378320E+05 8.834446E+05 7.917186E+05 6.755515E+05 5.515182E+05 4.341975E+C5 3.326561E+05 2.502038E+05 1.858518E+05 1.385S7SE+C5 2.354474E+C5 1.766377E+06 5.934748E+4C6 8.3843C5E+05 8.524022E+05 8.869792E+C5 9.476352E+05 1.032024E+06 1.129638E*C6 1.22196EE+*C 1.2E6C98E+06 1.300738E+06 1.253270E+06 1.145221E+06 9.927960E+05 8.210281E+05 6.545376E+C5 5.096433E+05 3.92239CE+05 2.S595S2E+C5 2.035455E+05 2.c,35527E+05 2.178164E+C6 5.935200E+C6 9.331344E+05 9.511 18E+05 9.981679E+C5 1.C81479E+06 1.198567E+06 1.36762E+Cf 1.473;57E+06 1. E25CgEO06 1.635404E+06 1.613903E+06 1.514497E+06 1.351832E+06 1.154685E+C6 9.567557E+05 7.809468E+05 6.322913E+05 4.852531E+05 2.985756E+C5 5.400638E+C5 2.655721E+06 5.S3520CE+C6 1.C10045E+C6 1.031529E+06 1.088544E+06 1.188565E+06 1.328293E+0C 1.4c4C79EC6 1.663620E+Cf 1,8C974gE+06 1.905952E+06 1.932941E+06 1.884170E+06 1.766743E+06 1.601143E+06 1.41786CE+06 1.233248E+06 1.043744E+06 7.85~1053E+05 3. 939107E+05 1.116459E4C6 3.1937C3E+C6 5.935200E4+06 1.080623E+06 1.1035c0E+C6 1.164412E+C6 1.269827E+06 1.416635E+C6 1.9323?E+0C I.781239E+O0 I.S5EC6-E+06 2.100835E+06 2.191275E+06 2.219856E+06 2.186750E+06 2.098170E+06 1.958723E+C6 1.766119E+06 1.520481E+06 1.C55600E+06 1.179014E+06 2.C36:CCE+C6 3.798874E+06 5.93520CE+C6 1.154258E+06 1.175066E+06 1.22S921E+CE 1.3245E5E+C6 1.45806CE406 1.f24365E+CE 1.813165E+C6 2.ClC3qCE+Oe 2.199476E+06 2.363194E+06 2.485310E+06 2.551067E+06 2. 546392E+06 2.459478E+06 2.2S6918E+06 2.131CS'E+06 2.142720E+06 2.4225C8E+C6 3.136242E+C6 4.410170E+06 5.93520CE+C6 1.245637E+C6 1.25912CE+06 1.2946C2E+C6 1.356413E+06 1.44f475E+C6 1.565582E8C6 1.713427E+0(E 1.E8E371E+06 2.087346E+06 2.305938E+06 2.538417E+06 2.777359E+06 3.012465E+06 3.228357E+C6 3.401934E+06 3.5C0320E+O6 3.5S48311E+06 3.825617E+06 4.271677E+C6 5. C14033E+06 5.935200E406 1.406871E+06 1.406E7lE+06 1.406871E+06 1.406871E+06 1.406871E+06 1.406871E+Cf 1.4C68?1E+CE 1.4C6871E+06,,C) 1.406871E+06 1. 406871E+C6 1.406871E+06 1.406871E+06 1.406871E+06 1.406871E+06 1.4C6871E+06 1.406871E+06 1.4C6871E+06 1.406871E+06 1.406e71E+06 1.406871E+06 8.C098e308E+C6 2.757322E+06 2.761148E+06 2.773917E+0C6 2.782206E+06 2.790452E+06 2.ECC415E+C: 2.E13239E+C6 2.E2c,46E+06 2.851 579E+06 2.879268E+06 2.914367E+06 2.958854E+06 3.016393E+06 3.CS4863E+C6 3.212158E+C6 3.408915E+06 3.798092E+06 4.64G331E+06 6.451624E+06 1.013947E+07 1.626798E+07 3.207775E+06 3. 210762E+06 3.22C469E+06 3.230277E+06 3.243468E+C6 3.262375E+C6 3.2S392E+CE 3. 327132E+06 3.378359E+06 3.445959E+06 3.533225E+06 3.644864E+06 3.789244E+06 3.S82990E+06 4.260600E+06 4.6S4864E+06 5. 45e455E+a6 6.502468E+C6 S.6M1520E+06 1.442168E+07 2.107286E+.C7 3.327484E+06 3.33025'E+06 3.339277E+C6 3.353029E+C6 3.374766E+06 3.40E081E+C: 3.457124E+06 3.526258E+06 3.619502E+06 3.7401'7E+06 3.891572E+06 4.078997E+06 4.314502E+06 4.6241C1E+0C 5.059841E+06 5.722315E+06 6.808566E+06 8.683136E+C6 1.1'C651E+C7 1.704966E+07 2.369652E+07 3.3605C4E+06 3.364249E+06 3. 376527E+06 3.398260E+06 3.43421EE+C6 3.490326E4C06 3.5i334CE+C6 3.68'745E+06 3.844240CE+06 4.038626E+06 4.272488E+C6 4.547213E+06 4.874233E+06 5. 286S79E+C6 5.855696E+06 6.7C73S4E+06 8o0421C1E+C6 1.019902E+07 1.364171E+C07 1.e74366E+07 2.503187E+07 3.375419E+06 3.38C844E+06 3.39E86CE+C6 3.431861E+06 3.487733E+06 3.577.,1E+4C6 3.7C6410E+C6 3.889290E+06 4.129902E+06 4.427300E+06 4.774237E+06 5.162592E+06 5.595781E+06 6.107355E+C6 6.782237E+06 7.777661E+06 S.282379E+06 1.159067E+07 1.505635E+C7 1.989007E+07 2.566S46E+07 2. 386461E+C6 3.393946E+06 3.418422E+06 3.465375E+06 3.545480E+06 3.672663E+0C 3.E62455E+06 4.128e682E+C6 4.479012E+06 4.910935E+06 5.410870E+06 5.959756E+06 6.547943E+06 7. 199447E+C6 7. S9509E+06;.1C7521E+06 1.070154E+07 1. 3020(f5E+C7 1.631730Et+07 2.073170E+07 2.595989 E+07 3.. 394C1-3E+06 3.403259E+06 3.433345E+C6 3.492558E+06 3.594916E+06 3.75. 207E+C6 4..C68C4E+0d 4.351545E.06 4.824174E+06 5.406989E+06 6. C91371E+06 6.851277E+06 7.661646E+06 8.521216E+06 9.482237E+06 1.067147E+07 1. 263(:E+07 1.4462C9E+07 1.746120E+C7 2.139431E+07 2.609094E407

3.39665UE+C6 3.405487E+06 3.437397E+06 3.501316E+C6 3.613208E+06 3.755233E+0- 4.073756E+C6 4.475427E+06 5.021635E+06 5.722715E+06 6.573977E+06 7.555459E+06 8.636866E+06 9.788871E+06 1.100162E+07 1.230784E+07 1.3875C8E+C7 1.588638E+07 1.8523f6E+6C7 2.195796E+07 2.615204E+C7 3.38908eE+C6 3.397350E+06 ".4242e2E+C6 3.478352E'06 3.573698E+C6 3.73C77E+C6 3.975563E+O 4.3372C9E+06 4.844132E+06 5.5194COE+06 6.376798E+06 7.418684E+06 8.635669E+06 1.000723E+C7 1.150276E+07 1.308518E+07 1.479170E+07 1.677615E+C7 1.921145E+C7 2.231454E+07 2.617683E4C7 3.378C47E+06 3.383307E+06 3.400382E+06 3.433839E+06 3.491631E+06 3.585548E+CE 3.7:14C1E+06 3.S4E68fE+O6 4.259613E+06 4.687600E+06 5.255606E+06 5.984868E+06 6.894636E+06 8.CO3653E+C: 9.334601E+06 1.C92369E+7 1.283851E+07 1.518814E+07 1.811242E+C7 2.177403E+C7 2.617446E+07 SF{lvl...SF( 31t21).CCCOOCE+00 OOOOCOE+00.OOOOOOE+00.000000E+00.OOOOOCE+00.COOOOCE+OC.CCOOCCE+OC.CocGOcE+00.CCOOOOE+00.OOOOOOE+00.OOCOCCE+CC.COOOOOE+00.OOOOOCE+CO.COOOOOE+CC.CCOOCOE+00.CCCCCCE+O0.COOOOOE+00.OOOOCOE+00.CCCOCCE+CC.COOOOOE+00.OOOOOCE+CC.CCOOOOE+00 -2.763921E-02 -1. 108687E-Cl -2.505503E-01 -4.481759E-C1 -7.C5S832E-C1 -1.027141E+0C -l.416C17E+CC -1. 878569E+00 -2.42290 1E+00 -3.06CC85E+00 -3. 805385E+ 0O -4.680109E+00 -5. 114486E+CC -6.952208E+OC -8.457559E+O0 -1.032565E+01 -1.267874E+Cl -1.5427C4E+C1 -1.542543E+01.CCOOOE+400.COOOOOE+00 -5. 854513E-02 -2.347553E-Cl -5.300263E —C1 -9.467464E-01 -1.488467E+CC -2. 1601 3E+CC -2.S68442E+4 0 -3..92221CE+00 -5.032993E+00 -6.315781E+CO -7.789923E+00 -9.480276E+00 -1. 141859E+C1 -1.364505E+01 -1. 620:39E+01 -1.916731E+01 -2.253579E+C01 -2.5E87q29E+01 -2.488746E+01.OOOOO0E+OC.OOOOOOE+00 -8.76634'E-02 -3.512246E-CI -7.920587E-01 -1.412516E+00 -2.2160O6E+CC -3.2C7607E+OC -4.3S3541E+C -5.782291E+00 -7.384573E+CO -9.213820E+00 -1.128660E+01 -1.362300E+01 -1.624692E+Cl -1.918603E+01 -2.247054E+1 -2.612373E+01 -3.OO~E1lE+Cl -3.363662E+Cl -3.138792E+01.COCOCCE+CO.CCCCCOE+00 -1.149763E-01 -4.602233E-01 -1.036529E+00 -1.845290E+CO -2.E88671E4CC -4.169807E+CC -5.6?3040E+OC -7. 464178E+00 -9.490761E+00 -1. 178230E+01 -1.435046E+01 -1.720916E+01 -2.C37452E+C1 -2.386446E+01 -2.76S6E7 0 -:. 18752CE+01 -3.628284E+01 -3.993589E+C1 -3.644811E+01.OOOOOOE+00.OOOOOOE+00 -1.415132E-01 -5.658488E-01 -1.272593E+00 -2. 26123CE+00CO -3.5213_1E+CC -5.C82917E+00 -f.S16571E+0C -9.033966E+00 -1. 143808E+01 -1. 413350E+01 -1.712661E+01 -2.042570E+01 -2.404093E+01 -2.798392E+01 -3.226537E+O0 -3.6E7636E+01 -4.164970E+01 -4.535546E+C1 -4.C68C66E+01.CCOOOCE+CC.CCOOOCE+00 -1.6c4381E-01 -6.766260E-01 -1.519038E4+00 -2.692842E400 -4.193256E+CC -6.015022E+CC -8.152CS9E+CC -1.060256E+01 -1.336059E+01 -1.642554E+01 -1.979811E+01 -2.348138E+01 -2.748089E+01 -3.180431E+01 -3.645852E+01 -4.1426CCE401 -4. 4'5 C 7E+C1 -5.02Ce82E+Cl -4.439052E+01.000000CE+CO.COCOCCE+00 -2.023796E-01 -8. 067628E-01 -1.806871E+00 -. 193003E+00 -4.S52732E+CC -7.071844E+CC -S.535461E+OC -1.232928E+01 -1.544080E+01 -1. 88E034E+Cl1 -2.258183E+01 -2.660318E+01 -3.C02630E+01 -3.555648E+01 -4.C49845E+01 -4.572743E+C1 -5.C98893E+01 -5.461434E+C1 -4.767278E+01.OOOOOOE+00.OCOOCE+00 -2.4641 88E-01 -9.800253E-01 -2.187701E+00 -3.848957E+00 -5.S37307E+CC -8.422153E+CC -1.127120E+O1 -1.445323E+01 -1.794084E+01 -2. 171267E+01 -2.575478E+01 -3.006106E+0! -3.463280E+01 -3. S47742E+C1 -4.460383E+O1 -4.999032E+01 -5.536073E+01 -5.887678E+01 -5.083165E+0 1.0 OOOCCE+CO.0000OOOE+00 -3.104096E-01 -1. 231168E+GC -2.737149E+00 -4. 788409E+00 -7. 3257EE+CC -1.030850E+01 - 1. 365154E+C1 -1.730011E+01 -2. 120214E+01 -2. 531988E+01 -2.963211E+01 -3.413350E+01 -3.E8317CE+Cl -4.374248E+01 -4.888018E+1 -5.421573E+01 -5. S417C4E+01. -6.241799E+C1 -5.320110E+CI.COOOOE+00. COOOOOE+00 -4.040566E-01 -1.598828E+00 -3.540497E+00 -6.15730CE+00 -9.152361E+CC -1.30C1001E+C1 -1.700717E+O1 -2, 122856E+01 -2.558191E+01 -3.000868E+01 -3.448766E+0O -3.90 3i170E+01l -4.368009E+C1 -4.848979E+01 -5.352541E+01 -5. e882017E+01 -6.409276E+01 -6.707186E+01 -5.678192E+01.OOOOOOEuCC

.CCOOOOE+OO -5.302218E-01 -2.095428E+00 -4.6281C7E+00 -8.012929E+00 -1.208694E+C1 -1.664859E+O1 -2.147c:EE+O1 -2.637702E+01 -3.118519E+01 -3. 581875E+01 -4.026767E+01 -4.458589E+01 -4.886952E+C1 -5.323277E+Cl -5.778392E+O1 -6.256811E401 -6.720474E+C1 -6. 866632E+01 -5.301947E+01.OOOOCCE+CC.CCOOOCE+00 -6.757254E-01 -2.669671E+C0 -5.890747E+00O -1.017541E+01 -1.528Cl6E+Cl -2.C89230E+01 -2.5666C7E+O! -3.227836E+01 -3. 748813E+01 -4.217918E+01 -4. 637087E+01 -5.019523E+01 -5.385306E+Cl -5.156774E+01 -6.154656E+O0 -6. 591363E+01 -7.028526E+01 -7.103087E+Cl -5.2176C1E+C1.COOOOCE+00.CCCOOOE+00 -8.122385E-C1 -3.208744E+00 -7.C79095E+00 -1. 221475E+Cl -1.8 29152E+C1 -2.4E8261E+O1 -. 5s53EO1 -3.772277E+01 -4.319106E+01 -4. 776980E+O0 -5.152193E+01 -5.468336E+01 -5.759123E+C1 -6.060770E+01 -6.405967E+01 -6.814855E+01 -7. 241537E+01 -7.22q745E+Cl -4.885914E+O1.CCOOCCCE+CC.COOO0CE+OC -9.100826E-01 -3.594322E+00 -7. 2 8624E+00 -1.366842E+01 -2. C42525E4C1 -2.768425E+C1 -?.4866: E+OI -4.142758E+01 -4. 696853E+01 -5. 132882E+01 -5.462038E+01 -5.718602E+01 -5.949959E+C1 -6.2C5999E+01 -6.531188E+O1 -6.952848E+01 -7.411692E+01 -7.336196E+C1 -4.5400C9E+01.00000OE+O00.OOOOOOE+00 -9.494868E-01 -3.745327E+CC -8.270087E+00 -1.425C07E+01 -2.127383E+O1 -2.879453E+01 -3.61S53E+01 -4.291095E+01 -4.850327E+01 -5.280595E+Cl -5.595241E+01 -5.834629E+01 -6.C54203E+01 -6.310914E+01 -6.652848E+01 -7.10C3245E+01 7.563523E+C1 -7. 275680E+C1 -4.C4C076E+01.OOOOOOE+CC.OOOOOCE+OC -9. 207396E-01 -3.638183E+00 -8.031273E+00 -1.38541C0E+01 -2.C72142E+Cl1 -2. 813248E+01 -3.552S89E+1 -4.238735E+01 -4. 831056E+01 -5.312532E+01 -5.692545E+01 -6.005801E+01 -6.301024E+01 -6.626544E+01 -7.C16535E+O0 -7.46811SE+01 -7.7709C7E+01 -6.8956C8E+01 -3.377397E+ 01. OOOOOE+00.CCOOCOE+00 -8.170650E-01 -3. 234513E+CC -7.158698E+00 -1.239840E+C1 -1.865662E+C1 -2.554995E+0 1 -3. 2653C5E+1 ~-3. 95642qE+01 -4. 597961E+01 -5.175258E+C1 -5.691676E+01 -6.165507E+01 -6.621523E+Cl1 -7.077929E+01 -7.526680E+O1 -7.911785E+01 -7.6148C5E+01 -5.921128E+01 -2.523016E+01.000000E400.000000E+00 -6.281894E-01 -2.494376E+00 -5.545207E+00 -9.668041E+00 -I.468567E+CI -2.G36620E+01 -2.644740E+ 01 -3.267398E+01 -3.884043E+01 -4.482520E+01 -5.059757E+01 -5.617299E+01 -6.148169E+01 -6.609819E+01 -6.881643E+O1,,,n -6.733522E+01 -5.882383E+01 -4.C78633E+01 -1.539789E+01.OOOOOOE+CO. COOOOOE+00 -3.485405E-01 -1. 388866E+00 -3.103731E+C0 -5.452514E+00 -8.3674115E+CC -1.115423E+C1 -1.54 735 E+O1 -1.946813E+01 -2.353326E+01 -2.755756E+01 -3. 139438E+01 -3.485750E+01 -3.167026E+01 -3.939451E+O1 -3.935639E+O0 -3.661024E+01 -3.012C61E+01 -1.955S37E+01 -7.081754E+00.OCOOCCE+.CO.CCOOOOE+00.OOOOCOE+00.CCCOCOE+OO.COOOOOE+00.OCOOOCE+CC.CCCOCCE+CC.CCOOCCE+OC.CCOOCE+OO.COOOOE+00.OOOOOOE+00.OOOOOOE+00.COOOOOE+00.OOOOOE+00 OOOOCE+CC.CCOOOOE+00.CCOOCOE+OO.CCOCCCE+00.COOOCCE+00.OOOOOOE+CC. COOOOOE+00.OOOOOCE+OC.COOOOOE+00 -1.883043E-02 -7.547617E-C2 -1.704313E-01 -3.043983E-C1 -4.779258E-C1 -6.SC877aE-01 -'.417187E-01 -1.226607E+00 -1.538326E+00 -1.865124E+00 -2. 189573E+00 -2.487615E+00 -2.727929E+00 -2. El71741E+C00 -2. 873117E00 -2.680200E+00 -2.240558E+00 -1.5253,SE+CC -6.209369E-01.COOOOCE+CC.COOCOOE+00 -3.606675E-02 -1.447103E-01 -3.270738E-01 -5.846578E-01 -S.18605qE-C! -1.328631E+CC -I.Ell5E+CO -2. 359483E+00 -2.957655E+00 -3.582190E+00 -4.197869E+00 -4.756251E+00 -5.194448E+CC -5.434967E+0C.G 87527E-C -4.9s558E+OC -4.054426E+CO -2.669462E+00 -1.027512E+00.00000OE+400.COOOOOE+00 -5.012177E-02 -2.012935E-C1 -4.554421E-01 -8.150496E-C1 -1.282054E+CC -1.856192E+OC -2.532760E+O0 -3.299854E+00 -4.135272E+00 -5.003333E+00 -5.851959E+00 -6.610314E+00 -7.187457E+CO -7.472903E+00 -7.341283E+OO -6.666780E+00 -5.356405E+00 -3.437488E+CC -1.274887E+00.000000E+00.COOOCOE+00 -5. S57515E-02 -2.395523E-01 -5. 428205E-01 -9.731524E-01 -1.533723E+CC -2.224867E+OC -3.C01465E+OC -3.S67291E+00 -4.975048E+00 -6.018517E+00 -7.031312E+00 -7.924196E+00 -8.584228E+CC -8.876560E+O00 -8.652452E+00 -7.772203E+00 -6.1514C7E+CC -3.868283E+00 -1.395964E+CC.0000OCE+CO.OOOOOCE+00 -6.296487E-02 -2.536425E-01 -5.760442E-01 -1.035563E+00 -1.(372CCE+CC -2.382'845E+OC -3.2676C5E+OO -4.275376E+00 -5.374086E+00 -6.511241E+00 -7.610477E+00 -8.569703E+00 -9.261371E+CC -9.536089E+00 -9.233982E+0O -8.21582~E+C0 -6.417412E+00 -3.965267E+CC -1.399005E+00.OOOOOOE+00

.COOOOOE+00 -5.903224E-02 -2.3844e4E-C1 -5.432192E-01 -9.807094E-01 -1.557825E+C0 -2.2191C8E+C. -3.142242E+00 — 4.133248E+00 -5.221071E+00 -6.352386E+00.-7.447550E+00 -8.39885qE+00 -9.072272E+00 -9.314318E+00 -E.'68440E+00. -7.SlC6C5E+00 -6.105229E+00 -3. 14'7E+.CC -1.-286512E+.O.OOOOOOE+00..OOOOOOE+00 -4.731345E-02 -1.918993E-01 -4.3930:1E-Cl -7.975822E-01 -1.275682E+CC -1.EEC92SE+.C -2.6I:]43E+(C -3.470928E+00 -4.424108E+00 -5.429822E+00 -6.416303E,00 -7.282268E+00 -7.899079E+CC -8.120812E+00 -7.805601E0O. -6.852285E+00 -5.24E7S2E+CO -3.162845E+CC -1.08_864E+00.CCOOOCE+CO.CCCOOOE+00 -2.928635E-02 -1.196112E-Cl -2.756771E-01 -5.045935E-01 -8.149866E-Cl -1.215401E+OC -1.711;ICE+OO -2.304694E+00 -2.983648E+00 -3.722815E+00 -4.474210E+00 -5.162557E+00 -5.68421qE+CO -5.915314E+00 -5.734128E+-O -5.C59'9C7E+CO -3.887204E+00 -2.348073E+CO -8.C9l169E-01.00000CE400.CCOOOOE+00 -1. 003438E-02 -4.173250E-02 -9.726093E-02 - 1. 801'322E-01 -2.950050E-Cl -4.472661E-01 -(. 423445E-O1 -8. 846550E-01 -1.175938E+00 -1.512600E+00 -1.881630E+00 -2. 255088E+00 -2.584461E+00 -2.7S8708E+OC -2.813471E+00 -2.561618E+00 -2.026410E+OC -1.264174E+CC -4.560396E-01.COCOOCE+CC.COOCOOE+00.OOOOOOE+00.CCCOOOE+00.CCOOOOOE+00.OOOOO0E+OO.CCOOCCE4CC.CCOOOCE+CC.CCOOOOE+OC.OCOOOOE+00.OOOOOOE+00.OOOOOOE+00.COOOOOE+00.OOOOOOE+00.ooooooE+CC.OOOOOOE+OC.CGOCOOE+OO.000000E+00.COOOOOE+00.OOOOCOE+CC.COOOCOE+CO.OOOOOOE+00 VG(1)...VG(21) 1.884517E-37 4.180647E+01 8.375581E+01 1. 260429E+02 1.687858E+02 2.11452E+C2 2.552647E+02 2.S82064E+02 3.398825E+02 3.789929E+02 4.137776E+C2 4.420003E+02 4.60725EE402 4.676224E+C2 4.585959E+02 4.3C3634E+02 3.793518E+02 3.024777E+02 1.99'7301E+C2 8.398052E+01 1.884517E-37 UL( )...UL(21) \,0 ONx 5.1768C0E-05 4.222789E-05 1.678685E-05 -2.725246E-05 -9.268162E-05 -1.E68194E-C4 -3. 115143E-04 -4.727127E-04 -6.749414E-04 -9.217645E-04 -1.215196E-03 -1.554897E-03 -1.936269E-03 -2.344234E-03 -2.734328E-03 -2.c92454E-O3 -3.C67164E-03 -3. 2758CIE-03 -3.664973E-C3 -4.355241E-03.00000CCE4C HL(1)...HL(21) 1.S25877E-04 2. 701233E-04 2.840120E-04 2.860186E-04 2. 849335E-04 2.828131E-C4 2.8C2910E-04 2.71148CE-04 2.756607E-04 2.772220E-04 2.760033)E-04 2.478657E-04 2.378216E-04 2.107919E-C4 1.75S38CE-C4 1.403507E-04 1.144092E-04 9.C595C5E-05 6.992228E-05 5.33C813E-05 3.424333E-C4. HG(1) o..HG(11) 2.057539E-05 6. 322576E-06 4. 800074E-C6 4.094149E-06 3.589955E-06 3.196971E-Cf 2.881997E-O 2.621537E-06 2.393444E-06 2.250477E-06 2.493679E-06 NU'L = 44.99614C' NUG. 4.2954C4, RMAX. 9.00633EE-C39 Rl l I.2(:1686 NE1 = 2, NE2 6, RAL = 2.859151E+C8, RAG= 2.870898E +05 CL =.346032, CG=.185567, RP = 1.015559, RICIC =.261686 SFI 1)...SFI(21) 1.884517E-3.7 1.916388E-04. 7.665551.e-04 1.724749)E-03 3.066220E-03 4.790969E-03 6.898996E-03 9.390300E-03 1.226488E-02 1.552274E-02 1.916388E-02 2.318829E-02 2.759598)E-02 3.238695'E-C2 3.756120E-02 4.311872)-02 4.905S53E-02 5.538361E-02 6.205096E-02 6.918160E —02 1.884517)E-37 DQILDT = -1.127238E-0'3

REFERENCES L. Clark,, J. A. A Review of Pressurization, Stratification and Interfaciai Phenomena: International Advances in Cryogenic Engineering, Vol. 10,~ pp. 259-284, 1966. 2. Thomas., P.D. and. Morse,, P.H. Analytical Solution for the Phase Change in a Suddenly Pressurized Liquid-Vapor System: Advances in Cryogenic Enginee~ring., Vol. 8, pp. 550-562, 196)4. 5. Knuth., E.L. Evaporation and Condensation of One Component Systems: Journal of the American Rocket Society: Vol. 32, pp. 1)42)4-1426,~ 1.962. 4. Yang, W.J. Phase Change of One-Component System'in a Container: AIChE Paper No. 63-A)48, presented at the 6th National Heat Transfer Conference., Boston, Mass., August 1965. 5. Yang, W.J. and Clark, J.A. On the Application of the Source Theory, on Problems Involving Phase Change-Part 2: J. Heat Transfer., Trans. ASME, Ser. C. Vol. 86, pp. 4)45, 196)4. 6. Yang, W. J., Larsen., P. S. and Clark, J. A. Interfacial Heat and Mass Trans*fer in a Suddenly Pressurized,, Binary Liquild-Vapor System: Trans. ASME., J. Eng. for Ind.,. Vol. 8y, pp. 415, 1965. 7. O'Loughlin., J.R. and Glenn., H. Bulk Liquid Mass Transfer with Variable Ullage Pressure: Proceedings of the Conference on PropellantU Tank Pressurization and Stratification, NASA MSFC, Jan. 20-21,, Vol. 1, 1.965. 8. Epstein,, M.,, Georgius,, H. K. and Anderson, R. E. Generalized Propellant Tank Pressurization Analysis: International Advances in Cryogeniic Engineering, Vol. 10, 1966. 9. Huntley,, S.C. Temperature-Pressure-Time Relations in a Closed Cryogenic. Container: Advances in Cryogenic Engineering., Vol. 5, p. 54-2, 1,96'0.10. Leiben-berg,, D.H. and Edescuty,F.J. Pressurization Analysis o-f a Large Scale Hydrogen Dewar: International Advances in Cryogenic Engineering, Vol. 10, 1966. 11. Van Wylen., G.J. Thermodynamics: John Wiley and Sons, Inc. pp. 425, 1960.

REFERENCES (Concluded) 15. Jacob, M., Heat Transfer, Vol. I: John Wiley and Sons,, Inc. Chap. 29. l14. Swim,, R. T. Temperature Distribution in Liquid and Vapor Phases of Helium in Cylindrical Dewars: Advances in Cryogenic Engineering., Vol. 5, 196o. l5. Brenteri,. E. G., Giarratano,. P. J. and Smith,, R. V.: Boiling Heat Transfer for Oxygen., Nitrogen', Hydrogen and Helium: NBS Technical Note No. 517, Sept. 1965. 16. Barakat, H. Z. and Clark,, J. A. Transient Natural Convection Flows in Closed Containers: Technical Report O4268-lO-T., Office of Research Administration., University of Michigan., Ann Arbor., Michigan, Contract NAS 8-825., Marshall Space Flight Center,, Aug. 1965.