FORTRAN IV PROGRAM FOR ANALYSIS OF SATURN V PROPELLANT SYSTEMS E. B. Wylie NATIONAL AERONAUTICS AND SPACE ADMINISTRATION Contract No. NAS 8-20312 June, 1967

This report was prepared by The University of Michigan College of Engineering Department of Civil Engineering under Contract Noo NAS 8-20312 PROPELLANT LINE DYNAMICS for the George C. Marshall Space Flight Center of the National Aeronautics and Space Administratio-n. The work was administered under the technical direction of the Propulsion and Vehicle Engineering Laboratory, George Co Marshall Space Flight Center with John Ro Admire acting as project manager. ii

TABLE OF CONTENTS Page LIST OF FIGURES......................... aaI00........... iv I INTRODUCTION............. OOOO....................... o 1 II DESCRIPTION OF MISSILE SYSTEM.......................... 2 III BASIC THEORY............................................. 4 IV PROGRAM ORGANIZATION AND DETAILS....................oo 11 A. Program Options..........................o........ 11 B. Program Organization................................. 13 C. Initial Steady State Conditions................... 14 V BOUNDARY CONDITIONS...............o.............. oo..... 17 A. Series Pipe Connection..............o...............o 17 B. System Excitation at the Feed Tanks.o..............o. 19 C. Turbopump Representation............................. 19 D. System Excitation at the Pump... o o..oo.............. 23 E. Pulser............................. o........o... 24 F. Pressure Volume Compensator..... o.............OO....O 26 G. Feed Tank Surface............000...........000.....o 26 H. Engine for Combined System Analysis.....O............ 26 I. Engine When Systems are Isolated at the Combustion Chamber. o................................. oo........ 31 Jo System Excitation at the Burner Plate Orifice....oo o 32 VI THE COMPUTER PROGRAM o..... o...........oooooo..o..oo... o 33 VII SAMPLE PROBLEM........... o............................. 4 VIII SUMMARY............00..00000.000......Ooooooooooooooooooo. 45 APPENDIX A o o o... o... o o o o o 0 0 0 o... o o. o o o o o o o. o o o o o o o o o o 0 o o 0 ~ 0 o 0 0 5246 REFERENCES. o...o.....................,.......... oo o.o....... 52 iii

LIST OF FIGURES Figure Page 1. Schematic Diagram of Saturn V....o..o..... o o........ 3 2. Control Volume.........................e o o o o o o o........o o o o o o 5 3. Characteristics in xt Plane...oo........... o.......o 5 4. Characteristics at the Boundary Conditions.............. 5 5. Series Connection.....,.......o.......... o o o o 18 6. Connection at Feed Tank.............,.................. 18 7. Pump Characteristic Curve. Pressure Rise vso Suction Pressure...................o o............................. 22 8. Pump Characteristic Curve. Pressure Rise vs. Discharge...............................................o o o o o o o o o o o o o o o o o o o o 22 9. Location of Pulser..............o.......... o 26 10. Schematic Diagram of Engine Combustion Chamber, Burner Plate Orifices, and Pump Discharge Lines.,,,o........... 26 11. Listing of Saturn V Program................,o......... 35 12. Sketch of System Showing Location of Program Output Quantities P, Q o...... o................ o o o.. o o... o o o o o o o o o o 40 13. Input Data for Example Problem.......................... 42 14. Output from Analysis with Oscillating Pump....o o.. o,, 43 15. Pressure at Pump Inlet of Oxidizer System of Titan II Due to Pump Oscillation, o = 62.8 rad/sec....... o..... 44 iv

I. INTRODUCTION The digital computer program to be used to analyze the propellant system dynamics of Saturn V is presented in this report, The basic theory underlying the method of analysis is developed and the program details are described, The objective of the program is to provide an analysis of the transient flow condition that may develop in the propellant feed system of a liquid fueled rocket when subjected to longitudinal structural vibrations. The details described herein refer to Stage I, Saturn V, but the program is equally applicable to the other stages of the missile, or to other liquid-propelled missiles of similar configuration. In this report the missile propellant system is first described, the basic theory behind the method of analysis is set forth, then the program is described and details are presented. The report is concluded with the results from an example problem, -1

II, DESCRIPTION OF MISSILE SYSTEM A schematic view of the propellant feed system to one engine of Stage I, Saturn V, is shown in Figure 1o Only the component parts considered to be an important influence on the fluid system response are showno The turbine gas generator and turbine are assumed to have no affect on an oscillation in the feed system. Stage I has five engines, each with separate supply systems leading from the lox and fuel tanks. These are assumed to act independently in parallel without any fluid flow or pressure interaction at the feed tanks, Any differences that exist between the outboard and inboard systems can be treated by using the data appropriate to the system being considered. Figure 1 shows the component parts that are believed to be important in the dynamic response study, namely, the feed tank, the suction lines, the PVC, the pump, the parallel discharge lines, and the engiie which includes the injector plate and combustion chamber. The fuel and lox systems are separate units with the combustion chamber beyond the orifice injector plate being the only point in the fluid flow and pressure field that is common to both systems, With this configuration a meaningful analysis can be performed in either system alone, assuming the other steady, or in the two systems together wherein an interaction is permitted at the combustion chamber, If a singlefueled-engine missile were being analyzed only one system would be considered,

-3 1 TANK LOX I I FUEL ~~~~~I I TANK 1 LOX SYSTEM K= SYSTM = 1 FUEL SYSTEM K= 2 SYSTM = 2 COMBINED SYSTEMS SYSTM = 3 I Pulser or PVC Lox Pump Fuel Pump Lox Dome L I l Fuel Monifold Burner Plote Combustion Chamber Figure 1. Schematic Diagram of Saturn V

III. BASIC THEORY The equation of motion and continuity equation that describe one-dimensional unsteady flow of a compressible fluid in an elastic pipeline are first developed. The equations involve the actual distributed parameters that describe the system and include the nonlinear viscous terms. These partial differential equations are transformed into particular total differential equations by use of the method of characteristics. A finite difference method is then used to place the total differential equations in a form suitable for numerical solution on a digital computer. The transient conditions in each pipeline in the system are treated in accordance with the above descriptiono Boundary conditions are used to transfer the effect of one pipeline to another, and to introduce the effect of the pump, the PVC, the engine, etc. In a typical computer analysis of a specific problem, an initial steady-state operating condition is imposed upon the system; then a numerical solution is obtained at finite points throughout the entire system at finite time increments during the life of the transient. The equation of motion is developed by examining a control volume that is moving with the vehicle, Figure 2. p A Ax + T iD LAx - pg Ax A sin a + pA Ax d - 0 x 0 dt The density of the fluid is given by p and the pressure is given by po The angle a is measured between a plane normal to the flight direction -4

iFLIGHT DIRECTION X-S AREA A t oA rD g pA Ax,, x p to at R S I I -- I -_ I \In x Figure 2. Control Volume Figure 3. Characteristics on the XT Plane p p S x R Figure 4. Characteristics at the Boundary Conditions

-6 and the positive x axis, and g is the acceleration at the particular time in flight. Other variables are defined in Figure 2. The equation (1) can be simplified to the following form: Px fV2 - + -+ VVX + Vt - g sin a = 0 p 2D The subscripts x and t indicate partial differentiation. The condition of continuity applied to the same control volume at any instant yields: (pAV)x Ax + (pA Ax)t = 0 which can be simplified to x + 1 dA + 1 dp = O A dt p dt The introduction of the elasticity of the pipe walls and compressibility of the fluid allows the equation to be rewritten in the form: 2 VPx Pt a Vx + -- + - = 0 P P The pressure pulse wave speed in the pipeline is denoted by a, and, for liquid flow in an elastic tube, is defined by a = /p 1+_ Ee The bulk modulus of elasticity of the liquid is given by K, the modulus of elasticity of the pipe wall by E, and the pipe wall thickness by e. The equation of motion and continuity equation can be combined

