ARL 65-58 A DIGITAL COMPUTER PROGRAM FOR CONDENSATION IN EXPANDING ONE-COMPONENT FLOWS LEONARD J. HARDING UNIVERSITY OF MICHIGAN ANN ARBOR, MICHIGAN MARCH 1965 Contract AF 33(657)-8867 Project 7116 Task 7116-01 AEROSPACE RESEARCH LABORATORIES OFFICE OF AEROSPACE RESEARCH UNITED STATES AIR FORCE WRIGHT-PATTERSON AIR FORCE BASE, OHIO

FOREWORD This technical report was prepared by the College of Engineering, The University of Michigan, Ann Arbor, Michigan, on Contract AF 33(657)8867 for the Aerospace Research Laboratories, Office of Aerospace Research, United States Air Force. The research reported herein was accomplished on Task 7116-01, "Internal Flow Research, " of Project 7116, "Energy Conversion Research," under the technical cognizance of Everett D. Stephens of the Thermo-Mechanics Research Laboratory of the Aerospace Research Laboratories. This report represents a continuation of the work, "Digital Computer Analysis of Condensation in Highly Expanded Flows, " by James L. Griffin reported in ARL 63-206, November 1963.

ABSTRACT This report describes a digital computer program for calculating vapor condensation processes that occur in rapidly expanding flows, The treatment emphasizes the program logic required to satisfy the requirements of the mathematical model Among the most important of these requirements are the search for the onset of nucleation, in an isentropically expanding flow, and the iteration procedure necessary for the joint solution of the nucleation and growth equations and the diabatic flow equations. The program features flexibility in allowable input conditions, short execution time, and a convenient output format. 111

TABLE OF CONTENTS Page 1o INTRODUCTION 1 2, MATHEMATICAL MODEL 2 2. 1 NOMENCLATURE 2 2, 2 ISENTROPIC EXPANSION 5 2. 3 NOZZLE GEOMETRY 6 2. 4 VAPOR SATURATION DATA 6 2. 5 LATENT HEAT AND LIQUID DENSITY DATA 7 2. 6 SURFACE TENSION DATA AND CRITICAL DROP SIZE 7 2. 7 SATURATION POINT DETERMINATION 8 2. 8 ONSET OF CONDENSATION DETERMINATION 8 2. 9 CONDENSING FLOW 10 2.10 ISEN SUBROUTINE CALCULATIONS 12 3. PROGRAM NOMENCLATURE 15 3,1 PROGRAM COMMON VARIABLES, GROUP 1 15 3, 2 PROGRAM COMMON VARIABLES, GROUP 2 16 30 3 MAIN PROGRAM INPUT PARAMETERS 17 3, 4 SUBROUTINE VAPOR INPUT PARAMETERS 17 3. 5 SUBROUTINE NOZZLE INPUT PARAMETERS 18 30. 6 SUBROUTINE CONDEN INPUT PARAMETERS 19 3, 7 SUBROUTINE ISEN OUTPUT VARIABLES 19 4, PROGRAM STRUCTURE 21 4. 1 PROGRAM INPUT AND INITIALIZATION 21 4. 2 PROBLEM INPUT AND INITIALIZATION 22 4. 3 SATURATION POINT COMPUTATION 22 4, 4 CALCULATION OF THE INITIAL POINT 23 4. 5 ITERATIVE CALCULATION OF THE FLOW CONDITIONS 24 4. 6 PROBLEM TERMINATION 25 50 STRUCTURE OF THE SUBROUTINES 27 6. I/O STRUCTURE 37 6,1 INPUT STRUCTURE 37 6. 2 OUTPUT STRUCTURE 41 7, REFERENCES 44 APPENDICES A, PROGRAM LISTING 45 Bo EXAMPLE INPUT 67 C. EXAMPLE OUTPUT 69 iv

1. INTRODUCTION In Ref. 1, Griffin reported on a digital computer program for analyzing condensation processes occurring in rapidly expanding flows of pure vapors. The present report describes, from the computer programmer's point-ofview, the current version of that original program. Additional details on the basis for changes made in the mathematical model and on results obtained with the current program are presented in Ref. 2. As was the original program, this version is coded in the Michigan Algorithm Decoder (MAD) language, Ref. 3, for running on an IBM 7090 with The University of Michigan Executive System. (Information on the Michigan Executive System is available through the SHARE organization.) In actuality the program herein described represents, except for a few subroutines, a recording of the original program that allows considerably more flexibility and at the same time a substantial decrease in computer usage time per problem. Although the mathematical model has not been altered significantly, the logical structure has been changed considerably. The program has been tested and extensively use; however, no guarantee as to its accuracy or functioning can be made. Manuscript released by the author September 1964 for publication as an ARL Technical Documentary Report.

2. MATHEMATICAL MODEL (by Ko R Sivier) The following listing of equations, along with brief comments about their use and the interrelationships between them, forms the basis of the mathematical model for the computer program and is included for completeness. The treatment has been kept brief since the mathematical model is discussed in detail elsewhere. Reference 1 presents a detailed discussion of the mathematical model of the original program. Most of this material continues to apply to the present modified program. Reference 2 presents discussions of Tolman surface tension correction, the change in the method of dealing with the vapor saturation data, and the equations involved in the ISEN subroutine. The foreign nuclei and variable step modifications involve program logic and operation only and, hence, require no new equations. They are, however, also discussed in Reference 20 Much of the following material has been taken directly from Section V-B of Refo 1 with little or no change. The treatment there is sufficiently concise to recommend its inclusion here. 2. 1 NOMENCLATURE The following list contains only those symbols used in the equations in the present description of the mathematical model. To be consistent with the computer program itself, all units are in the cgs system. A Nozzle cross-section area (cm 2) A* Nozzle throat area (cm ) C Vapor specific heat at constant pressure ~~P ~(dyne-cm/gm- K) D Tolman constant used to correct surface tension for finite drop size (cm) g Mass fraction of the mixture that is in the condensed phase 2

J Specific drop nucleation rate (drops/sec-cm3) k BoltzmannTs gas constant (1. 379 x 10-16 dyne-cm/0K) L Latent heat of vaporization (dyne-cm/gm) (L) o o o (L)1 Empirical latent heat constants, see Eq. (12) M Mach number m0 Rate of mass flow (gm/sec) N Drop nucleation rate (drops/sec) NA Avogadro's number (6. 027 x 1023 molecules/gmol) p Free-stream static pressure (dyne/cm2) p Supply (stagnation) pressure (dyne/cm2) p0 Stagnation pressure behind a normal shock wave (dyne/cm2) PCO Saturation vapor pressure for a plane surface of liquid (dyne/cm2) (PSAT) o.a Empirical saturation curve constants, (PSAT)4 see Eqs. (10) and (11) R Universal gas constant (8. 314 x 10 dyne-cm/gmol-~K) (RHOL)... Empirical liquid density constants, (RHOL)1 see Eq. (13) r Drop radius (cm) r* Critical drop radius (cm) (SIGMA)... Empirical surface tension constants, (SIGMA)1 see Eq. (14) T Free-stream static temperature (OK) T Supply (stagnation) temperature (OK) T Temperature of a drop of condensate (OK) drop T Saturation vapor temperature for a plane 00a~ ~surface of liquid (OK) 3

U Velocity (cm/sec) x Nozzle station measured from the throat (cm) a Accomodation coefficient 2' Ratio of specific heats e Fractional deviation, in a distance Ax, of the static pressure change from the corresponding change for an isentropic expansion ET Fractional deviation, in a distance Ax, of the static temperature change from the corresponding change for an isentropic expansion Ecritical fractional deviation, in a distance Ax, of the change in static pressure or temperature, from the corresponding isentropic change, that occurs at the condensation onset point Molecular weight (gm/gmol) p Free-stream static density of the mixture (gm/cm3) p Supply (stagnation) density (gm/cm3) 3 PL Liquid phase density (gm/cm ) a Surface tension (dyne/cm) 9a Surface tension of a plane liquid surface (dyne/cm) i. Half-angle of the convergent section of the nozzle (deg) 0 out Half-angle of the divergent section of the nozzle (deg) Subscripts ( ) Refers to conditions existing at the condensation onset point as established by the program ( ). Refers to the incremental step in which the droplets were originally formed 4

Refers to the incremental step presently being evaluated ( ) Refers to conditions existing in an isentropically expanding flow ()sat Refers to conditions existing at the vapor saturation point Refers to an incremental step Ax, at temperature T, that is being tested in the search for the onset of condensation 2. 2 ISENTROPIC EXPANSION The expansion of the vapor from its supply condition to the saturation point and, if the condensation onset point is found, the expansion of the supersaturated vapor from the saturation point to the condensation onset point are governed by the following isentropic flow equations: PP = Po (1) RMSIr r~i2 0tT~- ~ 1(2) 1 1, Y M (3) Y+l A-A [ 1 1+ 1 M2 (4) U= [2C (T -T)]2 (5) For a given value of T, the above equations are readily evaluated for given values of p0, T0, p0, y, Cp% and A*,

2o 3 NOZZLE GEOMETRY In the case of a wedge nozzle (two-dimensional flow), the nozzle is assumed to have unit width and its area is given by A = A* - 2x tan 0in (6) upstream (x < 0) of the throat (x - 0), and A = A* + 2x tan 0out (7) downstream (x > 0) of the throat, For the conical nozzle (three-dimensional flow), the area is given by A=ir ( - x tan 0i) (8) upstream of the throat, and A = 7T ( + x tan sot) (9) downstream of the throat. 20 4 VAPOR SATURATION DATA Vapor saturation data are supplied to the program in the form of curve fits to empirical data, expressed in terms of In po and 1/To Below the triple point of the vapor, the linear (Clausius-Clapeyron) approximation 1 = (PSAT) + (PSAT)1 In p, (10) 0) is used and, above the triple point, the slight non-linearity in the data is accounted for by the quadratic approximation __ = (PSAT)2 + (PSAT) In Poo + (PSAT)4 (In Poo) (11) 09