-7 in a linear manner using an unknown multiplier, X. Thus, by multiplying the continuity equation by X, adding the simplified equation of motion, and rearranging, we have [ Px (V + + Pt + Vx (V + a 2) + Vt + - g sin a = 0 P 2D If dx is restricted to the values V + 1 and V + a2, then this dt X equation can be written as a total differential equation. X dp dV fV - +- + + - - g sin a = O p dt dt 2D The restriction imposed upon d requires that X = + When these dt -- a values of X are substituted, the following set of equations results dp dV fV (1) L- d + dV + f _ g sin a = 0 (1) ap dt dt 2D C+ dx = V + a J (2) dt 1 dp dV fV2 - 1 dp + dV + fV g sin a = (3) ap dt dt 2D X Cdt=V -a J (4) These are the characteristic equations. The first equation is valid only if the second is satisfied, and the third equation is valid only if the last is satisfied. By use of the first order finite difference approximation to these equations, they can be placed in a form which is suitable for numerical solution on the digital computer. The independent variable (xt) plane is shown in Figure 3. If conditions (p,V) are considered known at points R and S, then they

-8 can be determined at point P at time At later. In Eqs. (2) and (4), the velocity V is generally much less than the wave speed a, and can be neglected. Since a is a constant in any given pipeline, Eqs. (2) and (4) are straight lines on the xt-plane, RP and SP, respectively. It is also convenient to write the equations in terms of mass flow rate (Q = pAV, slugs/sec.). In finite difference form, Eqs. (1) to (4) become QP -QR + a (PP - PR) + 2t gpA At sin a = o (5) xp - xR = a (tp - tR) (6) A, \PS f At 2 QP Q PS pA - ( P) + 2D gpA At sin a = 0 (7) xp - XS = -a (tp - tS) (8) These four equations enable a solution to be obtained for the four unknowns, Qp, pp. Xp, and tp, if conditions are known at R and So An orderly computer solution is obtained if each pipe length is divided into N equal reaches, Ax = L/No This specifies the time increment that can be used since for stability and convergence of the solution it is necessary that At a < Ax in each pipeo A specified time interval grid is therefore established in each pipeline, Figure 3. If conditions are known along the system at the initial time, to, they can be determined at time At later by a systematic application of Eqs. (5) and (7) at each point P. Simultaneous solution of Eqs. (5) and (7) yields = [R + Q + (PR - PS) - DA (2 + QS2 + gpA At sin a (9)

-9 a af At PP = [PR + PS + A (QR - ) -A22 (R2 - S2] (10) Conditions can be found at points R and S by use of a linear interpolation between known conditions at points on the specified time interval grid A, C, and B. Equations (9) and (10) permit the calculation of all interior points of the grid for all time, however, the influence of the pipe end conditions must be introduced to yield the complete solution. At the inflow end (left side), Fig. 4, of the pipeline, Eq. (7) can be written in the form A fAt 2 A = Q - A PS - fAp 2 + gpA At sin a + A() or QP = Cl + C2 pp (12) where C1 is a variable that can be evaluated in terms of known quantities at each new time step, and C2 = A/a. Similarly, at the outflow end (right end) of the pipeline, Eqo (5) can be written ApR f At 2 A P = QR + QR + +gpA At sin a - pp (13) a 2DAp a or p - C3 - C2 pp (14) where C3 can be evaluated in terms of known values. When an external condition that relates Qp and pp, or defines one of the variables, is applied to a pipe end a solution can be obtained in combination with either Eq. (12) or (14) for the inflow or outflow end, respectively.

-10By combining the entire system in an orderly manner with the appropriate boundary conditions, a complete solution is obtained for points throughout the entire system for all time during the life of the transient condition.

IV. PROGRAM ORGANIZATION AND DETAILS Figure 1 displays a schematic diagram of the entire propellant system from feed tanks to combustion chamber. The component parts of the systems have been numbered to enable simple identification in the program. The results from any given execution of the program show the complete system response to a particular excitation. That is, an initial steady state operating condition is defined for the system; then, the transient condition produced by the exciter is followed for as long as may be desiredo A summarized discussion of the options in the program are presented below, followed by an itemized account of the major steps in the program. The details involved in the establishment of an initial steady operating condition are then presented. In general, the order of discussion of the program details follows the order of computation in the actual program. Ao Program Options The following major options, pertaining mainly to the system to be analyzed and the type of exciter to be used, must be decided upon by the program user and specified in the data. (a) System to be analyzed. Both lox and fuel systems can be analyzed in the same run, in which case the variable SYSTM = 3. Excitation can be in either or both -11 -

-12 of the systems. If the lox system only is to be analyzed, the variable SYSTM = 1; and if the fuel system only is to be analyzed, the variable SYSTM = 2. (b) Combustion chamber. If one system is excited, it may be desirable to isolate that system from the other in order to study its response alone. This is accomplished by setting the variable SINGLE = 1, which isolates the systems at the combustion chamber and essentially assumes steady-state flow from the second system. The more common situation, wherein one system may influence the other through the combustion chamber, is obtained by the variable SINGLE = 0. By using this option, one system can be excited, and the extent to which it influences the other system can be studied. Or, both systems can be excited and the resulting oscillations in the combustion chamber can be observed. (c) Pressure volume compensator. The PVC is included in the fuel and lox systems if the variable PVC = 1. If PVC = 0, the PVC is not in the system. Thus the effectiveness of the PVC for any given excitation can be evaluated. (d) Points of excitation. There are three possible excitation points in each system and one disturbance mechanism which is common to both fuel and lox systemso The desired amplitude and phase angle can be given for each of these elements at the frequency OMEGA. The table summarizes the variables involved. The variables pertaining to a particular location must be set equal to one to use the exciter and an amplitude must be giveno If the exciter is wanted in only one of the systems, the variable controlling

-13 the amplitude in the second system must be set to zero. In order to remove any exciter from the system the variable controlling that element is set to zero. Lox Lox Fuel Fuel Location Variable Amplitude Phase Amplitude Phase Junction of supply tank to suction pipeline TANK = 1 AMl(1) ALPHA(1) AM1(2) ALPHA(2) Pulser in suction line PULSER = 1 AM5(l) ALPHA5 AM5(2) ALPHA5 Pump PUMP = 1 AM7(1) ALPHA7 AM7(2) ALPHA7 Burner Plate Orifice ORIFIC = 1 AM9 ALPHA9 AM9 ALPHA9 B. Program Organization The computing procedure followed in solving a particular problem involving a particular set of data is as follows: a) Read in the data pertaining to the system. b) Calculate steady-state conditions throughout the system and initialize the various variables and constants required for the transient analysiso c) Set up the various column headings for the print out of the results of the calculations and print out the initial steady-state conditionso d) Increment the time variable, T, by At and the integer counter, U, by one. If the lox system is to be analyzed, proceed. If the fuel system only is to be analyzed, the variable K is set equal to 2, e) Calculate the interior points of all the pipelines in the system being analyzed.

f) Calculate the boundary conditions of each pipe in the system. This includes all the series pipeline connections, the pump, the pulser or PVC, and the feed tank. g) If the lox system was analyzed and the fuel system is going to be analyzed, the program returns to step (e) and repeats the steps for the fuel system. If the fuel system is not to be analyzed the program continues with the next step. h) Calculate the engine boundary condition. i) Substitute all newly computed values of pressure and discharge to a second storage location so they can be reused during the next time increment. j) Compute the pump characteristics and constants to be used for the computations of the next time increment. k) If a print out is desired, transfer to the print statement, otherwise, return to step (d) and repeat the calculations at the new time. A check is made to see if the time exceeds the desired duration of the transient calculationso If so, the program is stopped. C. Initial Steady State Conditions Data are given for the number of reaches in one particular pipeline, as well as its length and pressure pulse wave speedo These data are used to calculate the time increment, DT, to be used for the transient calculation. The number of reaches in each of the remaining pipes in the system are then calculated, based on this DTo A number

-15 of other constants that are needed in later calculations are then established and stored for future use. At any location in the system where identical parallel pipes exist, they are treated as a single pipeline having an area equal to the combined area of the pipes but the actual pipeline diameter is maintained. An example of such a situation would be the two pipes labeled number eight in the pump discharge section of the lox system, Fig. lo With this representation the velocities are correctly represented in the pipe and the frictional losses are also correctly described. This substitution provides a significant simplification in the program. The given data, in addition to the physical description of the system, that define the initial operating conditions are: the flow in each system, the pressure in each feed tank, and the pressure in the combustion chamber. The factor DPF, which corresponds to the change in hydraulic gradeline in each reach of each pipeline, is computed; then the pressure at the pump outlet is computed by beginning at the combustion chamber and working through the orifice and discharge piping (see Figure 1)o All steady state operating conditions (Q,p) are now stored for each point at which transient computations will be performed in the system. The pressures in the suction side are computed by beginning at the feed tank pressure and proceeding towards the pumpo When the pump suction pressure is calculated, the pressure rise through the pump (PPUMP) is computed as the difference between the suction and discharge pressures. Since it is unlikely that this value of pressure head rise and given initial discharge will agree exactly with the given pump characteristic curve data, the amount (DPO) which the head

-16 rise vs. suction pressure should be shifted to produce the desired agreement is computedo That is, the pump speed is assumed to be changed slightly to obtain the desired agreement. Also the value of pressure rise, P0, corresponding to the initial discharge value is computed from the given head rise vso discharge datao This is discussed in more detail under the pump boundary condition, Constants pertaining to the engine operation are also computed, These are discussed in detail under the engine boundary condition. Information dealing with the print out of results is next established in the program. This includes column headings and the print out of some of the newly stored initial conditionso The balance of the program involves the calculations during the life of the transient conditiono It involves repetitive calculations which are performed every Ato This is accomplished by performing all calculations to the end of the program, incrementing time and repeating until sufficient results are obtained, The first repetitive calculations involve the interior points along each pipeline. Pressures and flow rates are computed at each time increment by use of the interpolation equations and EqsO (9) and (10) as described in Section IIIo At the same time increment as calculations are made at the interior points computations are also performed at the various boundary conditions. These are discussed in the next section.

V. BOUNDARY CONDITIONS The various boundary conditions are now treated individually in some detail. A. Series Pipe Connection At a junction of two pipes, as pipes 2 and 3, Figure 1, Q2 = QP3 and HP2 = Hp3 at any instant, In order to handle all similar series connections the numbers of all adjoining pipes are stored in pairs in the linear array IND2. A simple iteration statement then takes each series connection in sequence and obtains the solution at the new time. Figure 5 shows a typical series connection. Equation (14) is written for pipe X and Eq. (12) is written for pipe J. QPX = C3 - Cx PP Qj = C1 + CJ Pp Since QPX = QPJ9 the solution for pp is C3 - C1 p = CJ + CX With pp known, Qp can be found from either of the previous two equations. The pump inlet and discharge connections are treated as special series connectionso These are discussed under the pump boundary condition. Also if the excitation mechanism in the system is -17

PIPE X PIPE Figure 5. Series Connection _JAM cos (wt+a,) I SUCTION PIPE Figure 6. Connection at Feed Tank

-19 either a structural displacement of the feed tank or pump, an alteration is made to those particular series connections. Bo System Excitation at the Feed Tanks A structural motion of the tank with respect to the feed line is assumed if the variable TANK is set equal to one. The amplitude of the sinusoidal motion is AM1, the frequency is OMEGA, and the phase angle is ALPHA1, Fig. 6, The continuity condition at the junction becomes + QP = 1 2 where QP1 = p (PROP * Al - A2) (AM1) w sin (wt + al)~ The variable PROP can be used to describe the amount of the tank area that actually contributes additional flowo It would depend upon the shape of the pipe entrance and upon the number of suction lines from each tanko When this continuity condition is used with the earlier boundary condition equations, the pressure at the junction can be determined. C3 - C1 + QP1 p = CJ + CX Co Turbopump Representation In the analysis the pump is considered to be a pipe with a length equal to the mean particle flow path through the pumpo The head rise through the pump is considered a linear rise along its length and frictional losses are charged to the pump head rise. A pressure pulse wave speed is assumed through the pump which gives rise to a small time delayo Thus, the pump is treated as any other

-20 pipe in the system with interior point computations performed as described above, and it is connected to the suction and discharge lines using the same procedures as described above for series connections. The equation of motion takes a slightly different form when it includes the head rise produced by the pump. If the total pressure head developed by the pump is PPUMP, then an additional force PP — A Ax L can be considered to act in the direction of flow on the control volume shown in Figure 2, Then the equation of motion becomes -PPUMP dV p A Ax -PP A Ax - pgA Ax sin a + pA Ax d = 0 x L dt which reduces to the following form: Px PPUMP -- g- - g sin a + + W + Vt = 0 P p L When this equation is combined with the original continuity equation in the same manner as in section III, and placed in finite difference form, the following equations are obtained. P - QR + a (PP - PR) -A At (gp sin a + PPUP)= 0 (15) Xp - XR = a (tp - tR) (16) QP -S - A (PP - PS) - A At (gp sin a + PPIP)= o (17) a L xp - xS = -a (tp - tS) (18) These equations are used for the pump in place of Eqs. (5) to (8). In order to use this set of equations it is necessary that PPUMP be evaluated. The procedure used for this evaluation is next

-21 described. It is assumed that the steady state characteristic curves for the pump are valid during the transient conditiono It is further assumed that the shape of the curves remains the same although the pump speed may vary slightly during the unsteady flow. Tabular data of pressure rise vs. suction pressure (PNP vso PS) are given, starting at PS = Pi, at equal increments of PS, Figure 7. Since initial steady state conditions in the system are not likely to match the given pump curve exactly, the pump speed is assumed to be altered slightly so that the curve is shifted DPoo The initial suction pressure and pressure rise conditions in the system then fall on the newly positioned curve. This shifted curve is used throughout the analysis. At any suction pressure, the constants for a parabola (B1, B2, B3), passing through three adjacent points of the given curve, can be determined, The equation of the adjusted parabola is PRISE - B1 +DP + B2 (PS - ) + B S Pi)2 Tabular data are also given for the second pump characteristic curve, pressure-head rise vs. discharge at a constant speed, Figure 8 (PQ vs. Q, data given at equal increments of discharge, DQ)o Again for any discharge Q, the constants can be determined for a parabola (Zl, Z2, Z3) that passes through three adjacent given data points. Thus for the initial discharge Qo, the parabolic equation is Po = Z1 + z2 (o - Q) + Z3 (o - 1)2 A slight change in the assumed pump speed is likely to be necessary to shift the curve so the actual initial operating conditions fall on the

-22 PRISE INITIAL OPERATING CONDITION oo__-, - 0 i IGIVEN DATA DP PNP vs Ps h- - SUCTION PRESSURE Figure 7. Pump Characteristic Curve. Pressure Rise vs Suction Pressure PPUMP Po QI QO Q DISCHARGE Figure 8. Pump Characteristic Curve. vs Discharge Pressure Rise

-23 pump curve. The amount the curve is displaced is AP = PRISE - Po During the transient calculations, for a given suction pressure, B1, B2, and B3 are found and PRISE can be evaluatedO The difference PRISE - Po is the amount that the pressure rise-discharge curve must be shifted. With values of Z1, Z2, and Z3 determined with the aid of the most current value of discharge, the relationship between the current head rise and discharge is giveno PPUMP + Z2 (Qp - Qi) + Z3 (QP - 4 )+ PRISE- Po It can be noted that this analysis assumes that the pressure rise vso suction pressure curve is a constant discharge curve as well as a constant speed curve. For convenience and simplicity in the analysis, values of suction pressure and discharge from the previously computed time increment are used to evaluate PPUMP. D. System Excitation at the Pump A structural motion of the pump with respect to the suction line is assumed if the variable PUMP is set equal to oneo The assumed sinusoidal motion has an amplitude, frequency and phase angle of AM7, OMEGA, and ALPHA7, respectivelyo The continuity condition at the suction line junction with the pump is QP6 + QP6 = QP7 where QP6 = pA6 c(AM7) sin (wt + a7) When this continuity condition is used with the boundary equations

-24 presented under A, the pressure at the pump inlet can be determined. C3 - C1 + QP6 PP =CJ + CX E. Pulser If the variable PULSER is set equal to one in the input data, a sinusoidal pulser is assumed at location 5 of Fig, 1o A piston is assumed in a branch at the junction of pipes 4 and 6 as shown in Fig, 9. The motion of the piston is controlled in amplitude, frequency, and phase by the variable AM5, OMEGA, and ALPHA5, respectively. The continuity condition at the junction is given by %4 + P 5 = 6 where Qp4 is defined by Eq. (14) p4 = C3 - C4 PP QP6 is defined by Eq. (12) QP6 = C1 + C6 PP QP5 is defined by QP5 = pA5(AM5) u sin (wt + a5) The resulting pressure at the junction is given by the relation C3 - C1 + QP5 pp = C4 + C6

-25 Figure 9. Location of Pulser LO) LOX 9 A - -.. - < DOME -r-^^FUEL MANIFOLD A L10 LI, \ /QF\U X J9p10; ^^FUEL. PPL | PF 9 QL 4: ) -Id BURNER PLATE COMBUSTION ORIFICES CHAMBER PC IQt Schematic Diagram of Engine Combustion Chamber, Burner Plate Orifices, and Pump Discharge Lines Figure 10.

-26 Fo Pressure Volume Compensator If the variable PVC is set equal to one in the input data, the volumetric response of the PVC is assumed at location 5 of Figo o1 This volumetric response can be visualized as a piston in a branch as shown in Fig. 9. that moves in response to the physical displacement of the pipeo The variables that control the PVC motion are the same as those controlling the pump motion, namely AM7, OMEGA, and ALPHA7. The equations at the junction are identical to those for the pulser except QP5 is defined Qp5 - pA6(AM7) n sin (aot + a7) It is not possible to use the pulser and PVC in the program at one timeo Go Feed Tank Surface The pressure in the fuel and lox supply tanks is given by PTANK and is assumed constant during the transient analysis. With the pressure defined at a boundary the discharge can be computed directly from Eq. (12)o Ho Engine for Combined System Analysis The action of the engine beyond the burner plate orifices has been greatly simplifiedo It is assumed that the combined mass flow rate of fuel and lox, Qt, to the engine nozzle can be related to the combustion chamber pressure, PC, by the relation Qt = C30 Pc (19) The variable C30 is dependent upon the nozzle throat area and the characteristic velocity C*, which in turn is a function of the mixture

-27 ratio Mr (flow rate of lox over flow rate of fuel)( 2) Since small variations in the mixture ratio may produce significant changes in the characteristic velocity, the variable C30 is evaluated during the transient as the mixture ratio changes. Tabular data, characteristic velocity vs. mixture ratio, are given for the particular fuel and oxidizer used in the system. Values of C*, CSTAR, are given at equal increments, DMR, of mixture ratio beginning at MRlo For the initial combustion chamber pressure, total mass flow rate, and mixture ratio, the constant ATH is evaluated. ATH = C30 C* = C* Qt/Pc During transient conditions C30 is evaluated from the constant ATH and the C* that corresponds to the current mixture ratio. Figure 10 displays a schematic view of important component parts of the engine that are needed to describe its dynamic response in connection with the fluid pipelineo The notation used to describe the variables in the engine is showno The fluid in each of the feed lines passes into a volume, then through the burner plate orifices to the combustion chambero The volume V and effective compliance K' of the space at the orifice entrance are specified by the data, VOL and KP, respectively. The effective compliance, which represents the combined effect of the fluid compressibility and elasticity of the container, can be expressed K' = Ap/(AV/V)o For a small time interval the mass flow rate going into storage in the volume is Ql = p A At By combining these two expressions q10 = C0 Ap