where (PSAT) o 0o (PSAT)4 are empirically determined constants. 2. 5 LATENT HEAT AND LIQUID DENSITY DATA Latent heat of vaporization and liquid density data are supplied to the program as linear functions of temperature; L = (L) + (L)1 T (12) and PL = (RHOL) + (RHOL)1 Tdrop 9 (13) where (L), (L)1, (RHOL), and (RHOL)1 are empirically determined constants and the temperature T is taken as the saturation temperadrop ture corresponding to the local ambient pressure of the vapor0 20 6 SURFACE TENSION DATA AND CRITICAL DROP SIZE Surface tension data, corresponding to an infinite plane liquid surface, is supplied to the program as a linear function of temperature (an E'otvos-type variation), io eo, v - (SIGMA) + (SIGMA)1 Tdrop s (14) where T is taken as the saturation temperature corresponding to drop the local vapor pressure, An estimate of the surface tension for very small drops is obtained from this value by use of the Tolman relation c= D 9 (15) r where D is taken as an empirically evaluated constant0 As shown by Eq. (29) below the critical drop size is a linear function of surface tension, iL e, 7

r* ( 2,4: PL RT in p Pl The joint solution of this relation with Eq. (15) for surface tension results in the following simple expression for critical drop size with the corrected surface tension value: r* = 2 ao - D (16) PL RT In P 2, 7 SATURATION POINT DETERMINATION The state of the vapor at its saturation point is established by the joint solution of the appropriate saturation data approximation, Eqo (10) or (11), and the isentropic relation y sasat~pot~! 0(17) Psat = P T! The values of psat and Tsat' thus found, are used with the isentropic relations, Eqso (2) through (5), to calculate Msat A sat'P p sat and UsatO The location of the saturation point, Xsat' is determined from the value of Asat and the appropriate nozzle equation, one of eqs. (6) through (9). The total mass flow through the nozzle is calculated from the equation m - Psat sat Asat 0 (18) 2. 8 ONSET OF CONDENSATION DETERMINATION The point of onset of condensation is determined by taking incremental temperature steps, AT, of supersaturation and testing for 8

significant condensation effects at the end of each successive step. To initiate this process, the isentropic values for the state of the vapor at the end of each AT step are introduced into the nucleation equations for critical drop size and formation rate; i. eo, 2I0 -D (19) PL RT In p- - ( 4X,r *2 P2 1 f ei 3kT T = kT pL e zN 0 (20) Note that the choice of a very small AT will cause the first step to remain very near the saturation point and that, for p/po near unity, r* approaches infinity. However, r* falls off very rapidly as supersaturation increases and for most cases a choice of AT > 50K is sufficient to avoid this problem. Since it is assumed that the critical drops experience no growth in the Ax increment in which they are formed and since the nucleation rate for this increment is NT = JTT AT XT (21) the mass fraction of condensate formed in this increment is 4PL 0 3 AgT= 3m NT rT (22) Next, the epsilon equations A - 1 AgT (23) P

A M L; CUT 1AgT (24) which test the effect of condensation on static pressure and temperature, respectively, are used to determine if this amount of condensate is sufficient to cause a significant deviation from the isentropic expansion. Because of their relative sensitivities, E is used for subsonic Mach p numbers and eT is used for supersonic Mach numbers, The value of ET, or p is compared with the value of Ecitical) which is assumed to represent a significant condensation effect and which is supplied as part of the initial data for the program. If the eT is less than critical the calculations are repeated at step T + ATo critical This process continues until ecritical is equalled or exceeded. When exceeded, a bracketing procedure is employed to improve the estimate. All quantities computed for Ag at values of e < Ecritical are discarded and it is assumed that there has been no condensate formed prior to e e.iti.cal Once the condensation onset temperature ound, the -critical' isentropic equations, Eqs. (1) through (5), and the appropriate nozzle equation, Eqso (6) through (9), are used to determine the following values at the condensation point: p XM A U and x con conP con' con' con' con 2. 9 CONDENSING FLOW The condensing portion of the flow requires a joint solution of the nucleation and growth equations and the diabatic flow equation. The calculations are performed for increments of Ax, starting from the onset of condensation~ The diabatic flow equations, in incremental form, are 10

AP + + 0 (Continuity) (25) p U A Ap _ UAU -p (1- R(Momentum) (26) P R P (1 - g) T Ap Ap AT Ag p p +T g) (State) (27) UAU + C AT - LAg = 0 (Energy) (28) The value of Ag, appearing in the above equations, is determined from the following nucleation and growth equations: r* =.- (29) )2 1(21/2 e3kT (30) N. =J. A. AX. (31) 1 11 1 Ar. P 2 1/2' (T T) x (32) (P) drop N LP( r. r(. - + Ar (33) ~~1)/2~~~~~~ Agj= N r A r + N r*j (34) 1 ij (35 g = gj-l A r j (3+ )

In the above, the subscripts i and j are numbered increments of x starting from the condensation point, i acts as a label for each group of droplets of a particular size and denotes the particular increment in which these drops originated as critical sized drops~ j denotes the increments presently under consideration. Thus a designation r13 denotes droplets which were formed as critical drops in increment 1, have undergone the growth Ar2, and are now undergoing the growth Ar3" 2. 10 ISEN SUBROUTINE CALCULATIONS The ISEN subroutine is used to evaluate the changes in the local flow conditions that have been produced by the condensation processes, To do this, the conditions computed for the condensing flow are compared with the corresponding isentropic flow conditions that would exist if condensation was rhot occurring. The isentropic flow conditions are based entirely on the geometric area ratio. An inversion of Eqo (4) is performed to determine the isentropic Mach number, Ms, corresponding to the local nozzle area. M is determined by a half-interval iteration technique, starting with the following first approximation for Ms; _____ Y+ 1) Ms = l' + t- Y* ~ 1) ( 1 + (36) Once M is determined, the isentropic flow equations yield the following isentropic ratios for static temperature, pressure, and density: 12

TsS = 1 + 2 M1 - (37) 0 PT S -'= 1+ 2.... P =,= 1+ 2i M (38) Ps PYo i 2 sY PS P =(1+ - 2 M (39) In addition, Ms is used to calculate the stagnation pressure ratio across a normal shock wave for the isentropically expanded flow; iL eo, 7' - 1 P (y+ 1) __ +s Y+J(40) The static ratios for the condensing flow are found simply from the relations T5 1 (41) 0 p P (42) p0 P P (43) In addition, the Mach number of the condensing flow is used to obtain the ratio of stagnation pressure downstream of a normal shock wave to the upstream static pressure (completely neglecting the effect of the condensate on the recovery process); i0 e, 13r

P L — i) Ir2 2'] 0 (44) ~ 2 ([y+ 1 ](4 Finally, the ratios indicating the extent of the non-isentropic condensation effects are found directly as: AT T T= T (45) T Ts /~ P = P (46) P P A p p p p5 (47) PS PS and ~ -{Do - ~ ____ L (48)

3, PROGRAM NOMENCLATURE This partial dictionary of symbols used in the main program and subroutines has been divided into seven groups. Two of these groupings are based on the universality of definition (special storage allocation) of the variables; four more are input lists for different parts of the program and the remaining group is a partial list of output designators. 301 PROGRAM COMMON VARIABLES, GROUP 1 These symbols have the meanings given below in the main program and the subroutines IFLOW, CFLOW, NUCLE1I, NUCLE2, CONDEN, ISEN and NOZZLEo ASTAR Throat area (cm ) PZERO Supply pressure (dyne/cm2) TZERO Supply temperature (OK) RHZERO Supply density (gm/cm3) GAMMA Ratio of specific heats MU Molecular weight (gm/gmol) CP Specific heat at constant pressure (dyne *cm/gm *0K) L Current value of the latent heat (dyne*cm/gm) SIGMA Current value of the surface tension (dyne/cm) RHOL Current value of the liquid density (gm/cm3) D Intermolecular distance in the liquid (cm) ALPHA Accomodation coefficient N Current value of the index. The main body of the program is concerned with computing the conditions at station N, based on the conditions at stations 0,. o., N-1. 15

3, 2 PROGRAM COMMON VARIABLES, GROUP 2 All of these variables arei Vectors of length 500. Except for RADIUS and DELRAD, the N-th entries of these vectors define the flow conditions at the N-th station. These symbols have the meanings given below in the main program and the subroutines CFLOW, NUCLEI, NUCLE2, and CONDEN. X(I) Nozzle coordinate at station I (cm) DELX(I) By definition DELX(I) = X(I) - X(I-1)0 The initial entry DELX(O) = DELX is an input parameter0 A(I) Area at station I (cm2) P(I) Pressure at station I (dyne/cm2) T(I) Temperature at station I (OK) RHO(I) Density at station I (gm/cm 3) U(I) Velocity at station I (cm/sec) M(I) Mach number corresponding to U(I) G(I) Mass fraction of condensate at station I DELG(I) Increase in mass fraction of condensate between stations I-1 and I NDOT(I) Number of drops per second 3rmed between stations I-1 and I RADIUS When computing the conditions at station N. RADIUS(I) is the radius, at station N-i, of the drops initially formed at station L Clearly, only entries 0, 0 0 *, N-1 are thus defined when computing the conditions at station N; the N-th entry is computed at station N (cm) DELRAD When computing the conditions at station N, DELRAD(I) is the increase in radius of the drops, initially formed at station I, since passing station N-1. When the conditions at station N have been computed, RADIUS(I) must thus be incremented by DELRAD(I) for I 0 o o0.09 N-lo 16

AI,::. IAE K 3.3 MAIN PROGRAM INPUT PARAMETERS These symbols have the meanings given below only in the main program. They do not represent the entire list of input parameters to the main program, but all of them are input variables0 PE Input parameter used as an initial approxi- 2 mation to the saturation pressure (dyne/cm ) XRANGE Length of the interval, starting at the initial point, in which the flow conditions are to be computed (cm) MAXG If the mass fraction of condensate exceeds this value, the calculations are halted MING When the mass fraction of condensate exceeds this value, then DELX(N) will be varied so that DELG(N) is greater than MNDELG but less than MXDELG MNDELG See MING MXDELG See MING MNDELX Minimum value that DELX(N) will be given, regardless of the computed DELG(N) (cm) MXDELX Maximum value that DELX(N) will be given, regardless of the computed DELG(N). This bound only applies when NDOT(N) is zero (cm) MXDELX(1) Maximum value that DELX(N) will be given, regardless of tie computed DELG(N)o This bound only applies when NDOT(N) is nonzero (cm) XPOINT Only those stations with subscripts which are integral multiples of XPOINT will be printed in the three output tables. 3o 4 SUBROUTINE VAPOR INPUT PARAMETERS These symbols have the meanings given below only in the subroutine VAPOR0 Three of these names are in the PROGRAM COMMON GROUP 1; however, VAPOR has no communication with PROGRAM COMMON and 17

the meanings are different. These variables must be supplied to VAPOR, as they provide the basis for the entries FSIGMA, FRHOL, FL, FTSAT and LNPSATo PTP The natural log of the pressure in dyne/cm2 at the triple point. TTP The inverse of the temperature in OK at the triple point0 SIGMA o o 0 Coefficients of the linear approximation to SIGMA (1) the surface tension in terms of the drop temperature, SIGMA + SIGMA(1)*TDRoP. RHOL o o o Coefficients of the linear approximation to RHOL(1) the liquid density in terms of the drop temperature, RHOL + RHOL(1)*TDRoP. L o o Coefficients of the linear approximation to L(1) the latent heat in terms of the vapor temperature, L + L(1)*TVAPOR. PSAT. o Coefficients of the linear approximation to PSAT(4) the saturation curve below the triple point, 1/T = PSAT + PSAT(1)*ln(P) and of the quadratic approximation to it above the triple point 1/T = PSAT(2) + PSAT(3)*ln(P) + PSAT(4)*[ln(P)]2 30, 5 SUBROUTINE NOZZLE INPUT PARAMETERS These five symbols have the meanings given below only in NOZZLE. These variables, together with ASTAR, are read by this subroutine and serve to define the nozzle geometry0 XMIN X co-ordinate of the nozzle inlet (cm) XMAX X co-ordinate of the nozzle exit (cm) INANG Half-angle of the convergent section of the nozzle (Deg) OUTANG Half-angle of the divergent section of the nozzle (Deg) 18

WEDGE Switch denoting whether the geometry is for a wedge nozzle (1B) or for a conical nozzle (OB)o This variable is of Boolean mode and may assume only the values lB,o true or OB ^o false. 3. 6 SUBROUTINE CONDEN INPUT PARAMETERS The symbols given below are defined only in CONDENo The input parameter EPSLON must always be provided. If EPSLON is zero, then any PROGRAM COMMON symbol is a legal input parameter, If EPSLON is non-zero, then a condensation point search employing the E(T) criterion is initiated; in this case, these input parameters have the meanings given below. DELT(1) Initial step in T from saturation point that will be used in searching for a maximum DELG point (OK) TRANGE The search for the maximum DELG point will be confined to the interval TSAT to TSAT - TRANGE. If TRANGE exceeds TSAT the lower end of the interval is taken as zero (OK) EPS Condensation is assumed to be started if E(T) exceeds EPS at some point EPS(1) Upper limit of the absolute error in the temperature at the maximum DELG point. EPS(2) Upper limit for the relative error of the temperature at the condensation point. 3, 7 SUBROUTINE ISEN OUTPUT VARIABLES Except for the symbol MI, these names are not defined in any part of the program. They are, rather, headings for some of the columns of Table 1. All are computed by the subroutine ISEN, and are returned to the main program via PROGRAM COMMON, 19

PHAT Ratio of the computed static pressure to the corresponding isentropic pressure, THAT Ratio of the computed static temperature to the corresponding isentropic temperatureo RHOHAT Ratio of the computed static density to the corresponding isentropic density~ PZEROtHAT Ratio of the computed total head pressure behind a normal shock wave to the corresponding value for an isentropic flow, MI Isentropic Mach number corresponding to the local geometric area ratio~ 20

4, PROGRAM STRUCTURE Before beginning the general description of the logical structure of the program, a few preliminary remarks are in order. To provide easy communication between the main program and some of the subroutines, many of the variables have been placed in PROGRAM COMMON; this storage assignment function is conceptually the same as FORTRAN COMMONo The function of the'READ DATA" statement referred to below is explained in the section on the data deck. Because the development of the program has been essentially experimental in nature, no attempt has been made to provide flow charts. The program has been so organized, however, that the particulars of the various iterations used in the calculations are easily obtainable from the program listing in Appendix A. To a large extent, this is due to the fact that the main program plays the role of an upgraded I/O monitor for a system of computational subroutines. Thus, this section, though intended to be a description only of the main program, makes repeated references to the subroutines and their functions0 401 PROGRAM INPUT AND INITIALIZATION The first action of the main program, hereafter referred to as MAIN, is the reading of the program initialization deck. This deck begins with a set of 72-column comment cards which are continuously read and printed until a card with ENDbbb(1) in columns 1-6 is encountered. Immediately following this, MAIN calls VAPOR which gives a?READ DATA" statement to obtain the SUBROUTINE VAPOR INPUT PARAMETERS. It should be (1)A small T'b represents a blank. 21

noted that these parameters uniquely specify the vapor and, as a consequence, all problems in a given run are based on the expansion of a single pure vapor. This input represents the end of the program initialization deck, and MAIN proceeds to the processing of the first problem of the run0 4. 2 PROBLEM INPUT AND INITIALIZATION The data deck for each problem also begins with a set of 72-column comment cards, in this case terminated by a card with DATAbb in columns 1-6. Following this comment card processing, MAIN gives a'READ DATA' statement to obtain the PROGRAM COMMON variables, TZERO, PZERO, GAMMA, MU, CP, D and ALPHA, and the entire list of MAIN PROGRAM INPUT PARAMETERSo After printing these input data, MAIN calls NOZZLE with a first argument of $THROAT$(2) Because of this value for its first argument, this subroutine gives a'READ DATA' to obtain the PROGRAM COMMON variable ASTAR and the SUBROUTINE NOZZLE INPUT PARAMETERS0 NOZZLE then prints these data and does some initialization for its computational entries which are specified by first arguments of $AREA$ or $INVERS$O 4, 3 SATURATION POINT COMPUTATION With these data, MAIN proceeds to the calculation and printing of the conditions at the saturation point. The saturation temperature is computed by performing an iteration on temperature and pressure that leads to the point of intersection of the saturation curve specified by FTSAT and the isentropic expansion curve given by (T ZRO GAMMA/(GAMMA-1) P PZ TZEROTZERO (2)Hollerith arguments to MAD subroutines may be given directly in this form, e. go, NOZZLE. ($THROAT$, X, A0 ) 22

The conditions at the saturation point are then supplied by IFLOW and NOZZLE ($INVERS$); IFLOW provides the pressure, density, area, velocity and Mach number corresponding to an isentropic expansion to the saturation temperature, and NOZZLE provides the coordinate corresponding to this area and Mach number. The saturation temperature and rate of mass flow, MDOT = ASAT * RHO AT* U then serve SAT SAT SAT' as the arguments to CONDEN. 4, 4 CALCULATION OF THE INITIAL POINT The initial point X z X(O) and the associated flow conditions A, P, T, RHO, U, M, RADIUS, NDOT, G and DELG are supplied to MAIN by the subroutine CONDEN. An explanation of how these initial conditions are computed is given in the description of CONDEN and will not be repeated here. It should be noted, however, that CONDEN gives a'READ DATAt in order to decide how to compute these initial conditions, prints these conditions, and may give an error return to MAIN if it finds linsufficient' nucleation. If this error return is given, MAIN immediately proceeds to the next problem. Having executed CONDEN with a successful returns MAIN proceeds with the calculation of the flow conditions at stations 1, 2, 3. o o The flow conditions to be computed at station N consist of: X(N), DELX(N), A(N), P(N), T(N), RHO(N), U(N), M(N), NDOT(N), G(N), DELG(N) and RADIUS.. RADIUS(N). In passing from the conditions at station N-1 to those at station N, the vector elements RADIUS... RADIUS(N-1) are incremented by the growth these drops experience between these stations and RADIUS(N) becomes defined. Thus the entire set of conditions at station N-1 are no longer available. 23

4. 5 ITERATIVE CALCULATION OF THE FLOW CONDITIONS Given the flow conditions at station N-1 and taking DELX(N) = DELX(N-1), the computation of the conditions at X(N) = X(N-1) + DELX(N) is accomplished via an iteration on the value of G(N). Starting with DELG(N-1) as an initial approximation to DELG(N), CFLOW is used to solve the diabatic flow equations for the pressure, temperature, density and velocity corresponding to X(N) and C(N)o The initial approximation to U(N) required by CFLOW is 2 * U(N-1) - U(N-2) on the first iteration; thereafter it is supplied automatically by the previous iteration. These values of temperature, pressure and velocity are then used by NUCLE2 to calculate the corresponding DELG(N) and, hence, a new G(N). This alternating use of CFLOW and NUCLE2 is continued until two successive values of G(N) differ by no more than 0 00001 or until ten iterations have been performed. If this iteration fails to converge, then DELX(N) is halved and the computations are restarted0 Consistent failure to converge results in termination of the problem; however, only rarely has this iteration failed and then convergence was obtained after halving DELX(N) only once. Assuming the iteration has converged, the values of X(N), DELX(N), 0. o are accepted as specifying the conditions at station N if any one of the following is true. 1) It was necessary to halve DELX(N) in order to obtain convergence. 2) G(N) does not exceed MING. 3) DELG(N) is greater than MNDELG but less than MXDELGo 4) DELG(N) is less than MNDELG but 2 * DELX(N) exceeds MXDELX or MXDELX(1), whichever applies. 5) DELG(N) exceeds MXDELG but o 5 * DELX(N) is less than MNDELXo 24

The decision to accept or reject DELX(N) as specifying the nozzle coordinate at station N is reached by testing these conditions in the order given. If DELX(N) is rejected, a new value is obtained by halving it when DELG(N) exceeds MXDELG and doubling it when DELG(N) is less than MNDELGO Note that conditions (4) and (5) insure that this new DELX(N) lies in the range specified for it. Having thus obtained a new value for DELX(N), MAIN returns to compute the corresponding flow conditions, Consideration of these five conditions shows that a slightly more sophisticated technique is required, since currently the program can enter an infinite loop by alternately halving and doubling DELX(N). The essential point is that, at present, the program does not remember whether or why it rejected the last value of DELX(N)O If a procedure for remembering this information was incorporated, the algorithm would be improved and effective. To date there are no positive indications that this loop has been encountered in our use of the program; however, its occurrence could be forced with proper values for the parameters MNDELG, MXDELG, MNDELX, MXDELX and MXDELX(1)O When MAIN accepts X(N), DELX(N), 0 as specifying the conditions at station N, the internal subroutine IOCTRL is executed, If N is not an integral multiple of XPOINT, this subroutine returns without taking any action. Otherwise, IOCTRL calls on ISEN to compute the SUBROUTINE ISEN OUTPUT VARIABLES, prints the Table 1 output line for station N, computes MDRAD(N) and RHODRP(N) and saves RADIUS(N) in SVRAD(N). These last three quantities are printed in Table 3. 4. 6 PROBLEM TERMINATION Termination of these calculations occurs when any one of the following conditions is found to hold. 1) N exceeds 500, 2) The next point at which the flow conditions are to be calculated exceeds X(0) + XRANGEo 25

3) The subroutine NOZZLE ($.AREA$) indicates that the next point at which the flow conditions are to be calculated is outside the nozzle. 4) The iteration on the diabatic flow equations in CFLOW has failed to converge. 5) The iteration to balance the diabatic flow and nucleation equations has failed to converge0 6) The mass fraction of condensate exceeds MAXGo The cause of termination and the output Table 1 line for this last station are printed0 Finally, Tables 2 and 3 are printed and MAIN continues to the next problem. The output lines, in Tables 2 and 3, for the last station are not printed unless the index of the last point is a multiple of XPOINT. The contents of all three output tables are fully described in the section on I/O structure0 26

5. STRUCTURE OF THE SUBROUTINES The following sections provide an explanation of both implicit (PROGRAM COMMON) and explicit arguments, length, transfer vectors and functions of the subroutines coded for this program. Subroutines dealing with I/O are omitted from the transfer vectors. Of those system subroutines listed, all but ZERO are fully explained by their names; this subroutine simply sets its arguments to zero. The formulas underlying the computations performed by CFLOW, CONDEN, NUCLEi and NUCLE2 have been incorporated, while those for the remainder may be easily obtained from the program listing, CFLOW Arguments: PROGRAM COMMON VARIABLES, Group 1o PROGRAM COMMON VARIABLES, Group 2o DELA Area increment corresponding to N, A(N) - A(N - 1). LOC Location'to be given control if the iteration fails to converge. Length: 277 octal Transfer Vector: FL, SQRT Purpose: Compute values of P(N), T(N), RHO(N), U(N) and M(N) that satisfy the diabatic flow equations DELRHO DELU DELA RHO(N) + U(N) + A(N) (Continuity) DELP MU * U(N) * DELU P(N) = -(1 - G(N)) * R * T(N) (Momentum) DELP DELRHO DELT DELG(N) P(N) RHO(N) + T(N) -1 - G(N) (State) U(N) * DELU + CP * DELT - L(T(N)) * DELG(N), 0, (Energy) 27

for the given values of A(N), DELA, G(N) and DELG(N)o The symbol R is the universal gas constant. The symbols DELRHO, DELU,. O O denote the differences RHO(N) - RHO(N-1), U(N) - U(N-1), 0 0. The constants MU and CP are obtained from PROGRAM COMMONO CFLOW uses FL to compute the latent heat L, and leaves the final value corresponding to T(N) in PROGRAM COMMON. An initial approximation to U(N) must be supplied in U(N)o Iteration on the value of U(N) is accomplished by solving the energy equation for DELT and, hence, T(N) and then using the equations of continuity, state and momentum to compute the corresponding DELUo The iteration is continued until the relative error between two successive values of U(N) is less than 10 or until fifty iterates have been computed. This iteration appears to converge quite rapidly, the iteration limit of fifty being much too higho Having obtained values for U(N) and T(N) from this iteration, P(N) is obtained from the momentum equation using P(N-1)? RHO(N) from the equation of continuity using RHO(N-1) and M(N) = U(N) * SQRT R *GAMT(MA) Finally, for those interested in recoding this progra.. it should be mentioned that convergence problems occurred whenever attempts were made to'improve' this iteration. Experience indicates that there is a fine line here between a convergent and a non-convergent iteration. For example, some of the attempted'improvements' did not seem to alter the logic, but were aimed only at reducing the amount of computation per iteration. CONDEN Arguments: PROGRAM COMMON VARIABLES, Group 1o PROGRAM COMMON VARIABLES, Group 2o 28

TSAT Computed saturation temperature. MDOT Rate of mass flow at the saturation point, ASAT SRHOAT USATU Length: 1124 octal Transfer Vector: FL, FRHOL, FSIGMA, IFLOW, NOZZLE, NUCLE1, ZERO Purpose: Supply to MAIN an initial point from which the calculations are to be started. The conditions at station zero are specified by X, A, P, T, RHO, U, M, NDOT, RADIUS, G and DELGo In order to permit the introduction of an externally computed initial point CONDEN first sets this point to the saturation conditions with no condensate or foreign particles0 It then gives a'READ DATAt statement which must be supplied with a value for EPSLONo If EPSLON = 0, the saturation conditions with no condensate or foreign particles and any overrides to the conditions at the initial point due to this'READ DATAT statement are returned to MAIN. If a non-zero value of EPSLON is obtained, CONDEN computes either the point of onset of condensation according to the c(T) criterion or, if this criterion cannot be satisfied, it finds the maximum point of DELG as computed by NUCLEI, This task is accomplished via a technique similar to the half-interval method and is most easily explained as follows0 1) Initialize DELG to zero, T to TSAT and DELT to DELT(1)o 2) T = T - DELT. 3) Use IFLOW and NOZZLE to obtain the flow conditions corresponding to an isentropic expansion to temperature To 4) Save DELG in DELG(1)o 29

5) Use NUCLE1 to compute the critical drop size, the nucleation rate for the volume element specified by A and DELX and the resulting mass fraction of condensate DELGo 6) If DELG is zero, return to step (2), otherwise compute A L(T) EPSLON DELA * (f(M)CP (T) ) * DELG where GAMMA- 1/M2 f(M) = MAX GAMMA - 7) If EPSLON exceeds EPS, the onset of condensation occurs between temperatures T + DELT and T and is obtained via a half-interval technique which begins at step (9)~ If EPSLON does not exceed EPS, continue the maximum DELG point search with step (8)~ 8) If DELG exceeds' DELG(1), return to step (2), otherwise the maximum DELG point has been passed. Set T = T + 2 * DELT and recompute the flow conditions (including DELG) for this temperatureo If DELT is less than or equal to EPS(1) print these conditions as the maximum DELG point and give an error return to MAIN, otherwise divide DELT by 10 and return to step'(2). 9) Initialize LOW to T and HIGH to T + DELTo 10) Set T = 5 * (HIGH + LOW) and duplicate the calculations of steps (3)9 (5) and (6) to obtain EPSLONo 11) If EPSLON exceeds EPS set LOW to T, otherwise set HIGH to To 12) If the absolute value of (1 - LOW/HIGH) exceeds EPS(2) return to step (10)o When this convergence criterion becomes satisfied, accept the temperature T and the associated flow conditions as those at the zeroth station, print them and give a successful return to MAINo 30

IFLOW Arguments: PROGRAM COMMON VARIABLES, Group 1 A Area corresponding to an isentropic expansion to temperature To M Mach number corresponding to an isentropic expansion to temperature To P Pressure corresponding to an isentropic expansion to temperature T. RHO Density corresponding to an isentropic expansion to temperature T. T Temperature for which the isentropic conditions are desired. U Velocity corresponding to an isentropic expansion to temperature T. Length: 167 octal Transfer Vector: SQRT Purpose: Compute the area, pressure, density, velocity and Mach number corresponding to an isentropic expansion to the given temperature To ISEN Arguments: PROGRAM COMMON VARIABLES, Group 1o A Area P Pressure T Temperature RHO Density M Mach number Length: 503 octal Purpose: ISEN is called from the internal subroutine IOCTRL of MAIN that controls the printing of the first output table. This subroutine computes PHAT, THAT, RHOHAT, PZERO3HAT and MI on the basis of the five 31

arguments and returns these values via the PROGRAM COMMON vector ZQ, ZQ(1) o.. ZQ(5) respectively. The first four of these quantities are easily computed from the value of MI, which is obtained via a modified half-interval technique so that the area ratio corresponding to MI lies in the interval (A/STAR - o 001, A/ASTAR +. 001) NOZZLE Arguments: PROGRAM COMMON VARIABLES, Group 1o Al Hollerith valued switch with the values'THROAT",'INVERSt or'AREA". This argument actually provides three entry points for NOZZLE. Throughout this report this subroutine is referred to both by name and by the value of the first argument enclosed in dollar signs0 A2 Unused when Al = $THROAT$o Nozzle coordinate for which the area is desired when Al = $AREA$. When Al = $INVERS$, A2 must be the Mach number on entry and will be returned as the nozzle coordinate with area A30 A3 Location for returning the throat area when - = $THROAT$, the area at A2 when Al = $AREA$, and when Al = $INVERS$ it is the area for which the nozzle coordinate is desired0 Length: 541 octal Transfer Vector: COS, SIN, SQRT Purpose: The origin of the nozzle coordinate is the throat and the coordinate increases positively in the direction of flow. In this system, the Xcoordinate of the intake is always negative and that of the exit always positive. Since the meaning of the throat area is basically different when using a wedge rather than conical nozzle, the type of nozzle is given as a Boolean input variable. 32

Al = $THROAT$o Compute the area at the X-coordinate in A2 and return it in A3. If A2 is not in the interval (XMIN, XMAX), a zero area is returned. Al = $INVERS$o Compute and return in A2 the nozzle coordinate corresponding to an area A3 and Mach number A2. If A3 is less than ASTAR, A2 is set to zero and A2(1) to 1. The nozzle coordinate will be computed in the diverging portion of the nozzle if the flow is supersonic and in the converging portion if it is subsonic. If the computed nozzle coordinate does not lie inside the nozzle, A2(1) is returned as zero; otherwise it is given the value 1. NUCLE1 [f tfl/: Arguments:' I' " ",'PROGRAM COMMON VARIABLES, Group 1. ". ~: // PROGRAM COMMON VARIABLES, Group 2. MDOT Rate of mass flow at the saturation point. Length: 252 octal Transfer Vector: ELOG, EXP, FSPILL, LNPSAT, RSPILL, SQRT, ZERO Purpose: Using the current values of D, SIGMA, RHOL, T(N), P(N), A(N) and DELX(N) compute: the critical drop size, a value for surface tension corrected for the intermolecular distance, the nucleation rate for the volume element specified by A(N) and DELX(N), and the resulting mass fraction of condensate due to the formation of these drops. The nucleation quantities are returned in RADIUS(N), NDOT(N) and DELG(N) and the corrected value of SIGMA replaces the old value. The values of SIGMA and RHOL should correspond to P(N) via FSIGMA and FRHOLo The subroutine FSPILL is used to control floating-point traps. This subroutine sets underflows to zero and transfers to the argument location if 33

an overflow occurs. All quantities are returned as zero if an overflow occurs, if RADIUS(N) exceeds 10 6 or if NDOT(N) is less than 1. The computations are based on the following equations. KBAR 2*MU RHOL * T(N) * R * (In P(N) - In P SAT RADIUS(N) = KBAR * SIGMA - D SIGMA 5 SIGMA - D/KBAR 4TEXP= SIGMA * RADIUS(N)2 3 K*T(N) NDOT(N) RHT) * 1 SIGMA * MU KK T(N),RHOL NA * exp (TEXP) * A(N) * DELX(N) DELG(N) = 3 * RHOL * NDOT(N) * RADIUS(N)3/MDOT where R is the universal gas constant, K is Boltzmann's constant, NA is Avogadro's number, and In PSAT is computed by LNPSAm Anild is the natural log of the saturation pressure corresponding to T(N). NUCLE2 Arguments: PROGRAM COMMON VARIABLES, Group 1o PROGRAM COMMON VARIABLES, Group 2. MDOT Rate of mass flow at the saturation point. Length: 157 octal Transfer Vector: FTSAT, NUCLE1, SQRT Purpose: Using the current values of ALPHA} D, L, RHOL, SIGMA, P(N), T(N) and U(N)3 compute the mass fraction of condensate formed in the 34

volume element specified by A(N) and DELX(N) due both to formation of tnew' drops and growth of'oldt drops. The values of RHOL and SIGMA should correspond to P(N), via FRHOL and FSIGMA, and L should correspond to T(N) via FL. This subroutine first calls NUCLEI to obtain RADIUS(N), NDOT(N) and the mass fraction of condensate due to the formation of these'new' drops. The radial increment of the drops initially formed at station I (DELRAD(I)) is not a function of the radius of these drops at station N - 1 (RADIUS (I)); the radial increments for all'old' drops is the same0 The equation for this uniform increment is 2 * K * NA 1/2 ALPHA * P(N) * DELX(N) * (T - T) * MU * T(N) L * RHOL * U(N) DROP where T DROP is the saturation temperature corresponding to P(N) as computed by FTSAT. The mass fraction of condensate due to this growth of the'old' drops is N-1 4 iDr__ * L NDOT(I) * DELRAD(I) * RADIUS(I)2 I=0 The value of DELG(N) returned is the sum of the two mass fractions of condensate. VAPOR FL FRHOL PSIGMA FTSAT LNPSAT Argument: ARG The entry VAPOR uses no argument. For the entries FL and LNPSAT, the argument is a temperature; for FRHOL, FSIGMA and FTSAT it is a pressure. Length: 637 octal Transfer Vector: ELOG, SQRT 35

Purpose: The entry VAPOR is used to supply, via a tREAD DATA' statement, the linear and quadratic approximations that form the computational basis for the other entries, iL eo, the SUBROUTINE VAPOR INPUT PARAMETERS. The entry FL computes the latent heat corresponding to the temperature ARG from the equation L + L(1) * ARGo The entry FRHOL computes the liquid density from the linear approximation RHOL + RHOL(1) * T where T is the saturation temperaDROP' DROP ture corresponding to the pressure ARGo The entry FSIGMA computes the surface tension from the linear approximation SIGMA + SIGMA * TDROP where TDROPis the saturation temperature corresponding to the pressure ARGo The entry FTSAT computes the saturation temperature corresponding to the pressure ARG from the linear approximation. 1/T = PSAT + PSAT(1) * In P, (L) if In ARG is less than or equal to PTP, and from the quadr"tic approximation 1/T = PSAT(2) + in P * PSAT(3) + in 2 * PSAT(4) (Q) otherwise. The entry LNPSAT computes the natural log of the saturation pressure corresponding to the temperature ARGo When the inverse of ARG exceeds TTP the result is obtained by solving (L) for in Po Otherwise, the roots of the quadratic (Q) must be computed. The proper root is chosen by taking the smallest root that exceeds PTPO 36

6. I/O STRUCTURE The following section describes how the input and output were implemented and organized in the program. The inflexibility of this aspect of the original program, Ref0 19 was one of the chief motivations for the rather drastic changes that have been made. The two most important changes concerning the I/O abilities of the program are reflected in greater flexibility in choosing the initial point and better organization in the output0 6. 1 INPUT STRUCTURE Since the M. A. D.'READ DATA' statement is used for obtaining the numerical input, an understanding of the data deck requires a knowledge of the function of this statement. The description of this statement given in the MAD manual is: "This statement causes information to be read from cards; no list of variable names or format specification is necessary. The values to be read and the variable names are punched in the data cards in a sequence of fields of the form: V1 1, V2 n2, * 0 o VK =nK* The V1, o o., VK are the variable names and nl o o0. nK the corresponding values. Reading is continued from card to card until the terminating mark * is encountered. Only the first 72 columns of a card may be used for data; o o0 However, as a convenience, the end of the card is treated as an implied comma 0. o For convenience in entering values of array elements, it is possible to designate only one variable name and have successive numbers, written without names, interpreted as the consecutive values of the array, io eo, V(J) n no 2 n2. o, nK 37

would be the same as V(J) = n1, V(J + 1) = n2 o o V(J + K - 1)=nK =" The use of this statement for input allows considerable flexibility, since no list or format is necessary, Thus for example, if a number of problems with the same nozzle geometry are to be run the nozzle parameters need only be read the first time. The remaining data sets for this subroutine could then consist only of a card with an * in column lo This situation is in fact true for all inputs, except EPSLON which is read by the subroutine CONDEN; that is, once they are read they maintain that value for the remainder of the run unless specifically changed by a subsequent tREAD DATA'O Finally, it should be remarked that many of the input parameters are assigned to storage in PROGRAM COMMON. Those variables assigned to PROGRAM COMMON have been divided into two groups0 Some of the subroutines have no access to these variables, some have access only to the first group and some have access to both groups0 In those subroutines that have access to some of PROGRAM COMMON and that also give'READ DATA' statements, care should be taken as to what is pun.nAed on the data cardso An example of an input deck and the resulting output have been included in Appendix B. The input deck was listed on an IBM 407. The problem required 9. 5 seconds of execution time on the IBM 7090. The input deck for a run consists of the program initialization deck followed by any number of problem decks. The program initialization deck structure is as follows: 38

(1A) Any number of 72-column comment cards0 For those familiar with FORTRAN, each card is read according to a format of 12A6 and immediately printed with a double space. (1B) A single card with'ENDbbb', where lb? denotes a blank, punched in columns 1-6. This card is used to terminate the comment card deck. (2) A deck of cards containing the SUBROUTINE VAPOR INPUT PARAMETERS. These parameters are obtained by a'READ DATA' statement. The variables that must be supplied at this point are PTP, TTP, PSAT o o o PSAT(4), SIGMA. o SIGMA(1), RHOL. o o RHOL(1) and L 0o o L(1). Only these variables may be given values in this deck. A program deck consists of essentially four parts0 (1A) Any number of 72-column comment cards, These are processed exactly as the comment card deck in the program initialization deck. They are double-spaced and are immediately preceded by the heading'BEGINNING OF PROBLEM NUMBER', which always begins a new page. (1B) A single card with'DATAbb' punched in columns 1-6o This card is used to terminate the comment card deck. (2) A deck of cards containing the PROGRAM COMMON parameters TZERO, PZERO, GAMMA, MU, CP, D and ALPHA and the entire set of MAIN PROGRAM INPUT PARAMETERSo These variables are obtained by a'READ DATA' statement. If any of these expected parameters do not appear, then they maintain their current value in the computer, The following is a list of the parameters that are loaded initially with the program and, hence, need never appear in the data deck unless these values are unsatisfactory: 39

XPOINT = 1 MNDELX =. 001 MNDELG = o0005 MXDELX = 10 MXDELG = 0015 MXDELX(1) = 1 Note that XPOINT is an integer and hence its value must be given without a decimal point. The remaining parameters are not preset by the program and must appear at least in problem 1o Since all problems of a run must deal with the same vapor, normally GAMMA, MU, CP, D and ALPHA will appear only in the first problem. (3) A deck of cards containing the PROGRAM COMMON parameter ASTAR and the SUBROUTINE NOZZLE INPUT PARAMETERS, again obtained by a'READ DATA' statement. Only these parameters should be included. Note that the two half-angles are specified in degrees, that XMIN is negative and that XMAX is positive. Since the parameter WEDGE is of Boolean mode, it has only two possible values; iL eo, OB.o true and 1B =. false0 (4) The final deck of cards is obtained via a'READ DATA? statement in the subroutine CONDEN and must supply a value for EPSLONo If EPSLON is given a non-zero value, the only other parameters that should occur are DELT(1), TRANGE and EPS 0 0 0 EPS(2)o These five parameters are preset, however, and thus need be supplied only if the values DELT(1) = 5o 0 EPS = o001 TRANGE = 100. 0 EPS(1) =. 01 EPS(2). 00001 are unsatisfactory0 If EPSLON is given a value of zero, however, then any of the parameters specifying the flow conditions at the initial point are legitimate inputs; i. e. X, DELX, A, P, T, RHO, U, M, RADIUS, NDOT, G, or DELGo If not supplied, these values will correspond to those at the saturation point0 40