-28 where C10 = p V/(K'At). As a first order approximation, Ap = pp - p, where pp is the unknown value of the pressure at the end of the time increment and p is the known value of the pressure at the beginning of the time increment. This expression is used in Eqs. (24) and (28) below. In addition to Eq. (19) and the continuity condition in the combustion chamber QL + QF = t (20) an equivalent set of equations must be written for the fuel and lox sides of the engine. For the fuel side these equations are as follows: The orifice relationship: F = (AOCD)F 2PF(PPF - P) (21) Continuity at the volume: ^F-^+ ^Fl0 = na(22) Characteristic equation from pipe 9: QF9 = C3F C9F PPF (23) Volume-compliance relationship: F10 = ClOF PPF - ClOF PF (4) The corresponding equations for the lox system are q = (AoD)L 4P (L(PPL - Pc) (25)

-29 QL + QL10 = QL9 (26) QL9 = C3L - 9L PPL (27) QL10 = C1L PPL - C1OL PL (28 At each new time increment Eqs. (19) to (28) must be solved for the ten unknowns, Qt, P c QF, QFl QF9, qF QL QL109 QL99 PPF, PPL~ Inasmuch as two of the equations involve quadratic terms, a direct solution is not possible. The equations are first combined, then an iteration procedure is used to obtain a solution. Equations (19) and (20) are combined to eliminate Qt 30 Pc = F + L (29 Equations (26) to (28) are combined to eliminate QL9 and QLiO L = L3 - L4 PPL (30) where L3 = C3L + CIOL PL and L4 =C9L + ClOL. Equations (22) to (24) are combined to eliminate QF9 and QF10 QF = F3 - F4 PPF (31) where 3 = 3F + COFF and F4 = C9F + C1OF Equations (21), (25), (29), (30) and (31) now involve five unknowns Pc, QF, QL, PPF9 and ppLo The method of solution involves the assumption of a small perturbation on each of the three pressures, a linearization of the orifice equations, then iterating and correcting the pressures until the perturbation is arbitrarily small. The final values are true solutions of the

-30 original nonlinear equations. Each pressure is expressed as the sum of a base pressure plus a variation from the base pressure. Pc Pc + P PPF = PF + PPF P =P +P, PL PL PL The orifice equation for the lox system, Eq. (25) can be written QLOR = ORL PL c +PLE P' L 4PPL - PC L PL - PcC + PPL' " PC' where ORL = (CDAo) l2pL o If the first two terms of a binominal expansion are used the orifice equation can be written QL' L 1 +Pp - c ) 2(FPL TcI or QL = L1 + L2 PpL' - L2 PC (25a) where L = ORL and (30) can be /PpL - Pc and L2 combined and PpL L1L4 + L2(L3 QL = L2 + = Li/(2 (PpL - Pc))o eliminated - L4 PPL) L2L4 L4 L2 + L4 Equations (25a).p I c (32) A similar manipulation of the equations for the fuel system produces a comparable equation.

-31F1F4 + F2(F3 - F4 PpF) F2F4 QF = - -— _ —- PC' (33) F2 + F4 F2 + F4 Equations (29), (32), and (33) are combined to give the solution for Pc PC' = L6/L7 (34) where L1L4 + L2(L3 - L4 PpL) F1F4 + F2(F3 - F4 PpF) L2 + L4 F2 + F4 and L2L4 F2F4 L7 = C30 + - L2 + L4 F2 + F4 Equations (25a) and (30) are combined to find PpLt; and PpF is found similarly. New base values for each of the pressures are now found by combining the previous base value and the newly computed variationo The process is repeated until the variations are very smallo In the program, four iterations are used to establish the new pressureso The four unknown discharges can be evaluated directly when the pressures are known, Io Engine when Systems are Isolated at the Combustion Chamber When the variable SINGLE is set equal to one in the input data, either the lox or fuel system can be analyzed alone without the interference of one upon the othero The fuel system analysis assumes that the lox flow rate to the combustion chamber is constant, QLBo When this assumption is made, Eqs. (19), (20), (21), and (31) can be combined and a direct solution can be obtained for QFQF = (- F + F + F2)/20

-32 where F1 = ORF + C30 -L -L \l )C^ and 2 F2 = ORF (F3 QLB SF4 C 30 J With QF known, a direct solution is available for the other variables, Qt c Pc, PpF, and QF9' A similar procedure is used when it is desired to isolate the lox system. The fuel flow rate, QFBY is assumed constant for this condition, J. System Excitation at the Burner Plate Orifice A structural motion of the engine with respect to the pump discharge lines is assumed if the variable ORIFIC is set equal to oneo The amplitude of the sinusoidal motion is AM9, the frequency is OMEGA, and the phase angle is ALPHA9. In the fuel system, the continuity condition between the end of the pump discharge line and the orifice plate becomes QF9 = QF + QFlO - QPF9 (22a) where QPF9 = p A9 o(AM9) sin (ct + ca9) Similarly, for the lox system, Eq. (26) becomes QL9 = + L10 - L9 where (26a) QPL9 = P A9 ((AM9) sin (t + a9).

VI. THE COMPUTER PROGRAM A copy of a listing of the program is shown in Figure 11o It is written in FORTRAN IV-G Level and was executed on the IBM 360/67 system at the Computing Center, The University of Michigan. The storage required by the program is approximately 42,000 byteso The execution time of the program is dependent upon the particular data being processed. The data that is necessary for execution of the program is listed. In all double subscripted data, the first subscript refers to the system (LOX = 1, FUEL = 2) and, except for the pump characteristics, the second subscript refers to the pipe number. D(l,l). o oD(l,9), D(2,1)...D(2,9) L(ll)oo oL(l,9), L(2,l),,.L(2,9) A(ll) o.. A(l,9), A(2,1) o o.A(2,9) F(l,)o.. oF(1,9), F(21)o o oF(2,9) NO(1 7l)o o NO(1, 9)) NO(2,l)o o oNO(2,9) SLOPE(l,l). oSLOPE(l19), SLOPE(2,1).o SLOPE(2,9) RHO(l), RH0(2), FLOW(l), FLOW(2), CDORIF(l), CDORIF(2), KP(1), KP(2), PTANK(l), PTANK(2), PCO, DQ, Ql(l), Q1(2), VOL(1), VOL(2), PQ(1,1)..PQ(1,20), PQ(2,l) o PQ(2,20), DP, PI(1), PI(2), PNP(l,1)...oPNP(1,20), PNP(2,l)o...PNP(2,20), IND2(1),. oIND2(12), J2, DMR, MR1, CSTAR(l)..oCSTAR(20), TANK, PULSER, PUMP, ORIFIC, OMEGA, AMl(l), AM1(2), AM5(1), AM5(2), AM7(1), AM7(2), AM9, ALPHAl(1), ALPHA1(2), -33

-34ALPHA5, ALPHA7, ALPHA9, SYSTM, SINGLE, N(2,7), G, TSTOP, UU, PVC, PROP The program prints out pressure and discharge information at the points indicated by P and Q in Fig. 12o The pressures are given in psi although all pressure computations are performed in psf in the program. The discharge unit is mass flow rate, slugs/seco Definitions of all program variables are given in Appendix A.

-35 FORTRAN IV G.LEVEL O, MOD 0 MAIN CSATURN 5 PROPELLANT SYSTEMS INCLUDING PUMPS AND ENGINE.CSUBSCRIPT K OR KK=1 REFERS TO THE LOX SYSTEM CSUBSCRIPT K OR KK=2 REFERS TC THE FUEL SYSTEM CSYSTM=1 LOX SYSTEM ONLY IS ANALYZED CSYSTM=2 FUEL SYSTEM ONLY IS -NALYZED CSYSTM=3 BOTH SYSTEMS ARE ANALYZED CSINGLE=1 SYSTEMS ARE ISOLATED AT COMBUSTION CHAMBER CPVC=1 THE.PVC IS INCLUDED IN THE ANALYSIS CALL PRESSURES ARE IN PSFA, ALL FLOWS ARE MASS FLOW RATES, SLUGS/SEC CPTANK(K) IS PRESSURE IN SUPPLY TANKS IN PSF CRHO (K) IS THE MASS OENSITY CTHE VARIABLES FOR THE AMPLITUDES AND PHASE ANGLES AT THE VARIOUS CEXCITATION POINTS ARE AS FCLLOWS (FREQUENCY OF OMEGA) CLOX TANK TANK=1, AMI(1), ALPHA(I) CFUEL TANK TANK=1, AMA12), ALPHA(2) CPUMP PUMP=1, AM7(1), AM7(2), ALPHA7 CORIFICE ORIFIC =1, AM9, ALPiA9 CLOX PULSER PULSER=i, AM5(1),ALPHA5 CFUEL PULSER PULSER=1, AM5(2), ALPHA5 0001 PARA- (DC1,DC2 DC3,DTH)=DC2+.5*DTH*(DC3-DC1+DTH* DC3+DC1-2.*DC2)) 0002 QRK R(Q1, CQ2,DTt )=DQ2-OTH*( CQ2-DQ1 ) 0003 PRKDP1,DP2,DTH)=DP2-CTH*( DP2-DP1 ) 0004 QSS(DQ2, C-3,DTH)=DQ2-DTH*(CQ2-DQ3) 0005 PSS(UP2, P3,TH)=DP2-DTF* (DP2-DP3) 0006 INTEGER UUUX,PULSER,ZORIFICPUMP,TANKSYSTM,SINGLEPVC 0007 REAL KP,L,L2,L3,L4,L6,L7,LMR1, NO 0008 DIMENSION 0(2,10),F(2,10),L(2,10),A(2,10),AR(2,10),THA(2,10)iC(2,1 20),N(210),FF210)F(2 QO(2,10),DPF(2,10),NO(2,10),SLOPEl2,10OACCEL( 32,10),P4(2,10), PQ(2,20),PNP(2,20,) OC00 DIMENSION RHO(2), FLOW(2), AM5(2), AM112), AM7(2), CDORIF(2), ORIF 2(2), VOL(2), KP(2), C10(2),Ql(2J)P0(2), PI(2),PPUMP(2), DPO(2), PT 3ANK(2), ALPHA1(2), INC2(20), CSTAR(20) 0010 DIMENSION Q(2,10,40), QP(2,10,40), P{2,10,40), PP(2,10,40) 0011 NAMELIST/INSE T/U,L,A,F,NC,SLOPE,RHO, FLOW,CDORIFKP,PTANK,PCODQQ1 2,VOL,PQ,DP,PI PNP IND2t J2 DMR,MR1 CSTAR, TANK,PULSERPUMPORIFICtOM 3EGAAM1,AM5,AM7,AM9, ALPHA, ALHAALPHA A7,ALPHA9,SYSTMSINGLE,N,G, T 3STOP,UUPVC, PROP 0012 10 READ(5, INSET) 0013 WRITE (6,INSET) 0014 U=0 0015 T=0. 0016 D T=L(2,7)/(N(2,7)*A12,7)) CINITIAL STEADY STATE CONDITICNS IN SYSTEM 0017 DO 60 K=l,2 00.18 PC=PCO 0019 DO 20 J=1,9 0020 N(K,J)=L( K,J)/(DT*AIK,J)) 0021 IF(J.NE.5) THAK, J)=N(K,J)*DT*A(K,J)/LtK,J) 0022 AR(K,J)=.7854*Di)K,.J)** 2*NO(K,J) 0023. C( K,J)=AR(KJ )/A( K, J) 0024 FF(K,J)=F(K,J)*DT/1 2.*DIKJ)*AR(K,J)*RHO(K)) 0025 QOiK,J)=FLOW(K) 0026 IF( J.EQ.5)O( K, 5)=0. 0027 IF( J.NE.5 )DPFIK,J )=F(K,J )*L(K,J )*QO(K,J)**2/(2. *D(K,J )*N(K,J)*AR 2K, J)**2*RHO(K))-RHO (K)*G*L(K,J)*SLOPE(K,J)/N(K, J) Figure 11. Listing of Saturn V Program