6. 2 OUTPUT STRUCTURE The output structure is most easily understood by studying the example output given in Appendix Co The first page for each run consists of the comment cards in the program initialization deck and a list of the input supplied to VAPOR. Each problem is headed by the comment'BEGINNING OF PROBLEM NUMBER' at the top of the first page, the program supplying the problem number. This heading is followed by any comment cards for the problem. Following these comments is a complete list of the program parameters read by MAIN and NOZZLE ($THROAT$), the saturation conditions printed by MAIN, and finally a comment about and the conditions at the initial point, both printed by CONDEN. The remainder of the output consists of the three output tables: the first one printed while the calculations are being done and the last two after the problem has been terminated. Each station for which the index is an integral multiple of XPOINT produces one line of output in each of the three tables, Each page of a table consists of a heading for the columns followed by at most 27 doublespaced lines of output, The headings for the columns of TABLE 1 are: X Nozzle coordinate. (cm) PHAT Ratio of the computed static pressure to the corresponding isentropic pressure. THAT Ratio of the computed static temperature to the corresponding isentropic temperature, RHOHAT Ratio of the computed static density to the corresponding isentropic density, PZERO'HAT Ratio of the computed total head pressure behind a normal shock wave to the corresponding value for an isentropic flow, 41

M Computed Mach number. MI Isentropic Mach number for the corresponding area ratio. U Computed velocity, (cm/sec) G Computed mass fraction of condensate multiplied by 100o This scale factor means that G actually represents the percent condensate0 The headings for the columns of TABLE 2 are: X Nozzle coordinate0 (cm) P Computed static pressure. (dyne/cm ) PBAR Ratio of the computed static pressure to the supply pressure0 T Computed static temperature0 (OK) TBAR Ratio of the computed static temperature to the supply temperature0 RHO Computed static density. (gm/cm3) RHOBAR Ratio of the computed static density to the supply density0 A Computed area. (cm2) ABAR Ratio of the computed area to the throat area, io eo, the geometric area ratio0 The headings for the columns of TABLE 3 are: X Nozzle coordinate0 (cm) DELX Nozzle coordinate increment since the last station. (cm) G Same as the last column of TABLE 1. DELG Increase in mass fraction of condensate since the last station times 100o 42

RADIUS Critical drop size at this station. (cm) NDOT Number of drops per second formed between the last station and this one. RHODRP The number density of the drops. (drops/cm3) MEAN RADIUS The arithmetic average radius of all drops in the flowo (cm) 43

7. REFERENCES 1. Griffin, J. L., "Digital Computer Analysis of Condensation in Highly Expanded Flows", (Univ. of Mich. ) USAF Aerospace Research Labs, OAR, Report ARL 63-206, November 1963. 2. Sivier, K. R., "Digital Computer Studies of Condensation in Expanding One-Component Flows", (Univ. of Mich. ) in preparation, to be published as a USAF Aerospace Research Labs., OAR, report. 3. Arden, B., Galler, B., Graham, R., "The Michigan Algorithm Decoder, " Univ. of Mich. Press, November, 1963. 44

APPENDIX A PROGRAM LISTING 45

SCOMPILE MAD,PRINT OBJECT,PUNCH OBJECT DIRTY101 999069 06/24/64 5 59 45.5 PM MAO (12 MAR 1964 VERSION) PROGRAM LISTING......... AERONAUTICAL ENGINEERING NOZZLE CONDENSATION RESEARCH MAIN PROGRAM PROGRAM COMMON ASTAR,PZERO,TZERO,RHZERO, *001 1 GAMMA,MUtCPL,SIGMA,RHOL, *001 2 D,B,C,ALPHA, *'001 3 N,PPOINT,LINENC,PERPG,ZQ *CO1 DIMENSION ZQ(10) *C02 INTEGER N,PPOINT,LINENO,PERPG *0G3 PROGRAM COMMON X,CELX,A,P,T,RHO,U,M, *C04 1 G,DELG,RADIUS,NOOT,OELRAD *004 DIMENSION X(5CO),DELX(500),A(500), *005 1 P(5CO),T(500),RHO(500),U(5000),M(50), *005 2 G(500),DELG(5C00),RAOIUS(500),NDOT(500),DELRAD(500) *0C05 DIMENSION SVRAD(500),RHODRP(500),MDRAD(500) *'06 PROGRAM INPUT AND INITIALIZATION PRINT COMMENT $2NOZZLE CONDENSATION RESEARCH PRCGRAM$ *007 PCMNT REAC FORMAT VVA,CMNT(I)...CMNT(12) *C008 WHENEVER CMNT(1).E.$END $, TRANSFER TO PCMNTi *009 PRINT FCRMAT VVB,CMNT(1)...CMNT(12) *010 TRANSFER TO PCMNT *Ci1 PCMNT1 PROBNO=O *012 VAPCR.(Z) *013 PROBLEM INPUT AND INITIALIZATION START PROBNC=PROBNO+1 *014 REAC FORMAT VVA,CMNT(1)...CMNT(12) *015 PRINT FORMAT VVC,PROBNO *C16 PCMNT2 WHENEVER CMNT(1).E.$DATA $, TRANSFER TO START1 *C17 PRINT FORMAT VVB, CMNT(1)...CMNT(12) *018 REAC FORMAT VVA,CMNT(1)...CMNT(12) *019 TRANSFER TO PCMNT2 *020 START1 REAC CATA *C21 RHZERO=MU/TZERO*PZERO/8.314E+07 *022 THIS CATA SHOULD INCLUDE PZERO = SUPPLY PRESSURE (DYNE/CM**2) TZERO = SUPPLY TEMPERATURE (DEG. K) PE = ISENTROPIC EXIT PRESSURE. THIS IS USED AS AN INITIAL APPROXIMATION TO THE SATURATICN PRESSURE. GAMMA = RATIO OF THE SPECIFIC HEATS MU = MCOLECULAR WEIGHT OF THE VAPOR CP = SPECIFIC HEAT AT CONSTANT PRLSSURE 0 = INTERMOLECULAR DISTANCE IN THE LIQUID ALPHA = ACCOMODATION COEFFICIENT 46

DELX = IF XCO) IS THE FIRST POINT THEN X(1)=X(C)+DELX. THIS SAME DELX WILL BE USED TO OBTAIN X(K+1) FROM X(K) AS LONG AS THE OTHER PROGRAM PARAMETERS DC NOT INDUCE A CHANGE FOR DELX(K). XRANGE = LENGTH OF INTERVAL, STARTING AT THE CONDENSATION POINT, IN WHICH THE FLOW CONDITIONS WILL BE COMPUTED. MAXG = STOP COMPUTATION AT POINT WHERE G EXCEEDS MAXG. MING = WIEN G EXCEEDS MING THEN START VARIATION OF DELX. THIS DATA MAY ALSO INCLUDE XPOINT = TFE INITIAL POINT IS ASSIGNED A SUBSCRIPT OF ZERO, ONLY THOSE POINTS WITH SUBSCRIPTS XPOINT, 2*XPOINT,... WILL BE PRINTED. MNOELG = IF FOR STEP DELX FROM SOME POINT DELG IS LESS THAN MNDELG THEN DELX IS DOUBLED AND THE POINT IS RECALCULATED WITH THIS NEW DELX. MXDELG = IF FOR STEP DELX FROM SOME POINT DELG IS GREATER THAN MXDELG THEN DELX IS HALVED AND THE POINT IS RECALCULATED WITH THIS NEW DELX. MNDELX = MINIMUM VALUE THAT WILL RE GIVEN TO DELX REGARDLESS OF THE DELG PRODUCED MXOELX = MAXIMUM VALUE THAT WILL BE GIVEN TO DELX REGARDLESS OF THE DELG PRODUCED IF NDOT=0 -- OTHERWISE MXDELX(1) IS THE MAXIMUM ALLOWABLE VALUE. NOTE TIAT IF THESE PARAMETERS ARE ADJUSTED BY THE READ DATA STATEMENT THAT THEY W-ILL MAINTAIN THE ADJUSTED VALUE UNTIL CHANGED BY ANOTHER READ DATA STATEMENT EXECUTION. PRINT COMMENT $OCOMPLETE LIST OF PROGRAM PARAMETERS FOR THIS *023 1 PROBLEM.$ 0 23 PRINT FORMAT VIA,PZEROTZERORHZERO,PE *024 PRINT FORMAT VIBGAMMAMUCP *025 PRINT FORMAT VIC0D,ALPHA *026 PRINT FORMAT VIO,DELX,MNDELXMXDELXMXDELX(1) *027 PRINT FORMAT VIEMNDELGMXDELGMING *028 PRINT FORMAT VIFXRANGEMAXG *029 EXECUTE NOZZLE.($THROAT$,O.,ASTAR) *030 CALCULATE A SATURATION POINT BY MATCHING THE ISENTROPIC AND SATURATION CURVES. P=PE *031 Z=GAMMA/(GAMMA-1.) *032 S1 T=FTSAT.CP) *033 PSAT=PZEROC* (T/TZERO).P. Z ) *034 TEST=.ABS.i1.-PSAT/P)-.001 *035 WHENEVER TEST.L.O,,TRANSFER TO S2 *036 P=P+.5*(PSAT-P) *037 TRANSFER TO 51 *038 S2 TSAT=T *039 EXECUTE IFLOW. (ASATMSATPSATRHOSATTSATUSAT) *040 Z=MSAT *041 EXECUTE NOZZLE.($INVERS$,ZvASAT) *042 S3 XSAT=Z *043 MDOT=ASAT*RHOSAT*USAT *044 47