-36 0028 20 ACCEL(K, J)=G*RHO(K) *UT*SLOPE (K,J)*AR(K,J) 0029 DPf-(K,5)=0. 0030 OkF (K)=CDCGRIFU K)T*SRT(2.*RHO(K)) 0031 C 1O(K)=VUL(K) *RHO (K)/ ( T*KP(K) ) 0032 P(K,8,1 )=PC+(FLUI(K)/CRIF (K))**2+OPF(K,9)*N(K,9)+DPF(K,8)*N(K8) 0033 P ( K, 11 ).= PT ANK(K.) 0034 PP(K, 1, 1)=PTANK(K) 0035 WRITE i6,.311)DT,(N(K, I), =19), (AR(K,I),1 =1,9) 0036. 11 FORMAT (iHO,4H DT=,Fb.5,dH N(K,I)=,915,11H AREA(K,I)=,9F6.3) 0037 D050 J=1,S 0038 NN=N(K,J)+1 OC39 UC 30 I=1,NN 0040 P(,,J,I)=P(K,J,1)-(I-1)*UPF(K,J) 0041 30 t(K,J,I)=QO(K,J) 0042 IF(J-6) 50,40,50 0043 4C PPUMP(K)=P(K,d,l)-P(K,6,N(K,6)+1) 0044 DPF(K,7)=-PPUMP(K)/N(K,7) 0045 50 P(K,J+1,1)=P(K,J,N(K,J)+1) 0046 PP(K,9,N(K,9)+1 )=P(K,9,N(K,9)+1) CEVALUATION OF PUMP CHARACTEKISTIC CURVE CONSTANTS 0047 Z= ( (,6tNK,6 6) +1)-1 (K))/DQ 0048 P1=PQ(K,Z) 0049 P2=PQ(K,Z+1) 0050 P3=PQ(K,Z+2) 0051 L3=(P1+P3-2.*P2 )/(2.*L(*U') 0052 Z2=(P2-P1 )/DW-Z3*D*/( 2.*Z-1. ) 0053 Zl1=P2-DQ*Z**( Z2+Z3*DQ*Z) 0054 z=(( K,,N( K,6) +i)-PI (K) )/DP 0055 P1=PNP(KZ) 0056 P2=PNP(K,Z+I) 0057 P3=PNP(KZL+2) 0058 3= (P1+P3 —.*P2l/ (2. UP*CP) 0059 32= P2-P1 )/0P-b3*DP~ (2.*Z-1. ) 0060 bl=P2-DP*Z* ( B2+E3*UP* Z) 0061 PO(K)=Zi+Z2*((Q(K,6,N( K,)+.)-.1(K) )+Z3*(Q((K,6,N(K,6)+1)-Qi(K) )**2 0062 6C DPO(K)=PPUMP(K)-d1-32*(P(K,6,N(K,6)+1)-PI(K) )-B3*(P(K,6,N(K,6)+1 )2P I (K) )**2 OC63 (L=FLCW(1 ) 0064 QF=FLOW(2) 0065 (Q1 T=L+ QF 0066 LZ= (L/QF-Mki)/D'MR + 0067 TH=(QL/QF-tMR1-Z*DMR)/ECR 0068 CSTA=PAiAd((CSTA(<Z-1),CSIARI(Z),CSTAR(Z+1),TH) 0069 ATH=CSTA L^T/PC 0C70 C30=ATH/CSTA 0071 IF ( SI NGLE-1 )80,70,80 0072 70 QFi=QF 0073 QLb=QL 0C74 CF30=C30 0075 CL30=C30 00C76 0O F4=C10(2)+C(2,9) 0077 L4=C10(1)+C(1,9) OC78 JJ2'=2*J2 OC79 WRITE (6,630) OC80 630 FORMAT (IHO,8iH LIQUID OXYGEN SYSTEM 2 FUEL SYSTEM/111H TIME TANK OUT PULSER P 3UMP IN PUMP OUT ORIFICE TANK OUT PULSER PUMP IN PUMP OUT ORIF 4ICE COMBUSTION) Figure 11. (Concluded)

-37 0081 320 DO 321 K=l,2 0082 Do 321 J=1,9 0083 321 P4(K,J)=P(K,J,N(K,J)+1)/144. 0084 PC4=PC/144. 0085 NN=N41,4)+1 OC86 I'=N(2,4)+1 0087 WRITE (6,640)T,(P4(K,1),P4(K,4),P4(K,6,P4(K,7) P4(K,9),K=1,2)tPC4 2, ( 1,2, 1),Q( 1,4 NN),Q(1,7,1),Q(18,1 )QL,Q(2,2 1),Q(24 I),Q( 2 7'1 3), (2,8,1 ),QF,QT 0088 640 FORMAT (lHO,F6.4,3H P=,11F9.1/7X,3H Q=,llF9.2) 0089 330 T=T+DT 0090 U=U+i 0091 IF(T.GT.TSTOP) GO TO 240 0092 DO 120 K=i,2 0093 IF (SYSTM.EQ.2) K=2 CINTERIOR POINT COMPUTATIONS 0094 DO 100 J=1,9 0095 NN=N(K,J) 0096 DO 100 I=2,NN 0097 R=WR(Q( K, J I-1),{KJ, I ),THA(K,J)) 0098 Pk=PRR(P(KJ, -1) P(K,J, I ),THA(K,J)) 0099 QS=(SS(Q(KJ,I),Q(K,J, +1),THA(K,J)) 0100 PS=PSS(P(K,J,I),PK,J,I+1),THA(K,J)) 0101 Y1P(K,J, I )= 5*(QR+S++CK,( J)*PR-PS)-FF(K,J)*(QR*ABS(QR)+QS*ABS(QS} 2)+ACCEL(K,J) 0102 PP(K,J, I )=.5*PR+PS+(R-QS-FF(K,J)*(QR*ABS(QR)-QS*ABS(QS)))/C(K,J) 2) 0103 100 IF (J.EQ.') QP(K,7,I )=CP(K,7,I)+PPUMP(K)*AR(K,7)*DT/L( K,7) CSERIES PIPE CONNECTI CNS INCLUDING PUMP INLET AND DISCHARGE 0104 D0110 M=1,JJ2,2 0105 X=IND2(M) 0106 I=N(K,X)+I 0107 iR=4RR( Q(K,X, I-1),Q(K,X, I),THA(K,X)) 0108 C3= QR+C (K,X)*PRR(P(K, X, -1),P(K,X,I),THA(K,X))-FF(K,X)*QR*ABS(QR)+ 2ACCEL (K,X) 0109 IF (X.EQ.7)CC33+PPUMP(K)*AR(K,7)*OT/L(K,7) 0110 J=IND2(M+i) 0111 QS=QSS( Q(KtJ,1),Q(K,J,2),THA(KJ)) 0112 C1=QS-C(KJ)*PSS(P(K,J 1 ),P(K,J,2),THAK',J))-FF(KJ)*QS*ABS(QS)+AC 2CEL(K,J) 0113 IF(J.EG.7) C1=C1+PPUPFF(K)*AR(K,7)*UT/L(K,7) 0114 QP6=0. 0115 IF( PUMP.EQ.. AND. X. EC.6) P6=RHO(K)' AR(K,6)*AM7(K)*OMEGA*SIN(OMEGA* 2T+ALPHA7) 0116 QPI=O. 0117 IF( TANK.EQ.. 1AN. X.EQ.1 ) P1=RHO(K ) PROP*AR(K, )-AR(K,2)) 2*AM1(K)*OMEGA*SIN( OMEGA*T+ALPHAI(K ) 0118 PPIK,XI)=(C3-C1+QP6+CP1)/(C(K,J)+C(K,X)) 0119 PP(K,J,1)=PP(K,X, 1) 0120 QP(KX, I) =C3-C(K, X)*PP(K,X, I) 0121 110 P( K,, J 1) =Q P(K,X,I ) CPULSER OR PVC JUNCTICN 0122 I=N(K,4)+1 0123 QR=QRR( Q(K,4, I-1),Q(K,4,I ),THA(K,4)) 0124 C3='R+C(K,4)*PRR(P(K,4,I-1),P(K,4,I),THA(K,4) )-FF(K,4)*QR*ABSIQR)+ 2ACCEL (K,4) Figure 11. (Continued)

-38 0125 QS=QSS( QK,6,1),QI(K,6,2 ),THA(K,6)) 0126 C1=QS-C(K,6)*PSS(P(K,6t 1),P(K,6,2),THA(K,6))-FF(K,6)*QS*ABS(QS)+AC 2CEL<K,6) 0127 QPK,5, 1)=0. 0128 IF (PULSER.EQ. 1) QP K,5,1)=RHO(K)*AR(K,5)*AM5(K)*OMEGA*SlN(OMEGA*T+A 2LPHA5) 0129 IF( PULSER.EQ. O.AND.PVC.EG.1)QP4K,5, 1)=-RHO(K)*AR(K,6)*AM7(K)*OMEGA 2*SI N(OMEGA*T+ALPHA7) 0130 PP(K,4 I )=(C3-C1+QP(K,5, 1) )/(C(K,4)+C(K,6)) 0131 PP(K,6, 1)=PP( K4, I) 0132 PP K, 5,1)= PP K,6,1) 0133 (P (K,4 I) =C3-C( K,4)*PP(K,4 ) 0134 P (K,6 1) =C1+C( K,6)*PP(K,6, 1) CBOUNDARY CONDITION AT FUEL AND LOX TANKS 0135 QS=QSS(Q(K,1,1),Q{K,1,2)tTHA(K,1)) 0136 C1=QS-C(K,1)*PSS(P(K 1, 1),P(K t,2),THA(K-l))-FFK, 1)*QS*QS+ACCEL(K 2,1) 0137 QP(K,,l1)=Cl+C(K,1)*PP(K,1,l) 0138 120 IF (SYSTM.EQ.1) K=2 CBOUNDARY CONDI TIN A T ENGINE 0139 I=N(1,9)+1 0140 QR=QRR(Q(1,9,I-1),Q(11, 9I),THA(1,9)) 0141 C3-=QR+C(1,9)*PRR(P( 1t,I-1), P(1,9,I ),THA( 1,9)-FF( lt9)*QR*QR+ACCEL 2(1,9) 0142 J=N(2,9)+1 0143 QR=QRR(Q(2,9,J-1),Q(2,9,J),THA(2,9)) 0144 C2=QR+C(2,9)*PRR(P(2,9,J-1,P(2,9,J),THA(2,9))-FF('2,9)*QR*QR+ACCEL 2(2,9) 0145 L3=C3+C 10 (1 )*P( 1,9, I) 0146 F3=C2+C10(2) *P( 2,9, J) 0147 IF(ORIF IC-1)140,130,140 0148 130 QPL9=RHO( 1)*AR( 1,9) *AM9*OMEGA*SIN(OMEGA*T+ALPHA9) 0149 LPF.9=RHO(2)*AR 2,9) *AM9*OMEGA*SIN(OMEGA*T+ALPHA9) 0150 GO TO 150 0151 140 QPL9=0. 0152 QPF9=0. 0153 150 IF(SINGLE-1)190,160,190 0154 160b IF(SYSTM.EQ.2) GO TO 170 0155 LI=OR IF ( )**2*( 1./L4+1./CL30) 0156 L2=CRIF(1)**2*( (L3+QPL9)/L4-QFB/CL30) 0157 QL=.5*(-L1+SQRT(Ll*Li+4.*L2)) 0158 QT=,L+QFB 0159 PC=QT/CL30 0160 PP( 1 9,1I ) =(-QiL+L3 ) /L4 0161 QP(1,9,I)=C3-C( 1,9)*PP(1,9,) 0162 Z=(QL/QFB-MR1.)/MR+1 0163 TH=I QL/ QFB-MR1-.Z*MR) / MR 0164 CL30=ATH/(PARAB (CSTAR(- STARZ-1STAR(Z),CSTAR Z+1 ),TH ) 0165 170 IF(SYSTM.EQ.1) GO TO 210 0166 FI=ORIF(2)**2*(./F4+1./CF30) 0167 F2=CRIF(2)**2*{ (F3+Q(PL9)/F4-QLB/CF30) 0168 QF=.5*(-F1+SQRT(Fi*FI44.*F2)) 0169 QT=QF+QLB 0170 PC=QT/CF30 0171 PP( 2,9, J)=(-QF+F3 )/F4 0172 QP (2,9, J)=C2-C 2,9)*PP(2,9,J) 0173 Z=( QLB/jQF-MR1 )/DMR+1 Figure 11. (Continued)

-39 TH= (QLB/QF-MR 1-'Z*DMR) / D M CF30=ATH/(PARAB(CSTAR(Z-1),CSTAR(Z) CSTAR(Z+ l),TH)) GO TO 210 190 DO 200 M=1.,4 L =ORIF (1 )*SQRT (PP( 19, I) -PC) L2=L1/( 2. * PP(1,9, I )-PC ) F =ORIF(2)*SQRT(PP(2,g9,J)-PC) F2=F1/( 2.*( PP (2,9, J.)-PC) ) L6=-C30*PC+( L IL4+L2* ( L3-L4*PP(1,9,I) +QPL9)) / (L2+L4+ (F 1*F4+F2*(F3 2 -F4*PP 2,9,J)+QPF9))/(F2+F4) L 7=C30+L2*L4/( L2+L4 )+F2*F4/ F2+F4) PCP=L6/L7 PLP=(L3 -L4*PP 1,9, I )-L1+L2*PCP+QPL9 )/(L2+L4) PFP={F3-F4*PP(2,9, J)-F1+F2*PCP+QPF9)/(F2+F4) PP(1,9,I)=PP(1,9, I+PLP PP(2,9 J)=PP(2,9,JJ+PFP.200 PC=PC+PCP QP( 19, I)=C3-C(1,9)*PP(1,9,I) QP (2,9,J)=C2-C(2,9)*PP(2,9,J) QL=GRIF (1 )*SQRT'(PP( 1,9., 1 )-PC) QF=CRIF 2 )*SQRT (PP(2, 9 J )-PC) QT=QF+QL. Z=(QL/QF-MR1) /DMR+1 TH= ( QL/QF,-MR1-Z*DMR)/CMR C30=A'TH/(PARAB(CSTAR1Z-1),CSTAR(Z),CSTAR(Z+1),TH)) CSUBSTITUT1ON STATEMENTS FOR PRESSURE AND DISCHARGE 210 DO 230 K=1,2 IF (SYS'TM.EQ.2) K=2 DO 220 J=1,9 NN=( K, J) +1 DO 220 I=1,NN Q(K,JI )=QP(K,JI) 220 P(KJI )=FPP(KJ ) CPUMP CHARACTERIST-IC CONSTANTS AND PUMP PRESSURE RISE Z=(Q(K,6sN(K,6)+1)-Q1(K))/DQ Pl=PQ(KZ) P2=PQ(K,Z+1) P3=PQ(KZ+2) Z3=(P1+P3-2.*P2)/ (2.DQ*DQ) Z2=IP2-P )/DQ-Z3*DQ*( 2.*Z-1. ) ZL=P2-DQ*Z*(Z2+Z3*DQ*Z) =(P(K,6,N(K,6) +1)-PI (K))/DP P1=PNP K,.Z) P2=PNP(K, Z+) P3=PNP(K,Z+2) B3=( P l+P3-2.*P2 )/ (2. *DP*DP) B2=(P2-P1 )/DP-833*DP* 2.*Z-1.) B1=P2-DP*Z* ( 2+ B3*DP* Z) PPUMP(K)=Z1+Z2*(Q(K,6,N(K,6)+1)-QlK) )+Z3*4Q(K,6,N(K, 6)+1l-QI(K) ) 2**2+Bl+DPO(K) +B2* (P(K,6,N(K,6)+I)-PI K) )+B3*(P( K, 6,N(K,6)+11-Pi(K.3) )**2-PO(K) 230 IF (SYSTM.EQ.1) K=2 IF(U/UU*UU-U) 330-,320,330 240 STOP END Figure 11. (Continued)

-4o 1 LOX TANK + r 2 FUE L. L TAN K LOX SYSTEM K = 1 SYSTM =I FUEL SYSTEM K=2 SYSTM = 2 COMBINED SYSTEMS SYSTM = 3 I PULSER OR PVC LOX PUMP FUEL PUMP LOX DOME FUEL MANIFOLD BURNER PLATE COMBUSTION CHAMI Figure 12. Sketch of System Showing Location of Program Output Quantities P, Q.

VII. SAMPLE PROBLEM In order to demonstrate the use of the program, a sample problem is included which analyzes a transient condition in Titan II missile. The program analyzes the transients that develop in the oxidizer system as a result of a forced motion of the pumpo The motion is assumed to be sinusoidal at a frequency of 6208 rad/sec and half amplitude equal to 0,007 ft. Figure 13 gives a copy of the listing of the input data and Figure 14 shows a portion of the output from the calculations. The effect of the addition of a PVC in the system can be seen in Fig. 15, where a comparison is made between two runs. Both runs are subjected to the same excitation at the pump, one includes the PVC, the other does noto -41

EXECUTION BEGINS &INSET D= 10.000000 0.41700000 0.41700000 3.5999994 2.0000000 483.00000 1750.0000 2980.0000 0.11000000E-01, 0.0 NO= 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 0.0 f KP= 100000000. 5. COOOO0, 189000.00 166000.00 1.0000000 1.0000000 1.0000000 176000.00 176000.00 176000.00 176000.00 1.CO0000 IND2= 1, 10.000000 0.41700000 0.42999995 0.42999995 1.5000000 985.00000 3640.0000 4230.0000 0.11000000E-01, 0.0, 1.0000000 1.0000000 1.0000000 1.0000000 0.0 RHO= 2.8199997 100000000. VOL= 0.0 224000.00 188000.00 1.0000000 1.0000000 1.0000000 190800.00 200000.00 203000.00 204000.00 1.0000000 2, 0.57499999 0.500000000 0.57499999, 0.50000000 ~ 0.57499999, 0.50000000 0.57499999, 0.50000000 t 0.57499999, 0.50000000, 0.41700000, 0.42999995, 0.0, 0.0,L= 8.0000000 y 8.0000000, 23.000000, 0.39999998 1.5000000, 1.6699991, 0.0, 0.0, 0.50000000, 0.50000000, 2.5000000, 6.0000000, 3.0000000, 14.000000, 0.0, 0.0,A= 1750.0000, 138.00000, 1100.0000 ~ 138.00000, o50.00000, 138.00000 650.00000, 138.00000, 2980.00.00 3000000 2980.0000 4230.0000, 0.0 0.0,F=.79999976E- 0.79999976E- 02, 0.1100000 O E-01, 0.11000000E-01. 0.72999954E-01, 0.12099999, 0.11000000E-01, 0.11000000E-01, 0.72999954E-01, 0.12099999 1.0000000, 0.35999998E-01, 1.0000000, 1.4799995, 0.0 t 0.0, 1.000000, 1.0000Q00 1.0000000, 1.0000000 I 1.0000000, 1.0000000 1.0000000 1.0000000 1.000 0 0000000, 1.0000000 1.000000 10000000 1.0000000, 1. 0000000 SLOPE= 1.0000000.0000000 1.0000000 1.0000000 1.0000000, 1.0000000, 0.0 0.0 1.0000000 v 0.0, 0.0, 0.0, 0.0, 0.0 0.0, 1.7599993,FLOW= 17.000000, 8.1699991,CDORIF= 0.36699999E-01, 0.20999998E-01,,PTANK= 10500.000, 2460.0000,PCO= 115000.00,DQ= 1.0000000,Q1= 12.000000, 0.0,PQ= 197000.00, 231000.00, 193000.00, 229000.00, 184000.00, 218000.00, 179000.00, 210000.00, 173000.00, 200000.00 159000.00 174000.00, 150000.00, 170000.00, 141000.00, 165000.00 1.000000 00 1.0000000, 1.0000000, 1.0000000, 1.0000000 1.0000000 1.0000000 t 1.0000000, 1.0000000, 1.0000000, 1.0000000 1.0000000 1.0000000 1.0000000,DP= 250.00000,PI= 10500.000, 2000.0000 tPNP= 176000.00 194500.00, 176000.00, 196800.00, 176000.00 198800.00 176000.00, 201000.00, 176000.00, 201900.00, 176000.00, 202500.00 176000.00, 203400.00, 176000.00, 203700.00, 176000.00, 203900.00 176000.00, 204100.00, 1.0000000, 1.0000000, 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000, 1.000000 2, 3 3, 4, 6, 7, 7, 8, I rI 8, 9, O, O, 6, DMR= 0.99999964E-01,MR1= 1.7999992 5650.0000, 5610.0000, 5565.0000 0 0.0, 0.0. TANK= 0,PULSER= 1 PUMP= 0,,CSTAR= 5775.0000 0.0, 0.0 0.0, 0.0 1,ORIFIC= Ot O0, 5756.000OC y 0.0, 0.0 O,OMEGA= 62.7S 0, ) 5730.0000 t 0.0 I 0.0 99988,AM1= 0.0,ALPHAl= 0.0 1,SINGLE= O0 0, n. n 1i, 5695.0000, 0.0 t 0.0 y 0.0 t 0.0 ltN= O, 0, O, 0,G= 32.199997 2= 9 9 I AM5= 0.0, 0.0 ALPHA5= 0.0,ALI O, O, 0 0O TSTOP= 0.63999951E-01,UU= &END PHA7= 0.0 ( 0,AM7= 0.29999998E-02, 0.0 ALPHA9= 0.0 I) O0 O,,. O O0,AM9= 0.0,SYSTM= 0, O0 0, 0, 1 8 P VC= OPROP= 0.19999999 Figure 13. Input Data for Example Problem