PRINT COMMENT $OSATURATION POINT SPECIFICATIONS — $ *045 PRINT FORMAT VIH,PSATTSATRHOSAT *046 PRINT FORMAT VII,MSAT,USAT *047 PRINT FORMAT VIJ,XSAT,ASAT *048 PRINT FORMAT VIKMDOT *049 WHENEVER Z(1).NE.1. *050 PRINT COMMENT $4 THE SATURATION POINT FOUNO IS NOT INSIDE *u51 1 THE NOZZLE, THE PROGRAM WILL CONTINUE.$ *051 TRANSFER TO START *052 END OF CONDITIONAL *053 CGMPUTE AN INITIAL POINT VIA SUBROUTINE CONDEN N=O *054 Z=CCNCEN. (TSAT,MDOT) *055 WHENEVER Z.E.C., TRANSFER TO START *C56 XL I IT=X+XRANGE *057 EXECUTE IOCTRL. *058 COMPUTE POINT N GIVEN POINT N-1 S29 ITER(1)=l *059 N=N+i *060 WHENEVER N.G.5C0 *061 PRINT COMMENT $OTHE NUMBER OF MESH POINTS HAS REACHED THE PRO *062 1 GRAMED UPPER LIMIT.$ *062 N=N-1 *063 TRANSFER TO END *064 ENC OF CONCITIONAL *065 DELX(N)=DELX(N-1) *066 RESTART POINT FOR COMPUTING POINT N IF. DELX(N) IS ALTERED FROM ITS INITIAL VALUE OF DELX(N-1) S30 X(N)=X(N-1)+DELX(N) *067 WHENEVER X(N ).G.XLIMIT *068 PRINT COMMENT SOTHE NEXT MESH POINT WOULD EXCEEC THE RANGE DE *069 1 SIREC.$ *069 N=N-1 *0C70 TRANSFER TO END *071 END OF CONCITIONAL *072 NOZZLE.[($AREA$,X(N),A(N)) *073 WHENEVER A(N).E.O. *074 PRINT COMMENT $ONOZZLE SUBROUTINE ERROR INDICATION - MESH POI *075 1 NT COES NOT LIE WITHIN THE NOZZLE.$ *075 N=N-1 *076 TRANSFER TO END *077 END OF CONDITICNAL *C78 CELA=A(N)-A(N-1) *079 DELG(N)=DELG(N-1) *080 G(N)=G(N-1)+DELG(N) *081 WHENEVER N.G.1 *082 DELL=(N-1 )-U(N-2) *083 OTHER I SE *084 DELU=.CCO1*U * C85 END OF CONCITIGNAL * 086 U (N) = U(N- 1)+DELU *087 ITER=I *088 BASIC ITERATION TO BALANCE NUCLEATICN AND DIABATIC FLOW EQNS. 48

531 CFLCW.(DELAEND) *089 RHOL=FRHOL.(P(N)) *090 SIGMA=FSIGMA.(P(N)) *091 TEST=O(N) *G92 NUCLE2.(MOOT) *093 WHENEVER ITER.E.1 *094 WHENEVER (NDOT(N).G.l.).AND. (DELX(N).G.MXOELX(l)) *095 DELX(N)=MXCELX(1) *C96 TRANSFER TO S30 *097 END OF CONDITIONAL *098 END OF CONDITIONAL *099 G(N)=G(N-1)+DELO(N) *100 TEST=.ABS.(G(N)-TEST) *101 WHENEVER TEST.L. C00001, TRANSFER TO S32 *102 ITER=ITER+1 *103 WHENEVER I'TER.LE.GITER, TRANSFER TO S31 *104 ITER(l)=ITER(1)+1 *105 OELX(N)=.5*OELX(N) *106 WHENEVER ITER(1).LE.5, TRANSFER TO S30 *107 PRINT FORMAT VVOX(N),PROBNO *108 TRANSFER TO END *109 VARIATION OF DELX(N) IF THE CALCULATED POINT N SHOWS TOO MUCH OR TOO LITTLE CONDENSATION SINCE X(N-1) S32 WHENEVER (O(N).L.MINO).OR. (ITER(l).G.1), TRANSFER TO S33 *110 1=0 *111 WHENEVER NDOT(N).G.1., I=l *112 WHENEVER DELO(N).L.MNOELG *113 TEST=2.*DELX(N) *114 WHENEVER TEST.G.MXDELX(I), TRANSFER TO S33 *115 OR WHENEVER DELG(N).G.MXDELG *116 TEST=.5*DELX(N) *117 WHENEVER TEST.L.MNOELX, TRANSFER TO S33 *118 OTHERWISE *119 TRANSFER TO S33 *120 END OF CONDITIONAL *121 DELX(N)=TEST *122 TRANSFER TO S30 *123 THE CURRENT VALUES AT POINT N HAVE BEEN ACCEPTED, ADJUST RADIUS(O)...RADIUS(N-1) FOR GROWTH SINCE X(N-1) S33 1=0 *124 S34 RADIUS(I)=RAOIUS(I)+DELRAD(I) *125 1=1+1 *126 WHENEVER I.L.N, TRANSFER TO S34 *127 WHENEVER G(N).G.MAXG *128 PRINT FORMAT VVEX(N) *129 TRANSFER TO END *130 END OF CONDITIONAL *131 EXECUTE IOCTRL. *132 TRANSFER TO S29 *133 THE COMPUTATIONS HAVE IN SOME WAY BEEN TERMINATED. PRINT THE LAST TWO TABLES AND CONTINUE ON TO THE NEXT PROBLEM. END PPOINT=N *134 PRINT FORMAT VVFPRDBNO *135 49

EXECUTE IOCTRL. *136 TABLE 2 -- XP,T,RHC,A AND THE ASSOCIATED BAR QUANTITIES I=O *137 PRINT FORMAT VVR *138 LINENO=1 *139 IOF1 PBAR=P(I)/PZERO *140 TBAR=T(I)/TZERO *141 RHOBAR=RHO(I)/PHZERO *142 ABAR=A(I)/ASTAR *143 PRINT FORMAT VVS,X(I),P(I),PBAR,T()), TBAR,RHO(I),R(HOARA(I),ABAR *144 LINENO=LINENO+1 *145 WHENEVER LINENO.G.PERPG *146 PRINT FORMAT VVR *147 LINENO=1 *148 END OF CCNDITIONAL *149 I=I+XPOINT *150 WHENEVER I.LE.N, TRANSFER TO IOF1 ~151 TABLE 3 -- XDELX ANO NUCLEATION INFCRMATION I=O *152 PRINT FORMAT VVT *153 LINENO=1 *154 ICF2 PRINT FORMAT VVU,tX(I)X(I) (I),GC(I),DELO(I),SVRAD(I),NDOT(I),RHODRP(T) *155 1 MERAD(I) *155 LINENO=LINENO+1 *156 WHENEVER LINENO.G.PERPG *157 PRINT FORMAT VVT *158 LINENC=1 *159 END OF CONDITIONAL *160 I=I+XPOINT *161 WHENEVER I.LE.N, TRANSFOR TO IOF2 *162 TRANSFER TO START *163 THIS INTERNAL FUNCTION MONITORS THE PRINTING FOR TABLE 1 ANC -- SAVES IN SVRAD(N) THE INITIAL RADIUS CF DROPS FORMED AT X(N) -- COMPUTES MORAD(N), THE MEAN DROP RADIUS AT X(N) -- COMPUTES RHOCRP(N)? THE DROP DENSITY AT X(N). INTERNAL FUNCTION *164 ENTRY TO ICCTRL. *165 WHENEVER N.G.O, TRANSFER TO IOC1 *166 PPOINT=O *167 TRANSFER TO IOC2 *168 ICC1 WHENEVER LINONO.G.PERPG *169 ICC2 LINENC=1 *17C PRINT FCRMAT VVP *171 END OF CONCITIONAL *172 WHENEVER N.E.PPOINT *173 ISEN.(A(N), P (N),T(N),RHO (N), M (N)) 174 PRINT FORMAT VVQ,X(N),ZQ(1)...Z(4),M(N),ZC(5),(J(Ni),G(N) *175 ZERC.(ZQ...ZQ(3') *176 THRCUGH IC3, FOR I=0,1,I.G.N *171 ZQ=ZG+NDOT(I) *178 ZQ(4)=NCOT(I) *17)' ZQ(1)=ZQ(1)+ZQ(4) 180 ICC3 ZQ(2)=Z() )+ZQ(4)*RADIJS( I) *181 RHOCRP(N)=ZQ/A(N)/U(N) *182 50

MDRAO(N)=ZQ(2)/ZQ(1) *183 SVRAC(N)=RADIUS(N) *184 PPOINT=N+XPOINT *185 LINENO=LINENO+1 *186 END OF CONDITICNAL *187 FUNCTION RETURN *188 ENC OF FUNCTION *189 DIMENSION CMNT(13) *190 DIMENSION (10),ITER(2) *191 INTEGER ITER *192 INTEGER CMNT,PROBNO,I *193 VECTOR VALUES XPOINT = 1 *194 VECTOR VALUES MNDELX =.001 *195 VECTOR VALUES MXDELX = 1.,.1 *196 VECTOR VALUES MNDELG =.0005 *197 VECTOR VALUES MXDELG =.C015 *198 VECTOR VALUES LITER = 10. *199 VECTOR VALUES GITER = 10. *200 VECTOR VALUES PERPG = 28 *201 VECTOR VALUES VVA =$12C6*$ *202 VECTOR VALUES VVB = $lH0,S6,12C6*$ *203 VECTOR VALUES VVC = $1H1,S40,H*8EGINNING OF PROBLEM NUMBER*,I3*$ *204 VECTOR VALUES VVD=$1HO,H*ITERATION FOR MESH PuINT*,F9.5,H* HAS FAILED T *205 1 0 CONVERGE. THIS TERMINATES PROBLEM NUMBER*,13*$ *205 VECTOR VALUES VVE=$1HO,H*PERCENT CONDENSATE EXCEEDS THE UPPER LIMIT'MA *206 1 XG' AT MESH POINT*,F9.5*$ *206 VECTOR VALUES VVF=$1HO,H*FINAL OUTPUT IN TABLE 1 FCR PROBLEM*,I3*$ *207 VECTOR VALUES VVP=$1Hl,S4,1HX,S 12,4HPHAT,S11,4HTHAT,S11,6HRHOHAT,SE, *208 1 9HPZERO'HAT,S4,tHM,S9,2HMI,S14,1HU,S11, 1HG*$ *208 VECTOR VALUES VVQ=$1HO,F10.5,4F]5.8,2FlOS,F15.3,2PF12.8*$ *209 VECTOR VALUES VVR=$1HltS4,HX,S1 2,1HPtS12,4HPBAR,S13,1HT,S12,4HTBAR,S11 *210 1,3HRHO,Sll,6HRHOBAR,Sll,lHA,S13,4HABAR*$ *210 VECTOR VALUES VVS=$1HO,FlO.5,8E15.7*$ *211 VECTOR VALUES VVT=$1HI,S4,1HX,SlO,4HDELX,S12,1HG,S13,4HDELG,S O,6HRAOIU *212 I S,S10,4HNDOT,S10,6HRHODRP,S7, 11HMEAN RADIUS*$ *212 VECTOR VALUES VVU = $IHO,F10.5,F12.8,S3,2(2PF12.8,S3),5E15.7*$ *213 VECTOR VALUES VIA=$HQOPZERO(DYNE/CM**2)Q,'PE10.3,S5,H*TZERO(CEG. K)*, *214 1 F0O.4,S5,HQRHOZERO(GM/CM**3)Q,LPElC.3,SS,HQPE(DYNE/rM**2)Q, *214 2 IPE10.3*$ *214 VECTOR VALUES VIB=$6HOGAMMA,F6.3,S5,H*MU(GM/GMOL),F7.2,S5,H*CP( YNE-CM *215 1 /GM-DEG. K)*,lPE10.3*$ *215 VECTOR VALUES VIC=$H*OINTERMOLECULAR DISTANCE*,lPE1O.3,S5,H*ACCOMOATIO *216 1 N COEFFICIENT*,F7.3*$ *216 VECTOR VALUES VID=$H*ODELX*,F6.4,S 5,H*MIN(DELX)*,F6.4,S5,H*MAX(DELX) IF *21T 1 NDOT=0*,F7.4, S5,H*MAX (DELX)*,F6.4*$ *217 VECTOR VALUES VIE=$H*OWILL BEGIN CHANGING'DELX' TO KEEP DELG GREATER T *218 1 HAN*,2PF7.4,H* PERCENT AND LESS THAN*,2PF7.4,H* PERCENT WHEN *218 2 G EXCEEDS*,2PF7.4,H* PERCENT**$ *218 VECTOR VALUES VIF=$H*OXRANGE IS*,F7.2,H* BUT WILL TERMINATE IF G EXCEED *219 1 S*,2PF7.2,H* PERCENT.**$ *219 VECTOR VALUES VIH=$H*OPRESSURE*,1PE10.3,S5,H*TEP.*, 1PE10.3,S5,H*CENSIT *220 1 Y*,1PE10.3*$ *220 VECTOR VALUES VII=$H*OMACH NUMBER*,F7.3,S5,H*VELOC ITY(CM/SEC)*,F1O.2*$ *221 VECTOR VALUES VIJ=$H*ONOZZLE POINT*,F10.4,S5,H*AREA*,F9.5*$ *222 VECTOR VALUES VIK=$H*OTOTAL MASS FLOW(GM-CM/SEC)*,F1O.5*$ *223 END OF PROGRAM *224 THE FOLLCWING NAMES HAVE OCCURRED ONLY ONCE IN THIS PROGRAM. 51.

$CCMPILE MAC,PRPINT CJECT,PUNCH OBJECT CFLOw1O1 999066 06/24/64 11 11 0.7 AM MAC (12 MAR 1964 VERSION) PROGRAM LISTIN(G......... AERCNAUTICAL ENGINEERING NCZZLE CONDENSATION RESEARCH CONDENSATION FLOW EXTERNAL FLNCTION (DELA,LOC) *001 STATEMENT LABEL LCC *002 PRCGRAM COOVCN ASTAR,PZERO,TZERC,RHZERC, *003 1 GAMMA,MU,CP,L,SIGMA,RHCL, *003 2,HB,C,ALPHA, *003 3 N,PPCINT,LINENO,PERPG,Z *003 LIMENSION ZQ(1C) *CO4 INTEGER N,FPCINT,LINENC,PERPG *005 PROGRAM COMMCN X,DELX,A,PT,RHO,U,M, *006 1 G,CELG,RADIUS,NDT, DELRAC *006 ~IMENSION X(SCC),DELX(500),A(5CC) *007 1 P(5C0C) T(500),RHC(500),U(5C00) (500), *007 2 G(CC),DELG(500),RADIUS(5CC),NOT(50C),DELRAD(500) *007 USING THE CURRENT VALUE OF N AND THE SUPPLIED VALUES FOR A(N) CELA(N),G(N) AND DELG(N) CALCULATE T(N),P(N),RHC(N) AND U(N) SATISFYING THE DIABATIC FLOW EQUATIONS. THE CONCITICKS AT PCINT N-1 AND AN INITIAL APPROXINATICN TC L(N) aRE REQUIRED. ENTRY TO CFLCI. *C08 NM1=N-1 *009 CELU=t (N)-U(NM1) *010 UNMG=1.-G(N) *011 MCCG=EELG(N)/UNMG *012 CQA=CELA/A(N) *013 L=FL.(T(NM1)) *014 THRCUCH SI, FOR I=1,1,I.G.50 *015 CELT=(L*CELG(N-U)-(N)*DELU)/CP *016 T(N )=T(NM1 )+CELT *017 L=FL. (T(N)) *018 DELU=CQA+MCQG+CELU/U(N)-DELT/T(N) *019 CELL=CELU*R/MU*T(N)/U(N)*UNMG *020 T=L(NM 1)+DELU *021 TEST=.ABS.(1.-UT/U(N)) *022 L (N)=LT *023 S1 ~HENEVER TEST.L.EPSLON,TRANSFER TO S2 *024 ITERATION LIMIT HAS BEEN REACHEC WITH NC CONVERGENCE PRINT COMMENT $'CFLOI' - ITERATION FAILED TC CCNVERGE$ *025 TRANSFER TO LCC *026 ITERATION HAS CONVERGED S2 P(N)=P(NM1)/(1.+(U(N)*CELU*MU )/(T(N)*R*(1.-G(N) ))) *027 RHC(N)=RHC(NM1)/(1.+DELU/U(N)+DELA/A(N)) *028 M(N)=L(N)/SQRT.(R*GAMMA/MU*T(N)) *029 FUNCTION RETURN *0C30 VECTCR VALUES R=8.314E+07 *031 VECTOR VALUES EPSLCN=1.E-C6 *032 52

I rN TE. GE P, NM I *C33 END OF FUNCT ION *034 53

$COMPILE MAD,PRINT OBJECT,?UNCH OBJECT C(NONOC1 99910'6 06/26/64 2 29 48.3 PM MAD (12 MAR 1964 VERSICN) PROGRAM LISTING........ AERONAUTICAL ENGINEERING NOZZLE CONDENSATION RESEARCH INITIAL/CONDENSATION POINT EXTERNAL FUNCTION (TSAT,,MDOT) *'l01 PROGRAM COMMON ASTAR,PZERO,TZERO,RHZERO, *u02 1 GAMMA,MU,CP,L,SIGMA,RHOL, *002 2 D,B,C, ALPHA, *0 2 3 N, PPOI NT,L INENO, PERPPG ZC' *1 02 DIMENSION Zw( 10) *;)3 INTEGER NPPOINT,L INENO,PERPG * n04 -PROGRAM COMMON X,DELX, A,P,T,RHOUMv *0, M 1 G DELG, RA I US,NDOT, DELR.AD0 *0 5 DIMENSION Xi( 50),DELX 500),A(500), *006 i P(50), T ( 500),RHO (500),U(500 ), M (50) C 2 G( 5C ) tELG(500),RAD I US( 5!00),NO 500 ),ELRAD ( 50) * ENTRY TO COiDEN. *;'7 SET INITIAL POINT TO SATJURATION POINJT vITH NO DIRT OR DROPS. N=C *- 8 T=TSAT *0g IFLOW. (A,M,P,RHO,T,U) * 1 C X=M.c 11 NOZZLE. ( INVFRS$,X,A) *[!12 ZFRO. (NOOT~RADIUS G,OELo) *(:13 ThIS REA4 MAY ALTER ANY OF THE VALUES AT THE INITIAL PO!INT OR MAY CAUSE THE PRKiGRAM TO COMPUTE A'ClNOE\lISATI GrN PtI iNT EASED ON THE EPSILON(PT) CRITFRIA. REAL DATA C 14 WHENEVER< E PSLONi NF.C., TRANqSFER T O SZECt1 *;- 15 PRINT CO]MMENT $.0 THE CALCULATIONS ARE T( BE STAR<TED AT TH *'16 I E SATURATION O)(INT WITH ALL NUCLEATION QlUANTITIES ZEiO.s$ *C1O PRINT CO)MME'dT $ A'READ DATA' STATEiMENT FOLLuWS ALL TH-ESE CO" *-17 1 PUTATIONS CONSEGUENTLY AN3Y Ot ALL OF THE FC)LLflWING- MiAY Fte -3 *C 17 PRINT COVMMENT $ CHAN;,IGEO - X,A,P,RiIO,TtU,,M,R/AN)IUSt,! O[) tG,DELGL *'- 10 PRINT FOR MAT VVA,X,AP,RHO,T,U, M *; 1 PRINT FO'Ri. MAT VVo?RADIUS,lNDOTG *'02lFUNCTION RETURN X(1) *;21 SLERO1 PRINT COMrtMENT $C30 DETERMINE A CONDENSATITON O(gN TIlE bASIS O!F T *32' 1 HE VALUE OF EPSILGN(P,T), WHICH EVER APPLIES.$ *022 TM IN=T-TRANGE * -.2 3 WHENEVER TMIN.L.O., TMIN=OD. *C4 CELT=CELT(i) *2 25 LELG=..= 2 6 MAXIMUM OELG PlO IN.T SEARCH -- IF EPSILON((PT) EXCECOI)SPS, JUMP IMMEDIATELY TO HALF-INTERVAL FOR'COINDENSATION' POINI S1 T=T-CEFLT * 27 wHENE VER'T.L.TIN, I I TRANSFERR O S3 *-28 IFLOW. (A,M,P,RHOrT,U) *T27 54

NOZZLE. ($ INVERS, X, A) *3 vWHENEVLE3 X(1).E.C., TRANSFER TU S3 *332 KHGL=FRHOL. (P) * 33 S iG(MA=FS IGMA. (P) *334 CELG ( 1 ) =ELG NUCLEi. (MDOT) * 36 wHENEVER DELG.E.O., TRANSFER TO SI * -37 L=FL. ( T) 38 NOZZLE. ($ARZAI,X-DELX,A(I) ) LELA=A-l 1) * C EPSLCN=(GAMMA-I./M/M)/(GA MMA-1.) *41 WHENEVER EPSLON.L.1.,EPSLONl= 1o *42 EPSLON=A* (I.-L/CP*EPSLON/T) /DELA*I)ELG 4 WHENEVER.AoS.EPSLON.G.EPSv TRANSFER T(] S5 S2 wHEENEV'R DELG.(;..DELG(i), TRANSFER TO S1 4 C S3 T=T+CELT+OEL r *46 TEMPERATUE- T IS MAXIMUM POINT nF DELG. TESTS ON EPSLON(PT) INDICATE THAT CON!OErJSATION IS NOT SUFFICIENT TO CALCULATE IN THE MAIN PROGRAM. S4 IFLOW. (A,M.P,RHO,T,tJ) *47 X=tI * 48 NOZZLE. ($1IVERSl,X,A) *349 RHOL=FRHOL. (P) *35( SIGMA=FSIGMA. (P) *I51 NUCLEi. (MDOT) *i52 WHENEVER DELT.G.EPS(1) *53 LELT=CELT*.1 *54 TRANSFER TO S1 *355 END OF CONDITIONAL * 56 PRINT COMMENT $? INSUFFICIENT INITIAL CONDEcNSATICN TO INITI *357 ATE CALCULATION.$ *'57 PRINT FORMAT VVA,X,A,P,RHO,TU,M *C58 PRINT FORMAT VVR,RADIUS,NDOT,DELG *'59 FUNCTION RETURN O. *C6 HIGH IS A TEMPERATURE AT WHICH DELG.NE.O. AND WHERE EPSLON DOES NOT EXCEED EPS. TEMPERATURE T IS SUCH THAT EPSLGN EXCEED EPS. BEGIN HALF-INTERVAL SEARCH FOR'CONDENSATICN' PCINT S5 LOW=T *361 H IGH=LOW+DELT f*62 S6 T=.5*(HIGH+LOW) *)63 IFLOW. (A,M,PRHO,T,U) *364 X=M1 *065 NOZZLE. ($INVERS$,X,A) *066 L=FL. IT) *67 RHOL=FRHOL. (P) *368 SIGMA=FSIGMA. (P) *369 NUCLEI. (MDOT) *C73 NOZZLE. ($AREA$,X-DELX,A(1)) *0,71 CELA=A-A ( 1 ) *072 EPSLON= ( GAMMA-1./M/M ) / GAMMA-1.) *)73 WHENEVER EPSLON.L.1., EPSLON=I. *374 EPSLCN=A* 1.-L/CP*EPSLON/T)/DELA*DELG *075 WHENEVER.ABS.EPSLON.G.EPS *0,76 LOW=T *077 OTHERWISE *078 HIGH=T *079 ENO OF CONDITIONAL'080 55

WHENEVER.ABS.(1.-LOW/HIGH).G.EPS[2), TRANSFER TO S& *C81 CONDENSATION POINT HAS BEEN FOUND G= DELG *032 PRINT COMMENT $0 SUFFICIENT INITIAL CONDENSATION TO INITIAT *083 1 E THE CALCULATIONS. CONDENSATION POINT FOLLOWS.$ *083 PRINT FORMAT VVA,X,A,P,RHO,T,U,M *584 PRINT FORMAT VVB,RADIUS,NDOT,DOLG *085 FUNCTION RETURN 1. *0 86 VECTOR VALUES DELT = 0.,5. *087 VECTOR VALUES EPS=.001,.01,.00001 *C 88 VECTOR VALUES TRANGE = 100. *C89 VECTOR VALUES VVA = $1HO,2HX=,F8.4,3H A=,F9.4,3H P=,E1l.5,5H KtiO=,E11.5 *9c, 1,3H T=,F9.4,3H U=,F9.2,3H M=,F8.4*$ *0s90 VECTOR VALUES VVB = $1HO,7HRADIUS=,lPE12.5,6H NDUCT=,lPE12.5,3H G=, *91 1 1PE12.5*$ *091 END OF FUNCTION *92 THE FOLLOWING NAMES HAVE OCCURRED ONLY ONCE IN THIS PROGRAM. COMPILATION WILL CONTINUE. S2 S4

$CCVPILE MACPRINT CPJECTPUNCH ORJECT IFLOWIOL 999&b6 06/24/64 Ii 10 42.8 AM MAC (12 MAR 1q64 VERSICN) PRCORAM LISTING AERONAUTICAL ENGINEERING NOZZLE CONDENSATION RESEARCH ISENTROPIC FLOW EXTERNAL FUNCTION (AlA2,A3,A4,A5,A6) *001 PROGRAM COMMON ASTARPZERO,TZERORHZERC, *C02 I GAMMAMUCPLSlGMA,RHL, *002 2,IBC,ALPHA, *002 3 NPPCINTLINENCPERPG,Z. *002 CIMENSION ZC(IC) *003 INTEGER NFPPCNTINLNENCPERPG *004 ENTRY TO IFLC%. *005 T=A5 *006 PRESSURE EXP=GAVMA/(GAVMA-1.) *007 Z=(T/TZERC).P. EXP *008 A3=PLRCR*Z *009 VELOCITY AND MACH NUMBER Z=2.*(T ERO-I) *010 A6=SCRT.(CP*Z) *011 Z=Z/(GAMMA-1.) *012 A2=SCRT.(Z/T) *013 CENSITY AND NOZZLE AREA EXP=1./(GAPVA-1.) *014 Z=1.+.5*A2*A2/EXP *015 A4=RHZERO/(Z.P. EXP) *016 EXP=.5*(GAMMA+1.*E:XP *017 Z=2.*Z/(GAIMA+1.) *018 Z=(Z.P. EXP)/A2 *019 A1=ASTAR*Z *0,20 FUNCTION RETURN *021 END CF FUNCTION *022 57

$CCMPILE MAC,PRINT OBJECT,PUNCH OBJECT ISEN 101 999066 06/24/64 11 11 22.0 AM MAC (12 MAR 1964 VERSICN) PROGRAM LISTING. AERCNAUTICAL ENGINEERING NOZZLE CONDENSATION RESEARCH EXTERNAL FUNCTION (A,PtTEMP,RHO,M) *001 GIVEN THE ABOVE QUANTITIES, THIS SUBROUTINE CCVPUTES AND PLACES IN PROGRAM COMMON MI = ISENTROPIC MACH NUMBER CCRRESPCNOING TO GIVEN AREA RATIC - ABAR. PHAT = RATIC OF COMPUTED STATIC PRESSURE TC CCRRESPONDING ISENTROPIC PRESSURE. TI-AT = STATIC TEMPERATURE RATIO. RFCFHT = STATIC DENSITY RATIO. PZPIAT = PITOT PRESSURE RATIO. PROGRAM CCMMON ASTAR, PZERO,TZERO,RHZERC, *002 1 GAMMA,MU,CP, L,SIGMA,RHCL, *002 2 C,B,C,ALPHA, *002 3 N,PPO I NTI, LINENC, PERPG, ZQ *002 CINENSICN ZCQ(1) *C03 INTEGER N,PPC IN1,LINEN,PERPG *004 INTERNAL FUNCTION FOR MI IN TERMS OF ABAR. INTERNAL FUNCTION F.(ARG)=(T(1)/ARG)*((1.+T(2)*ARG*ARG).P. *005 T) - ABAR *005 ENTRY TC ISEN. *006 PBAR=A/AS4AR *007 CCIPUTE ISENTROPIC MACH NUMBER FOR GIVEN AREA RATIO T= (CAMMA+1.) / (2.*(GAMMA-1.) ) *008 T(1)=(2./(GANMA+1.)).P. T *009 T ( 2 ) =.5*(GAMMA- 1. ) *010 FI=((GAMMA+1.)/(GAMMA-1.)).P. T *011 I= ( 1.+(ABAR-1. )*MI).P. T(2) *012 II-ENEVER MI.L.1.4 *013 CCNTINUE *014 CR WI-ENEVER MI.L.4. *015 MI=.E*MI *016 CA WI-ENEVER I.L.6. *017 PI=.9*MI *018 CR hFENOEVER MI.L.1C. *019 MI=. 5*MI *020 ENC CF CCNDITIONAL *021 T(3) F. (I) *022 WHkENEVER.ABS.T(3).L. EPS, TRANSFER TC S3 *C23 CELM=. ( T ( 3 )/.ABS. T ( 3 ) *024 S1 T(4)=F. (MI-ELM ) *025 WHENEVER.ABS.T(4).L. EPS, TRANSFER TC S2 *026 WHIENEVER T(3)*T(4).G. 0. *027 M I=: I-CELM *028 CR.I-FENEVER T(3)*T(4).L. 0. *029 CELM=. 5*EL *030 58

ENC CF CCNCITIONAL *C31 TRANSFE-R IC Si *C-32 S2 I= I-CELV *033 FALF-INT —RVAL TECHNIQUE FOR THE ISENTROPIC MACH NUMPER HAS CCNVERGEC. CCMPLUT ISENTROPIC DENSITY, PRESSURE ANC TEMP. S3 T1I=1./(I.+T(2)*MI*M) *034 PI= TI.P. (GANlA/(GGAMMA-1.)) *035 R- I = TI.P. (1./(GAMMA-1. ) ) *36 CCMPUTE PRESSLRE RATIOS NECESSARY FCR CCMPLTING THE PITOT FRESSURE I2ATIC. I=CAMMA+1. *C37 T 1) =CAMMA-1 *038 1 2 )=2.*GCAMP *039 R=[(T/(T(2)*f*M-T(1))).P, (1./T(1)) *040 KM I =( T ( 2 ) I *M I- T ( 1))).P. ( 1./ ( 1 *041 T (2) =G AA A/ ( 1 ) *042 RNI=tRMI*((T/(T (1) +2./(MI*II ))).P. T(2) *043 KR:M=R*(.5*T**M).P. T(2) *044 CCMPUTE FINAL CUANTITIES DESIRED AND PLACE IN PRCGRAM COMMON LZ);( 1)=P/PZt'RE/PI *045 Z ( 2 )= TM P / T ER / T I *046 ZC( 3)=RHC/Rti ZCRc/HOI *047 L ( ) =P//2 ZERC*R /RM 1 *048 Z ( 5 ) I1 *049 FUNCTION RETLRN *050 VECTCR VALCUS EPS = 1.E-03 *051 LCIENSIF1\N T(7) *052 ENC CF FtUNCTION *053 59

$CCE" PI LE fv AC, PRlI N I (ir9J ECT, PU.'\Ch C' i JECCT F) Z LElI C 1 999066 0\6/24/64 11 1 30.9 AM MAA (12 MAR 1 64 RE SiC ) P iE AC L [STING. hERCN UT[CAL ENG I NE I RU NCLLLC CONDll-)ENSATION RESEARCH N CL LL E CRICIN OF THE COPDIREATE SYSTE~t FCR THE NCLZLE IS THE CENTER CF TH-E THRCAT, I.E., AlE) = ASTAr<. THE VARIABILE X INCRiEASES PCESITIVCLY IN THF CIRCCTION OF FLCWi AND IS NEGATIVE UPSTREAM F kCC THC IF<RCAT T hE NOZZLE PARAMETE7RS AR ASTAR = THRCAT AREA x YMI = X-CCORDINATC CF INTAKE XMAX = X-CECkOINATE CF EXIT I'NAX.:i = INTAKE HALF ANGLE IN CEGREeS CLTANG = FXIT SID HALF ANGLE IN DEGREES WECGL = COLEAN VARIABLE, OP CR 1, INDICATING wHETHER T[L NCLLLF IS CONICAL CR wOCGE RESPECTIVELY. THE ENTRY THREAT READS AND PRINTS THESE PARAMETERS. EXTERNAL FUNCTION (A1,A2,A3 * 001 PRCCRAM COPvMCN ASTAR,PZERO,TZERC,RHZERO, *0 02 1 GAMMA,MU,CPL, SIGMA RHCL, * CC2 2 DECALPtHA *002 3 NPPC[NT, LINENCPERPGZQ *002 CI MNSION Z I2(1C) *003 INTEGER NPPCINT,LINENCPERPG *G04 ENTRY TO NCLZLE. *CO5 WHENEVER Al.E. STHROAT$ *CC6 RETURN THREAT AREA AND INITIALIZE IF NECCESSARY REAC CATA *007 tO-ENEVER WECGE *008 PRINT COMMENT $OA WEDGE NOZZLE WITH THE FOLLOWING SPECS. IS H *0CC 1 EING USEC — $ *C C 9 CSTAR=ASTAR *010 0TH ER WISE *011 PRINT COMMENT $0A CONICAL NOZZLE ciITH THE FOLLOWING SPECS. 1S *12 1 BEING USED — $ *012 CSTAR= l.1283792*SQRT.(ASTAR) * 013 ENC CF CONCITIENAL *014 PRINT FORMAT VVAASTAR, INANGXMIN *015 PRINT FORMAT VVBOUTANGXMAX *016 A= RANG/RAEIAN *0 017 INTAN=SIN'.(A)/CS.IA) *01 8 A=0UT A NO/RADIAN *019 CUTTAN=SN'.(A)/CCS.(A) *020 A3=ASTAR *021 FUNCTION RETURN *022 OR WHENEVER Al.E. $AREA$ *023 CCMPUTE THE AREA AT A2 AND RETURN IN A3 WHENEVER (XMIN.L.A2).AND. (A2.L..0) *024 60

A=INTAN*(. AS. A2) *025 GR WFENEVER (A2.G.0.).AND. (A2.L.XMAX) *026 A=OUT TAN*A2 *C27 CTHER I SE *028 3=C. *029 hFENEVER A2.E.O.,A3=ASTAR *030 FUNCTION RETLRN *031 ENC OF CONDITICNAL *032 WhENEVER ECCGE *033 A3=ESTAR+2. * *034 CTIERhISE *035 A=DSTAR+2.* *036 A3=.8a*/ *. 78 5382 *037 ENC CF CONDITIONAL *038 FUNCTION RETLRN *039 OR WHENEVER Al.E. $INVERS$ *040 CCMPUTE THAT POINT ON THE X-AXIS OF THE NCZZLE THAT FAS AREA'A3'. THIS POINT DEPENDS CN THE MACH NUHPER WAHICH MUST BE IN A2 ON ENTRY, SEE THE PROGRAM FCR TFE PARTICULARS. IF A2(1) IS 1. CN RETURN THEN THE SCLUTION IS IN A2, IF IT IS ZERC THEN THE SCLUTION IS STOREC IN A2 BUT DOES NOT LIE INSIDE THE NCZZLE. WFENEVER hEGGE *041 A=.5 *(A3-ASTAR) *042 CTHER I SF *043 A=.564189s*SCRT. (A3)-.5*DSTAR *044 END OF CCONITIONAL *045 WIFENEVER A.L.C.,A=C. *046 X( 1 )=l. *047 A C = A 2 *C048 WFENEVER MACt-.L.1. *049 X=-A/ INTAN *+50 WIENEVER XMIN.G.X,X(1)=0. *051 CTHER1ISE *052 X=A/CUTTAN *053 tFENEVER X.G.XMAX,X(1)=O. *C54 ENC OF CCNCITIONAL *055 A2=X *056 42(1 )=X( ) *057 FUNCTICN RETLRN *058 ENC OF CONDITIONAL *059 INTEGER Al *060 PCCLEAN WEDGE *061 CIVENSION X(1) *062 VECITC VALUES RADIAN=57.2957795 *063 VECTCR VALUES VVA=$H*OTHROAT AREA*,F9.5,S55,H*INPUT riALF-ANGLE(DEG.)*, *064 1 FS.5,S5,H*INPUT LFNtTH*,F9.4*$ *064 VECOTCR VALUES VV E=$ I C,S24,H * 1UTPUT HA L F-ANGLE(DEG.)*,F9.5,S4,H*.1UTPUT *0C65 1 LENC-TF*, F.4*$ *065 ENC CF FLNCTIrN *C66 61

$CCEPILE vAC,PRkINT OBJECT,PUNCH OBJECT NUCLllOl 999066 06124/64 11 11 8.1 AM MAC (12 MAR 1%64 VERSICN) PRCGRAM LISTING......... AERCNAUJTICAL ENGINEERING NOZZLE CONDENSATICN RESEARCH NUCLEATICN THEORY - 1 EXTERNAL FLNCTICN (MOOT) *001 PROGRAM COMMON ASTAR,PZERO,TZERO,RHZERC, *002 ~~1 ~ GAMMA, MU CP,L, SIGMA, RHCL, *002 2 CB,C,ALPHA, *002 3 NP PO I NT, L I NENC,PERPG tZQ *002 CIMENSIGN ZC(10) *003 INTEGER N,FPC INT, L INENC, PERPG *004 PROGRAM CCIMVON X,CELX,A,P.T,RHO,U,M, *005 ~1 ~~G,DELGRADIUS,NDOTDELRAC *C05 CIMENSION X( 5CO),DELX( 500),A(500), *006 1 ~ P(5CC),T(500),RHO(500C),U(500),M(50C), *006 2 G(5CC),OELG(500),RADIUS(5CO),NCCT(500),DELRAO(500) *006 ENTRY TO NUCLEi. *O07 AT=T(N) *008 CCMPUTE RADILS ( CV. ) OF THE DROPS THAT ARE CONDENSING FSPILL. (51) *009 KPAR=R*(ELOG. (P(N)) - LNPSAT. (AT)) *010 KBAR=2./RHCL*MU/AT/KBAR *011 h,-ENEVER KBAR.E.C., TRANSFER TO S1 *012 ARAC=KBAR*S IGMA-O *013 S IGMA=S IGMA-C/KBAR *014 WI-ENEVER ARAC.G.1.E-06, TRANSFER TO S1 *015 WFENEVER (SIGMA.L.C.).OR. (ARAD.L.C.) *016 S1 ZERO. ( S IGMA,RADIUS(N),NDOT(N),DELG(N)) *017 RSPILL. (C) *018 FUNCTION RETURN *019 ENC CF CONDITIONAL *020 COMPLTE THE NUMBER OF DROPS OF THIS SIZE THAT ARE CCNCENSING PER CENTIMETER CUBED PER SECOND. TEXP=-(4.1E879*ARAC/K*ARAO/AT*SI[GMA) *021 ANDOT=SQRT.(SIGMA*MU/NA(1))/K(I)*P(N)/AT *022 ANCOT=EXP.(TEXP)/AT*ANDOT/RHOL*P(N)/K *023 NDOT(N )=.7S78846*ANDOT*A(N)*DELX(N) *024 WHENEVER NDOT(N).L.i., TRANSFER TO S1 *025 COMPUTATION OF PERCENT OF LIQUIC MASS, DELTAG LELG(N)=4.18E7S*RHOL*NDOT(N)*ARAD*ARAD/MCCT*ARAC *026 RAC ILS (N)=ARAD *027 RSPILL.(C) *028 FUNCTION RETLRN *029 62

VEC ICR VALUES NA=6.027E+23,6.G27E+C3 *030 VECTCR V'ALUES K=1.379E-16,1.379E-C6 *031 VECTCR VALUES R=8.314E+C7 *032 ENC CF FUNCTION *C33 63

$CCt ILE- MaC, tFP.I\T I J\: CTP CHi ( i EJLCT;U'CL21O1 999.; 6, 06/24/64 11 11. 15.5 AM MAC (12 VAt 1.64 vE V I' ) IN CC, tA' LITI(....... ARCFr.A UT I CsL ENGINER- ING NCZZLL' CCNl,)EN[SATI, N RESECAiCh NUCLEATICN THECRY - 2 LXTE;rNAL FLNCl'-L(N (ECIT) CSC1 PROGCRA.Mr CfiN CN,\STA,PZtRC, TZEfC C,hZEL RC, *CC2 1 (C;ANFAMNMUCPLSIGM(MARH L, *C02 2 I),P,C,ALPhA, *C02 E IECNSICF LZU( 1,) C.3 I NT E >i-N N, F PP I N 1,L INNEC r PtEPG *CO4 PRICGRit CA FNV N X,CLLX,.A,P,T,RHL,U,R, *C05 1 G-,CGFLG, RAD IUS,Nt CT, CELRAC *005 C I EN ES O X (5C 2), ELX ( 50), A(5CC), *A06 1 P (5CC),T ( SC ), HCO (500),U( ( CO), (500), *C 2 (C5C ),I:, CLG(5 C),RiADIUS(5EC),NCc r'(50c),O DC LrL (500t) *C06 ENTRY TO NLCLE2. *007 N UCL l1. ( 1T r C) 00 hE HAnE NC'. CNLY TC COMPUTE THE RADI1AL INCREEN IS AND THEI' CONTIRIFUTICN T THE CONDENSATE MASS INCREASE. TCIRC=FTSAT. (P(N)) *Co0 CEL C = SRT. ( 1. 379/ T( ().6027/U)*79 7 8.846*(TCROCP-T(iN) * l C C ELRA=CLX ( N ) *ALPHA/L*P ( N ) /ROL*ELRAC/U ( N ) * I=1 *(012 S1 CELRA;( I )=CELRAb *0C13 I=I+1 *G14 -wFENEVER I.L.N,TRANSFCR TO S1 *015 CCMPlTATION CF 1RESLLTANT TOTAL CONDENSATE MASS INCREASE CRoh=RAO IUS*R.'OC ILS*N.I*FELRAD *C 16 I=1 *017 S2 GRCk=CROW+RACIUS(I )*RADIUS(I)*NCCT(I)*CELKRA(I) *C81 I=I+1 019 FEN -E.ER I.L.N, TR4NSFER TC S2 *C20 OEL ( N ) =EL; ( N ) + 12. 56637*RHOL/MDCT*GRG *02 1 FUNCTIN R E T LRN *022 INTEGER I *023 ENC CF FLNCTIGCN *024 64

$CCMPILE MAC,PRINT CBJECT,PUNCH OBJECT VAPOR001 999066 06/24/64 11 10 33.5 AM MAD (12 MAR 1964 VERSICN) PROGRAM LISTING. AERONAUTICAL ENGINEERING NCZZLE CONDENSATION RESEARCH VAPOR PARAMETERS EXTERNAL FUNCTION (ARG) *cO1 ENTRY TO VAPOR. *02 TFE FCI LCWING READ DATA STATEMENT MUST SUPPLY THE QUANTITIES PRINTEC BELOW. REAC CATA (VAPER PARAMETER APPRCXIMATICNS AND TRIPLE POINT) *003 PRINT FORMAT VOCC,PTP,TTP *004 PRINT FORMAT VOL,PSAT,PSAT(1) *005 PRINT FCRMAT VOE,PSAT(2)...PSAT(4) *CO0 PRINT FORMAT VOF,SIGMA,SIGMA((1) *007 PRINT FORMAT VCG,RHOL,Rt1L( 1) *00C8 PRINT FORMAT VOH,L,L(1) *C09 FUNCIION RETURN *C10 CCMPLTE TEE LATENT HEAT -- ARG IS VAPOR TEMPERATURE ENTRY TO FL. *011 FUNCTION RETLRN L+L(1)*ARG *C12 CCVPUTE TFE TEMPERATURE CORRESPONOING TO PRESSURE ARG ENTRY TO FTSAT. *C13 TP=C *C 14 SS1 P=ELCC.(ARG) *015 WhENEVER P.LE.PTP *C16 TST=1:./(PSAT+PSTT( 1)*P) *G17 CTFER ISE * 018 TSAT=1. /(PSAT(2)+P* (PSAT(3)+P*PSAT(4))) *019 EN, CF CCNCITIOJNAL *0C2 TRANSFER TC SS2(TP) *021 SS2(C) FUNCTION RETLRN TSAT *022 SS2(1) FUNCTION RETLRN SIGMA+SIGMA(1)*TSAT *023 SS2(2) FUNCT I RETLRN RHOL+RHOL(1)*TSAT *0G24 CCMPUTE SURFACE TENSICN -- ARG IS VAPCOR PRESSURE ENTRY TO FSIGMA. *025 T P= 1 *C26 TRANSFER TC SS1 *C27 CCOPUTE LICLIC CENSITY -- ARG IS VAPCOR PRESSURE ENTRY TO( FRHCL *628 TP=2 *029 TRANSFER TO SS1 *030 CCMPLTE LOG CF PRESSURE CORRESPONDING TE TEMPERATURE ARG ENTRY TO LNPSAT. *031 T= 1. /aRG 032 WhENEVER T.CE. TTP *033 LNP= (-PSAT)/PSAT (1) *C34 SS3 FLNCTION RETLRN LNP 3*35 65

CR vHENEVER PSAT(4).NE.Ol. *36 REAL=-PSATM(31/2. *37 LISCRN=RCAL*REAL-PSA1(4)*(PSAT(2)-T) *038 WFENEVER CISCRr.L.C. *039 PRINT CEfvVENT $CSATURATICN CURVE APPRCXIPATIUN APOVE THE TRIP *C4U 1 LE POINT GIEVS IVACINARY LCG.$ *04C E RCR. *041 NENC CF C[NUITIC":AL *042 U ISCRf=S RT, (CI SCRM) 1 *43 L P( P R E A L - CER M P S A T ( 4) *44 kFAL=(REAL-jS Cs 8 /PSAT 1) (r 44 5 i0ENEVER LNP.L.PTP, LNP=REAL *046 TRANSFER TO SS3 *047 CTFE% (SE 48 LNP= (T-PSAr(2))/PSAT(3) *C49 TRANSFER TU S53 *050 ENE CF CONUITICNAL *0-1 INTECER 1,TAPLTP *05 2 UIEN IF'N SICIA(I),"H;L ( 1 L(1),PS AT( 4,LMTX(IC t)OVTX(19),Tp( *83 1 ) *053 VECICR VALLUF VlA=M 2F2C,Cl*$ C * 4 VLCTC(R VALUOESV(".;A=-H011'SAI'UKATICON~l CURVE EATA - NATURAL LOG CF P(DYN\IE/CM** *055 1 r )'J[aSUS 1/lW K)T *05 VEUTO)i ilLLtE VO v?=118,51 0j,06,2(1PE20.71*$ *0:56 F C T C if! V A L L F_ S V 1 =B I i- S PL f I G, C, 2 1 P E 2' O. 7 56 E~C5 VLEC lRk V fLU,-S V()C= $H*C'TR\I PLE. PC INT*, S4 tk(lIPE20. 7) *05 7 VOCTCR VALULS VU~ IL*ELINER APPRCXIVATICN TC I/T HELCW TRIPLE POINT -- *056 1 * S, 1P 1- I + L (P ), P E12.5 5 *058 98018? VAL LS 9f V SH CADR AT I C A PP RUX IPA T ION TOC 1/ I ABOV E T ILPLLE POINT T9 - * S 5 1-25, +LN(P)*C,1E12.5,H2U +LN(P) P.2*CiPE12*5*; *5 V ET(:R I:ALUL- V -F=$H,*L(NEAR APPROXI"ATICN TO SURFACE TENSION — *,S5, *060 t)P IE 1 2 H + ( I q C P TM P. I) 1 PE 1 2. * * I 0 61 V UTC TI PV A L U, E S V9%76 S I L I NOEA R A P P COX I TOAT I C n T LICUIL CONSITY —, 5, *01 I 1P 1E2 5, +II,(P TP FivP LPE112.5 *0,61 V EO TICK V'L~UF V01F $ H *C L I N EA k AP P CX ItIA T ION TO LATENT H-EAT — I,' 55, *A6 l 1 PF 12. -, +IVAPCR TE'MP ) *Q, I PE 1 2. 5*$ *362 ENC CF LUNUTICN 6263 66

APPENDIX B EXAMPLE INPUT 67

COMMENT CARDS FOR THIS RUN EXAMPLE OF AN UNDEREXPANDED SONIC FLOW IN A FREE JET. THE VAPOR IS NITROGEN. END SIGMA=23.94, -.1933, RHOL=1.181, -.0048, L=2.9E+099 -.0117E+09, PTP=1.173883E+01, TTP=1.5832805E-02, PSAT=3.00396E-02, -1.21024E-03, 3.20E-02, -1.3548E-03, -1.913E-06 COMMENT CARDS FOR THE FIRST PROBLEM THE OUTPUT FOR THIS PRObLEM IS TO BE REPRODUCED IN THE PROGRAM WRITE-UP AS AN ILLUSTRATION OF THE OUTPUT STRUCTURE. ALSO, A LISTING OF THE DATA DECK IS TO BE INCLUDED. DATA PZERO=6.895E+06, TZERO=298.. PE=3000. GAMMA=1.4, MU=28.02, CP=1.0364E+07, D=.65E-08, ALPHA=1. DELX=.05, XRANGE=5., MAXG=.1, XPOINT=5 MING=.00(1, MNDELG=.0005, MXDELG=*0015 MNDELX=.00001, MXDELX=.2,.2 * WEDGE=OB, ASTAR=*0507 INANG=20., OUTANG=39.57, XMIN=-10., XMAX=5. * EPSLON=1., DELT(1)=5., TRANGE=100* EPS=.1,.OC1, 0001,.OO001 68

APPENDIX C EXAMPLE OUTPUT 69

$EXECUTEFULL DUMPI/0 OUMP )9 923 7/164/64 $DATA 9 9 2.3 4 7/ 16/64 1 21 21.4 P4 MAP SCARDS 30600* SPEEK EODGC * SYSTEM 2306* ERROR 2062K* s i Pb 2 6 Si'NinT T~t:2* (M A IN ) 100CC FL 32353 FTSAT 32353 FS I MA 32323 FRHnL 32353 L!N;-S AT 32353 VAPOR 32353 IFLOW 33212 CO NON 33421 CFLCW 34D24 NUCLEI 35223 NLUCL_2 32~275o ISE N 35454 NOZZLE 36157.IOB 3o72 2* DFDP 43 o5* [FMP 43 65!*.E 432`4*.03311 43322*.P RINT 43337*.READ 4342 RD*.ROATA 432(2I*K.PC`3-MT 4467 S*0 44 h7 7.EXIT 44766* ELOG 4 5 C, C*.21 30 4 5 1622 EX P 452D>1* L E Ro 45363* C (U s 454 17* SIN 45417* FSPILL 45616* RSPILL 45616* (PROG) 45070 SU7i) 73676 (E R A JI 7774 1 32043 LOCS. CAN BE SAFELY USED IN EXPANDING FROG. (OCTAL) 26000 LOCS. CAN HE SAFELY USED IN EXPANDING FROG. (OCTAL) BEFCRE FULL CUR' LOADING PROCEDURE IS USED NOZZLE CONDENSATION RESEARCH PROGRAM COMMENT CARDS FOR THIS RUN EXAMPLE OF AN UNDEREXPANOED SONIC FLOW IN A FREE JET. THE VAPOR IS NITROGEN. TRIPLE POINT 11738830E 01 1.5832805E-02 LINEAR APPROXIMATION TO I/T BELOW TRIPLE POINT -- 3.00396E-02 +LN(P)*-1.21024E-03 QUADRATIC APPROXIMATION TO l/T ABOVE TRIPLE POINT -- 3.23000E-G2 +LN(P)*-1.3548-E —3 +LN(P).P.2*-1."913CE-f6 LINEAR APPROXIMATION TO SURFACE TENSION -- 2.3940GE D1 +(DROP TE —MP.)*-1.93320-."i-rE-C1LINEAR APPROXIMATION TO LIQUID DENSITY 1.181. 40E 26; +(DROP TEMP.)*-4.8BSIEO2E-23 LINEAR APPROXIMATION TO LATENT HEAT -- 2.9EOOOE 29 +(VAPOR TEMP.)*-1.1760GE 27 70

BEGINNING OF PROBLEM NUMBER I CCMMENT CARDS FOR THE FIRST PROBLEM THE OUTPUT FOR THIS PROBLEM IS TO HE REPRODUCED IN THE PROGRAM WRITE-UP AS AN ILLUSTRATION OF THE OUTPUT STRUCTURE. ALSO, A LISTING OF THE LATA DECK IS TO BE INCLUDED. COMPLETE LIST OF PROGRAM PARAMETERS FOR THIS PROBLEM. PZERO(DYNE/CM**2) 6.895E 26 TZFRU(DEG. K) 298.0000 RHOLERD(6M/CM**3) 7.798E-3 3 PE(DIYNL/CM**2) 3.00DE C3 GAMMA 1.400 MU(GM/GMOL) 28.02 CP(DYNE-CM/GM-DEG. K) 1.036E 17 INTERMOLECULAR DISTANCE 6.500E-09 ACCOMODATION COEFFICIENT 1.00u DELX.050C MIN(DELX).CCMC MAX(OELX) IF NiOT=C.2000 MAX(DELX).2000 WILL BEGIN CHANGING'DELX' TO KEEP DELG GREATER THAN.0500 PERCENT AND LESS THAN.1500 PERCENT WHEN, G EXCEEDS.C1uC PERCENT XRANGE IS 5.uC- HUT'WILL TERMINATE IF G EXCEEDS 10.20 PERCENT. A CONICAL NOZZLE WITH THE FOLLOWING SPECS. IS BEING USED -- THROAT AREA.05070 INPUT HALF-ANGLEIDEG.) 20.00GGO INPUT LENGTH -10.0000 OUTPUT HALF-ANGLE(OEG.) 39.57JO0 OUTPUT LENGTH 5.0002 SATURATION POINT SPECIFICATIONS -- PRESSURE 1.881E 04 TEMP. 5.516E 01 DENSITY 1.149E-64 MACH NUMBER 4.692 VELOCITY(CM/SEC) 70947.86 NOLZLF PO1NT.5242 A iR EA.986C9 TOT-,L MASS FL4()M-CNl/SEC) 8.04159 DETERMINE A CO')NDENSATiO!,N ON THE OASIS`OF THE VALUE OF EPSILON(PT), WHICH EVER APPLIES. SUFFICIENCT IE I T NITIAL C OL E NILE A NSSAT ION\ F IN1I TIATE THE CALCULATI N S. CU)NbE NSAIIfIJ iN POINT FO0 LLUWS. X=.9257 A= 2.4993) P=.48695L 04 RHO=.43775E-C4 T= 37.4901 U=?s483.57 P= 5.8944 RADIUS= 4.10282E-C8-S \JCEOT= 1.29-92E 17 G= 4.3o6682L-06 71

X PH AT THAT RHHAT1 PiERC'HAI' H MI U u.92572 1. C OC165 l.rCOCOD46, 1 00125 1.G000901i4 5.89439 5.89439 73483.665.0C043668 1.22572 1.11637212 1. C4835926 1. 05232979 1.06049494 6.43497 6.60377 74347.828.3363U150 1.35072 1.17677504 1.10211183 1. 917951 1.062594(3 6.52466 6.86909 74598.18'.7974188 1.47572 1.24937637 1.16972242 1.06505322 1.06231128 6.50156 7.12033 74828.747 1.34568967 1.60072 1.32770474 1.24422354 1.07009847 1.06092466 6.57306 7.35937 74994.598 1.916o9C53 1.72572 1.40912810 1.32286535 1.07440567 1.05889918 6.57061 7.58760 75162.813 2.480'7035 1.85072 1.49146686 1.40331079 1.07829214 1.05659957 6.56237 7.80624 75318.182 3.02480909 1.97572 1.57300285 1.48366651 1.08126585 1.054283J8 6.55373 8.01632 75463.894 3.54369906 2.10072 1.65311074 1.563256C6 1.08399963 1.05200428 6.54629 8.21866 75601.555 4.r,3296441 2.22572 1.73123863 1.64142771 1.08635 16 1.04982834 6.54128 8.41397 75132.398 4.49093986 2.35072 1.80707572 1. 71777581 1.08839522 1.04779321 6.53920 8.60286 75857.274 4.91764873 2.47572 1.88039872 1.79199050 1.09016861 1.04592089 6.540i32 8.78587 75976.842 5.31425929 2.60072 1.95110263 1.86392486 1.09170452 1.04419923 6.54455 8.96346 76C91.563 5.68261015 2.72572 2.0193325C 1.93369046 1.09323683 1.0426C6o9 6.55139 9.13601 76201.710 6.02482647 2.85072 2.08551580 2.C00169969 1.09420782 1.04110962 6.55993 9.30393 76307.410 6.34296513 2.97572 2.14929032 2.66749055 1.09523097 1.03974678 6.57081 9.46751 76409.187 6.63914955 3.10072 2.21079704 2.13117468 1.09612484 1.03850290 6.58366 9.62703 76507.231 6.91541952 3.22572 2.27013130 2.19287246 1.09690614 1.03736521 6.59817 9.78276 76601.724 7.17364043 3.47572 2.38762268 2.3C990592 1.10038704 1.03767382 6.63233 10.08373 76779.726 7.63846189 3.72572 2.49839348 2.42076355 1.10330294 1.03799909 6.66970 10.37202 76945.393 8.05122972 3.97572 2.60334232 2.52622792 1.10576496 1,03832151 6.70912 10.64895 77099.935 8.42032981 4.22572 2.70310935 2.62686610 1.10785832 1.03863707 6.74988 10.91565 77244.466 8.75251126 4.47572 2.79834184 2.72325906 1.10964979 1.03894433 6.79133 11.17305 77379.929 9.05321681 4.72572 2.88977814 2.81615660 1.11119476 1.03921627 6.83271 11.42198 77507.104 9.32682598 NOZZLE SU8ROUTINE ERRCR INCICATION - MESH POINT DOES NOT LIE WITHIN THE NOZZLE. FINAL OUTPUT IN TABLE I FOR PROBLEM 1 4.92572 2.96298575 2.88817078 1.11312878 1.04029332 6.86559 11.61550 77602.939 9.52749002 72

A~~~AA ROAR.43fSE3 TBAR.24~99909E 01 PBAR,.43775IOE-0.561323E-0248E 01 p.379010 02 125go57E 00.3570823E-02.408257 ~ 8~ ~ ~.026E-3.084E0.2,784482E-04.4655958E 01.92572 486951E 04.213461 02 02E 04 -3009831E-rJ21253E0.sz57z.sssb~~389476E-03.15977E iO.2347.5696414E 01E 3.2686620E 04.3146813E 02 01055 2571434E-(91SE 0 1.22572.j2C~~0324018E-03 v.2005100.21E-04 9.1941 o g.2209171E 04,31050037E; 72804-0n.6646 E 3 1~35072 2~~~~~0207815E-03.3219.105157D0.1732804E 222215.7578463E 0.452.186703BE 04.313-3679E 02.11 56-04 0193931SE 02 I 052 2330lgBE-03.1057080E 00.1512256-ijZ.8620055E 01 1.6072.106664E 04.3150oq9E 02 ~31150E-04.1707C68E- 0z 1.60072 ~ ~~~~.203246SE-03.16123E 00.133'9128694E 015E C.752 1401386E 04 7108ia 7E 02.1064 bE-5140090OE-i 0 I *j7.790798E-03.31.1180606E3510604 043395E 025 I.502.1234755E 04.391766E 02.1071621 2.35106-0.109fLE0 1.85072 ~ ~~~~.1590026E-03.31 105412 1E242204 14711E 02E 9 1.9572.196323E 04.3o1705E 02.1214221E-02 3568 E 0 1.97572 ~ ~~~~~14?-0827E-03 320828.946 8344E2405 456228E 02 0.9796600E 045E30.32.71 0 8550551E-05,431 E0:::020.80197E03 1276.C-I 1087074E 00.9562-3.1483370E 02 2.22572 127E 03.3239480E 02.7759376E-05 *9950632E -0.7943914E 03.115.1090130E 00.q` 3.1627756E 02E 0 2.35072 119E 03 *3248589E 02 70-7 2606E-1)5.,9069916L 02 2.35072~ ~ 7 0.1044 io03.192013E 00.8 051778841E 029E 3 2.457 1719919 E 04 03254198E 02.64727OIE-05 #83(0599E 03 2.4757 6485C 0.9496809E-0.1092891E SG 4Z't1936643E 025 E 3 2.60072.866701E-04 3256815E 02.594565-le-05.641E0.2101143E 02 2.252 5975946E 03.3257353E 02.1093071E 05 07027750E 03 BE 03 0793533E-0.1092365E 30.548015.2 27234SE E Z 2.807.47169 103 325524SE 02.56 87-05.64979JOE 03 2~85072 1~~~~~0286685E-04 050669 0.0267E 2 55438~~ 2.752 5024169E 03.3250879E 02.109089E 00.02 S5ZE203 945 - 57 2.9~~~~5~~2 ~~6708541E 0.08789E Ui0 *.2826189E 0878 E 0 02.4625539E 03.3244590E 02.088 t0F-05 05232953E-U3 3-1007.6191794E-04.1-8617E 013.408059'322 894GE 0; 3 3.Z52.4269242E 03.5266-4.32261,99E 02.3576472E-05.4586471E ti.3669289E 03.2324E 02.107542E U.402350E03 *35850E 0; 3j.47572.4612733c-04.3?0390.664Z-3 9623E~ 2572.3180479E 0.3179117E 02 281L.3606 042810-E 3 3.7 3. 4028961E-04.057929E 00.2829946E-05,4598101E f, 3.752.2777969E 03.3152629E 02.1 1203E-05 *3229346E D3::9~ ~.4422E 3 3543471E-04.1048722E OD.518908538E-03 5014 4.22572 2396 03 03319-4.3'125192E 02 0341 0.22-68041 E-05 4.752.216 —6E0.3097611E 02 4.47575229 0.279192(E-04 4.72572.9 r) r_ A THE U~~r\ i~ F-~~ ~s ~ $t~~'~ S~A~_6b I~Is~aI::: k'~ ~ ~ ~ ~~7

X O EL X O ELG RADIUS NOO r RHOORP MEAN RADU.9 2 572.050-0C, 303.03343 66 8.330436'68.41G282lE —,.7.1293921E 18.73'-435660 12.413'2 82E7 1.2 2 572.05 CC0003"I3D.33630153.12 91 72 35.3GCb707 07-u7.5587431E 20~.5D74788E 15.349639-7 1. 0 7 2.2 503333C`0,k.79741886.13-2226842.2958053E0-7.25187610E 20.83198124E 15.373423-7 1.47 572.025COOC", 1.34563967.11275868. 3000i 62 3 6 E-7.1375553E 20.90.298850'- 15.43"-4426-7 1.60572.3~J2 5003%O 1.9 1613 90',5 3.11416137. 3 i,8 4 3 36W1.5671292E 19.8643484E 15.438714-7 1.7 257 2 C"2 5C C L 248 C: O7 l335.11154598.3193885E-J"1.21087420 19.780-,9273E 15.4 725 22-7 1.35072.302 5C0033 3. 248 959.170733834.3 3 16 2 6(1E-37.6782238Z: 18.6938583E 15. 503553E7 1.9 7572. 02 505033) 3.5430`991o'.10-149 754.34412020:-01.25"74347E 18.616C509`0 15.530957E7 2.10012.5~~ID2 5 2C C.23295441.09537515.35649%3E-C7.6181483E 17.5492370~E is. 5547 89-3 2.2 2 572 C 2 t5C 0%02 4.4909-)3986.0893C-7345.36836960-5J7.1853,8590-' 17.4923669C 15. 575496E7 2.3 5072.025C0C00- 4.'917648'73.0C8282397.37946400-07.5714726L 16.44375980 is'.59 3 5S7-0 2.4 757 2. ij2 5 03303- 5.331425929.2-7699641.38955670-0"7.891016.43195136 15. 6C946603 2.600O"724. C32 5 C 3 5.6 82 )I ~1(15.U7150978.3 98 5301-0-'ClE-'7.65595561: 1 5.36574900~L 15.6 23414-7 2.7 2 572. 62 s0%03 6. 2426247.36 64 52 43.40'64288E-211;7.24929840- 15.33427%80 15.6356876-0 2.85,072. J02 5 C33C 6. 342965 13.6 18;_23 32.4134879E-0-7.iC'96373- 1 5.30)6 54 570 1 5.647029-0 2.9 7 572. ~ —2 5 C36.63914955.0-5758389.41942110-z07'.44979570- 1 4.28216990 15.657 69-7 3.1,~J7 2.0i)2 53n030 6.91I541952.5,37 53 C01.4 24 3 0'2 2 E-,7.22125 56 1 4.26357s2l 10.6661574-0 3.2 2 572I.2 53003 7.113643-43.0s 2 8.37 6.42S2174E-U81.118 1359 14.241 35816 15).674426-7 3.4 1572.3s32 7.63846189.088490682.43281691c —~7.94 53"I4~ 1 3.27~876d20 15.668883207 3.7 2 572.351000 c 8.3512297k.0-7S61711. 4 34 89.3 E- 7. 4 942 8 2 1 3.182331kt3L 15. 7Cl1 1160 3.97 572 6 8. 4 20-i32''_981.-73-6-71 30. 4 3 5 4 7,8-31E. 3 2 1706 33 1 3.10631 i 11 17 3107 4.2 2 572 5, > 3, 8.752!81126.63763-64.433 7197E-H,~. 2 4,9 44 1 3 142-5233E 15 116L7 4.4 7 572 5 0,, I- O'- 9. 532 i68 1 C.05784621.4 313 149-.22 -391f0 1 3.12732260 -15. 72 9 1914.7 2 572 C 3b 29 8.5 27313.42626450-.1'.1 4,-''4 bL 1 3.11442193 1. 136472-7 ALL INPUT OATA HIVE i 309 2 E400005. ST 00 4,3552 74

JOB NO. 999234 MONDAY, JULY 06, 1964 1 27 17.5 PM 1 27 54.8 PM SIVIER (AERO. ENG.) S418F 0O2 200 C O 999234SIVIERS418FC070664484375 55 15 X 0 C 25 14 $ $ $ $ $ $ $ $ $ $ $ $ $ $ $ $: a $ CARDS PROCESSED 264, LINES PRINTED 144, PAGES PRINTED 7, CARDS PUNCHED ELAPSED TIME -- TOTAL PROCESSING EXECUTION LOAD I NG O' 33.0'' 0' 8.4'' C' 9.5'' O' 15.2'' SEE NOTICE ON BULLETIN BOARD DATED 6/8/64 CONCERNING REPLACEMENT OF MANY LINEAR SYSTEM SUBROUTIN.ES (6JRDT L-T AL). 75

Unclassified Security Classification DOCUMENT CONTROL-'DATA.- R&D (Security classification of title, body of abstract and indexing annotation must be entered when the overall report is classified) 1.- ORIGINATING ACTIVITY (Corporate author)- 2a. REPORT SECURITY C LASSIFICATION Aerospace Research Laboratories Unclassified Thermomechanics Research Laboratory 2b. GROUP 3. REPORT TITLE A Digital Computer Program for Condensation in Expanding One-Component Flows 4. DESCRIPTIVE NOTES (Type of report and inclusive dates) TDR ~ c ~ ~ A~s S. AUTHOR(S) (Last name, first name, initial) Leonard J. Harding 6. REPO RT DATE 7a. TOTAL NO. OF PAGES 7b. NO. OF REFS March 1965 73 3 8a. CONTRACT OR GRANT NO. AF 33(657) 8867 9a. ORIGINATOJR'S REPORT NUMBER(S) b. PROJECT NO. 7116 ARL 65-58 c. Task No 7116-0Ol 9b. OTHER REPORT NO(S) (Any other numbers that may be assigned this report) d. 1 0. A VA IL ABILITY/LIMITATION NOTICES 11. SUPPLEMENTARY NOTES 12. SPONSORING MILITARY ACTIVITY USAF 13. ABSTRACT This report describes a digital computer program for calculating vapor condensation processes that occur'in rapidly expanding flows, The treatment emphasizes the program logic required to satisfy the requirements of the mathematical model. Among the most important of these requirements are the search for the onset of nucleation,'in an isentroplically expanding flow, and the iteration procedure necessary for the joint solution of the nucleation and growth equations and the diabatic flow equations. The program features flexibility in allowing input conditions, short execution time, and a convenient output format.

Security Classification 14. LINK A LINK B LINK C KEY WORDS____ _ _ ROLE WT ROLE W.T ROLE W T T-~~~~~~~~~~''TV I, i r~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ i i I~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ I~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ I_ _ _ _ _ ___ _.. _I INSTRUCTIONS 1. ORIG!NATING ACTIVITY: Enter the name and address imposed by security classification, using standard statements of the contractor, subcontractor, grantee, Department of De- I such as: fense activity or otiher organization (corporate author) issuing (1) "Qualified requesters may obtain copies of this the report. | ~~~~~~~~~~(l) "'Qualified requesters may obtain copies of this the rt-:port. report from DDC." 2a. NEPORT SECURITY CLASSIFICATION: Enter the over-..REPORT SECUITY CLASSIFICATION: Enter the over-.(2) "Foreign announcement and dissemination of this ll security classification of the report. Indicate whether "Restricted Data" is included. Marking is to be in accord- ance with appropriate security regulations. (3) "U. S. Government agencies may obtain copies of this report directly from DDC. Other qualified DDC 2b. GROUP: Automatic downgrading is specified in DoD Di- serh reqe rou users shall request through rective 5200.10 and Armed Forces Industrial Manual. Enter the group number. Also, when applicable, show that optional markings have been used for Group 3 and Group 4 as author-. ized(4) "U. S. military agencies may obtain copies of this zed. report directly from DDC. Other qualified users 3. REPORT TITLE: Enter the complete report title in all shall request through capital letters. Titles in all cases should be unclassified., If a meaningful title cannot be selected without classification, show title classification in all capitals in parenthesis (5) "All distribution of this report is controlled. Qualimmediately following the title, ified DDC users shall request through 4. DESCRIPTIVE NOTES: If appropriate, enter the type of. report, e.g., interim, progress, summary, annual, or final. If the report has been furnished to the Office of Technical Give the inclusive dates when a specific reporting period is Services, Department of Commerce, for sale to the public, indicovered. cate this fact and enter the price, if known. 5. AUTHOR(S): Enter the name(s) of author(s) as shown on 11. SUPPLEMENTARY NOTES: Use for additional explanaor in the report. Entel tast name, first name, middle initial. tory notes. If military, show rank and branch of service. The name of the principal a.xthor is an absolute minimum requirement. 12. SPONSORING MILITARY ACTIVITY: Enter the name of the departmental project office or laboratory sponsoring (pay6. REPORT DAT-_- Enter the date of the report as day, 6. REPORT DAT.. Enter the date of the report as day, ing for) the research and development. Include address. month, year; or month, year. If more than one date appears on the report, use date of publication. 13. ABSTRACT. Enter an abstract giving a brief and factual summary of the documernt indicative of the report, even though 7a. TOTAL NUMBER OF PAGES: The total page count it may also appear elsewhere in the body of the technical reshould follow normal pagination procedures, i.e., enter the port. If additional space is required, a continuation sheet shall number of pages containing information. be attached. 7b. NUMBER OF REFERENCES: Enter the total number of { 7. NUMBER OF REFERENCES Enter the total number of It is highly desirable that the abstract of classified reports references cited in the report. references cited in the report. be unclassified. Each paragraph of the abstract shall end with 8a. CONTRACT OR GRANT NUMBER: If appropriate, enter an indication of the military security classification of the inthe applicable number of the contract or grant under which formation in the paragraph, represented as (TS), (S), (C), or (U). the report was written. | There is no limitation on the length of the abstract. How8b, 8c, & 8d. PROJECT NUMBER: Enter the appropriate ever, the suggested length is from 150 to 225 words. military department identification, such as project number, 14. KEY WORDS: Key words are technically meaningful terms subproject number, system numbers, task number, etc. subpojec numer yste numers asknumbr et. l14. KEY WORDS' Key words are technically meaningful terms subproject numbe' s m n s or short phrases that characterize a report and may be used as 9a. ORIGINATOR'S REPORT NUMBER(S): Enter the offi- index entries for cataloging the report. Key words must be cial report number by which the document will be identified selected so that no security classification is required. Identiand controlled by the originating activity. This number must fiers, such as equipment model designation, trade name, military be unique to this report. project code name, geographic location, may be used as key 9b. OTHER REPORT NUMBER(S): If the report has been l words but will be followed by an indication of technical conassigned any other repcrt numbers (either by the originator | text. The assignment of links, rules, and weights is optional. or by the sponsor), also enter this number(s). 10. AVAILABILITY/LIMITATION NOTICES: Enter any lim- I itations on further dissemination of the report, other than those Security Classification

UNIVERSITY OF MICHIGAN ~~~~~~~ I I/II I II 3 9015 02652 503I 91 5I 10265 1503 3 9015 02652 5033