DT= 0.00050 N(Kl)= 33 26 6 4 0 1 1 1 2 AREA(KI)=78.540 0.260 0.260 0.260 0.137 0.260 0.260 0.137 0.137 OT= 0.00050 N(K,I)= 16 5 6 24 0 7 1 2 6 AREA(KI)=78.540 (YGEN SYSTEM FUEL SYSTEM 0.196 0.196 0.196 0.137 0.196 0.196 0.145 0.145 LIQU[D OX TIME 0.0 P= Q= 0.0040 P= Q= 0.0080 P= 0= 0.0120 P= Q= 0.0160 P= Q= C.02CO P= Q= 0.0240 P= Q= 0.0280 P= Q= 0.0320 P= Q= 0.0360 P= Q= 0.0400 P= Q= 0.0440 P= Q= 0.0480 P= Q= 0.0520 P= Q= 0.0560 P= Q= 0.0600 P= Q= 0.0640 P= _~ TANK OUT 78.0 17.00 78.0 17.00 78.0 17.00 78.0 1-7.00 78.0 17.00 78.0 16.99 78.0 16.97 78.0 16.94 78.0 16.91 78.0 16.88 78.0 16.86 78.0 16.85 78.0 16.84 78.0 16.84 78.0 16.84 78.0 16.86 78.0 16.90 92.0 17.00 92.3 16.98 92.9 16.96 93.5 16.94 94.1 16.92 94.7 16.91 95.1 16.91 95.4 16.91 95.5 16.92 95.4 16.93 94.8 16.94 93.6 16.96 92.2 16.97 90.6 16.99 88.9 17.00 87.3 17.02 85.9 17.03 92.0 17.00 92.4 16.97 92.9 16.95 93.5 16.93 94.1 16.92 94.7 16.91 95.1 16.91 95.4 16.91 95.5 16.92 95.3 16.93 94.8 16.95 93.6 16.97 92.1 16.99 90.5 17.01 88.8 17.02 87.2 17.03 85.9 17.05 PULSER PUMP IN PUMP GUT ORIFICE TANK OUT PULSER 1314.4 17.00 1314.2 17.00 1314.0 17.00 1314.0 17.00 1314.1 17.00 1314.2 17.00 1314.4 17.00 1314.6 17.00 1314.8 17.00 1315.0 17.01 1315.1 17.01 1314.7 17.00 1314.3 17.00 1313.7 17.00 1313.1 16.99 1312.5 16.98 1312.0 16.98 1062.8 17.00 1062.7 17.00 1062.6 17.00 1062.5 17.00 1062.5 16.99 1062.6 17.00 1062.7 17.00 1062.8 17.00 1C62.9 17.00 1063.1 17.00 1063.2 17.01 1063.1 17.01 1C62.9 17.00 1062.6 17.00 1062.2 16.99 1061.8 16.98 1061.5 16.98 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 20.2 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 19.8 8.17 PUMP IN PUMP OUT 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 19.6 1401.2 8.17 8.17 ORIFICE COMBUSTION 1097.2 798.6 8.17 25.17 1097.2 798.6 8.17 25.17 1097.2 798.5 8.17 25.17 1097.2 798.5 8.17 25.17 1097.2 798.5 8.17 25.16 1097.2 798.5 8.17 25.17 1097.2 798.6 8.17 25.17 1097.2 798.6 8.17 25.17 1097.2 798.7 8.17 25.17 1097.2 798.7 8.17 25.17 1097.2 798.8 8.17 25.18 1097.2 798.8 8.17 25.18 1097.2 798.7 8.17 25.17 1097.2 798.5 8.17 25.17 1097.2 798.3 8.17 25.16 1097.2 798.2 8.17.25.15 1097.2 798.0 8.17 25.15 I I Figure 14. Output from Analysis with Oscillating Pump

110 WITHOUT PVC C) 100 L ii / \WITH PVC r ~90.5- -O.15 -" cc 90 Q0 (n 80 aQ. 70 0.05 0.1.15.20.25.30 TIME, SEC Figure 15. Pressure at Pump Inlet of Oxidizer System of Titan II Due to Pump Oscillation, w = 62.8 rad/sec.

VIIIo SUMMARY This report includes the digital computer program for the analysis of the propellant systems of Saturn Vo Much of the background material necessary to understand the computations performed by the computer during execution of the program is includedo Data for Saturn V have not been used in any trial runs of the program as complete data on the geometric properties of the missile are not availableo The program was cleared using Titan II datao Computed results have been checked against experimental results obtained on Titan II and reported in the literature (3) and favorable comparisons have been madeo It is our intention that the program can be used directly for any one of the three stages of Saturn V without further modificationo However, since the exact configurations are not available, minor changes may be necessary to accommodate the particular datao

APPENDIX A DEFINITION OF PROGRAM VARIABLES ACCEL(k,j) ALPHAl(k) ALPHA5 ALPHA7 ALPHA9 AML(k) AM5(k) AM7(k) AM9 AR(k, j) A(k,j) ATH BlB2,B3 ClO(k) C1 C2 C30 C3 CDORIF(k) Constant associated with vehicle acceleration in pipe j of system k Phase angle associated with the feed tank motion in system k Phase angle associated with the pulser motion Phase angle associated with the pump motion Phase angle associated with the engine burner plate motion Amplitude associated with the feed tank motion in system k Amplitude associated with the pulser motion in system k Amplitude associated with the pump motion in system k Amplitude associated with the engine burner plate motion Cross-sectional area of pipe j in system k Pressure pulse wave speed in pipe j in system k Constant associated with the engine Variables associated with the pump head rise-suction pressure curve Constant in system k describing the capacitance of fluid and chamber at the engine orifice plate inlet Variable associated with the C- characteristic Variable associated with the C+ characteristic Variable relating combustion chamber pressure and total propellant mass flow rate Variable associated with the C+ characteristic Constant used to describe the discharge coefficient and area of opening of the engine orifice in system k -46

-47 CF30,CL30 C(k,j) CSTAR(i) CSTA DMR DPO(k) DPF(k,j) DP DQ D(kj) DT F1,F2,F3,F4 FF(kj) FLOW(k) F(k,j) G IND2(i) I J2 Variables used with combustion chamber boundary condition Constant for system k, pipe j, AR(k,j)/A(k,j) Variable used for storage of tabular data of the engine characteristic velocity vso the mixture ratio Variable identifying the characteristic velocity at a particular mixture ratio Interval of mixture ratio values at which the characteristic velocities are stored beginning at MR1 Adjustment of the tabulated pump data, head rise vs. suction pressure, to match initial flow conditions Steady-state frictional head loss per section in pipe j in system k Interval of suction pressure values at which the pump head rise values are stored beginning at PI(k) Interval of discharge values at which the pump head rise values are stored, beginning at Ql(k) Pipe diameter, pipe j, system k Time interval at which computations are performed Variables used in connection with the fuel supply to the engine Constant associated with the friction losses in pipe j of system k Given steady-state mass flow rate in system k Darcy-Weisbach friction factor in system k, pipe j Acceleration in the missile axial direction at the time in flight when the transient is being investigated Set of integer constants that identify the pipes, in pairs, by number, that are connected in series Integer variable, generally refers to the pipe section number under consideration Constant identifying the number of series connections in each system

JJ2 J KP (k) K L1,L2,L3,L4 L6,L7 L(k,j) MR1 M NO(k,j) N(k,j) NN OMEGA ORIFIC ORIF(k) PO(k) P1,P2,P3 P4(k,j) Constant = 2(J2) Integer variable, generally refers to pipe under consideration Constant defining the effective bulk modulus of volume at the engine orifice plate entrance Integer variable, generally refers to the system, K = 1 refers to lox system, K = 2 refers to fuel system Variables used in connection with the oxidizer supply to the engine Variables used with the boundary condition at the engine Pipe length, pipe j in system k Lowest value of mixture ratio at which characteristic velocities are stored Integer variable, used in series connections and engine boundary condition Integer constant identifying the number of identical pipes connected in parallel that can be treated as a single pipe j in system k Number of reaches into which pipe j of system k is divided Integer variable used in place of N(kj)+l Angular frequency of excitation Integer used to include or exclude excitation at the engine orifice plate Constant used with the orifice in system k Pump head rise associated with initial steady state discharge, as determined from the given head rise vso discharge data in system k Variables associated with the pump head rise-discharge curves Pressure at outflow end of pipe j, system k, in psi; used for print out

-49 PCO PCP PC PFP PI(k) PLP PNP(k,i) PP(k,j,i) PPUMP(k) PQ(k, i) PROP PR, PS P(k,j, i) PTANK(k) PULSER PUMP PVC QO(kj) Ql(k) Initial steady state pressure in the combustion chamber Variation in combustion chamber pressure at new time increment Pressure in combustion chamber Variation in pressure upstream of orifice plate in fuel system at new time increment Lowest value of suction pressure at which pump head rise values are given in system k Variation in pressure upstream of orifice plate in oxidizer system at new time increment Values of pump head rise stored for equal increments of suction pressure, DP, in system k Newly computed values of pressure at section i, pipe j, system k Pressure head produced by the pump in system k Values of pump head rise stored for equal increments of discharge, DQ, in system k Constant used to define the proportion of the feed tank area to be used when the system excitation is occurring at the tank Interpolated value of pressure head in the pipeline Pressure head at section i, pipe j, system k Pressure level maintained in feed tank of system k Integer used to include or exclude excitation of the system using the pulser Integer used to include or exclude excitation of the system using the pump motion Integer used to include or exclude the pressure volume compensator in the system Initial steady state flow in pipe j of system k Lowest value of discharge at which pump head rise values are given in system k

QFB QF QLB QL QP1 QP6 QPF9 QPL9 QP(k,j,i) QRg QS Q(kgj i) QT RHO(k) SINGLE SLOPE(k,j) SYSTM TANK THA(k,j) TH T Constant flow rate of fuel to combustion chamber when lox system alone is being analyzed Flow rate of fuel to combustion chamber Constant flow rate of lox to combustion chamber when fuel system alone is being analyzed Flow rate of lox to combustion chamber Variable used when the system excitation is considered to be at the supply tank Variable used when the system excitation is considered to be at the pump Variable used in the fuel system when the system excitation is considered to be at the engine orifice Variable used in the lox system when the system excitation is considered to be at the engine orifice Newly computed values of flow rate at section i, pipe j, system k Interpolated value of flow rate Flow rate at section i, pipe j, system k Total flow rate to the combustion chamber Fluid mass density in system k Integer constant used to isolate the systems at the combustion chamber Constant used to indicate whether a pipe is in the axial direction of the missile or normal to the axial direction Integer constant used to specify what systems are being studied Integer used to include or exclude excitation of the system using the motion of the supply tank Constant used in pipe j, system k; pertains to the characteristics grid-mesh ratio Variable used with parabolic interpolation Time

-51 TSTOP U UU VOL(k) X Z1, Z2, Z3 Z Constant identifying the time at which the transient calculation should stop Integer counter, incremented by one when time is incremented by DT Integer constant controlling the print out Constant defining the volume at the engine orifice plate entrance Integer used to identify the inflow pipe number at series connections Variables associated with pump head rise-discharge curve Integer used with the pump characteristics Function References PARAB Function used for parabolic interpolation PRR Function used for linear interpolation of pressure at R PSS Function used for linear interpolation of pressure at S QRR Function used for linear interpolation of discharge at R QSS Function used for linear interpolation of discharge at S

REFERENCES 1. V. L. Streeter, and E. Bo Wylie, Hydraulic Transients, McGraw-Hill, Inc., 1967. 2. R. H. Fashbaugh, "Liquid Propellant Missile Longitudinal Oscillation," The University of Michigan, NASA Report, Nov., 1966. 3o V. L. Streeter, and R. H. Fashbaugh, "Resonance in Liquid Rocket Engine Systems," Trans. ASME, Jour of Basic Engineering, Vol 87, (Dec., 1965), pp. 1011-1017. -52

UNIVERSITY OF MICHIGAN 39015 035249559