011C FILE' .

NUSC Technical Report 6150

25 May 193

Transient Compressible Flow in a

Piping Network: A Solution Method
and Computer Simulation
1
00 Paul J. Lefobvre
V-" Richard F. Hubbell

(0 Launcher Systems Department

03

EDTIC
JUN 0 3 1988

H
a6

Naval Underwater Systems Center

Newport, Rhode Island / New London, Connecticut

Approved for public release: distribution unlimited.

88 6 . 4l
PREFAO;E
L

This report documents work performed under

Independent Research and Independent Exploratory

Development (IR/IED) Project No. A43115,

"Compressible Pipe Flow Simulation." Principal
investigator was Paul J. Lefebvre (Code 3711).
Associate investigator was Richard F. Hubbell
(Code 3711). Technical reviewer was William G.
Fenneil (Code 3634).

REVIEWED AND APPROVED: 25 May 1983 L

W. A. Von Winkle

Associate Technical Director for Technology 0

.
|a
,4
1. IMPORT MNJE
REPORT DOCUMENTATION PAGE
2. GOVT ACCESSION NO.
(7
3.
READ I s-MUCTIMOS
BEFORE COMPLETN FOR."

R CMIENTS CATALOG NUMIER

TR 6150
4. TL1 kad SubauLer S. TYPEOF REPORTA MUDD0COVERED

TRANSIENT COMPRESSIBLE FLOW IN A PIPING NETWORK:

A SOLUTION METHOD AND COMPUTER SIMULATION ,. ERo .om NUMBER

7. AIIINOR61 I. CONTRACT 001 GRANT RUl,

Paul J. Lefebvre and Richard F. Hubbell

S.PWORMS OfGMATMO "ME WODANARES TAWN

t&. PROGRAM SLAMNT. PROJECT".
&WOW UNT NUMIR
AMUA
Naval Underwater Systems Center A43115
Newport, Rhode Island 028I
11.
CONTROLLN OFCE M AN ADDRESS 12. REPorT DATE

Naval Underwater Systems Center 25 May 1983

0
NMBEROFPAGES
Newport, Rhode Island 028I 1

AGENCY NAME I ADDRESS lit dufom fem CawoUi.

14. kIONITORING OfWc I . SECURITy CLAS. lf thi pe")

UNCLASSIFIED
1S OE9CLASS1FICATIONI
DOWNGRADING
SCHEDULE

I. OISTIEITION S'MhiffT W d Rampj

Approved for public release; distribution unlimited.

1?.. I TRUUTIN
STATENUT Wt he .6wowe.,. BLbehA3 if dffw'.w jr'. R-Po 0

IL SUPP EENTARY NOTES

I S. IKY WORDS #Conmaa - mean* o.deif sm,.ary *ad .. nfy by Wch nansbe,

Fluid Dynamics,
Compressible Pipe Flow. Computer Program

Z&.
ABSTRACT fCanow.. on mune side if na..wsev an ideann.by bdael.aw..har

-" A general method of calculating the transient state of compressible

fluid flow through complex piping networks has been developed and a
computer program implementing this metho" has been written. This report
presentz the theory and logic upon which the solution method and program
are based. It serves also as a basic introduction to the program, which
is written for persons having minimal knowledge of compressible flow
theory.

IDLO, 1473
4.%
TABLE OF CONTENTS

Page

LIST OF ILLUSTRATIONS .......... .. .................... ii

LIST OF TABLES .......... ... ........................ ii
INTRODUCTION ........... ... ......................... 1
DISCUSSION ............ ... .......................... 3
Theory ............. .. .......................... 3
Assumptions ........ ........................ 4
Strategy ......... ........................... 4
Runge-Kutta Integration Scheme ...... ... .............. 7
Convergence Technique .... ..... .. .................. 9
Throttling Device; Forward Pass ... ............. ... 13
Throttling Device; Backward Pass ... ............. ... 15
PROGRAM STRUCTURE................................. .. 17
Main Program ........ ....................... ... 17
Main Subroutines ....... ..................... ... 18
Real Function Subroutines .... ................ .19
DEFINING THE NETWORK ....... .................... ... 20
Organization of Input Data .... ................ ... 20
Throttling Device and Tank Calculations . ........ .22
PROGRAM EXECUTION ......... ..................... ... 24
OUTPUT ........... ............................ .. 24
SUMMARY AND CONCLUSIONS ...... .................... ... 25
REFERENCES. . .......... ... ........................... 26
APPENDIX A -- FLOWCHARTS ...... ................... ... A-i
APPENDIX B -- SAMPLE INPUT DATA FILE ... ............. ... B-i
APPENDIX C -- SAMPLE OUTPUT DATA FILE ... ............. ... C-i
LIST OF ILLUSTRATIONS p

Figure Page

1 Compressible Flow Network ........ ................ 2

2 Runge-Kutta Integration Scheme ....... ............. 8
3 Ratios of Mass Flow Rates ..... ............... .... 12
4 Throttling Device in Pipe Section ...............
.... 13

LIST OF TABLES

N
Table Page

la Small y-Value Convergence .... ................ .. i

lb Large y-Value Convergence ..... ............... .. .

.-

DTIC TABE

Gr2_p/ I
-- - ---.
coP . } Ji .t S"pct i,,._

ii ~lst pec U..

Io' i"
TRANSIENT COMPRESSIBLE FLOW IN A PIPING NETWORK: -
A SOLUTION METHOD AND COMPUTER SIMULATION

INTRODUCTION

This report describes work conducted under the NUSC Independent

Research and Independent Exploratory Development project titled
"Transient Compressible Pipe Flow Simulation." The main objectives of
this project were to develop numerical techniques for calculating 0
transient compressible flow of a fluid through complex piping networks
and to create a user-oriented computer simulation sufficiently general to
be applicable to a wide variety of networks. The end result of this
effort is a computer program called COMP, which has been implemented at
the Naval Underwater Systems Center on a Digital Equipment Corporation
VAX 11/780 computer.

Several simulations for adiabatic flow of compressible gas in

relatively complex piping networks have been developed over the past 20
years (see references 1 through 5). Each, however, was written for a
specific application and consequently has serious drawbacks when applied
to other networks. Specifically, the simulations are not user-oriented
(making it difficult to define and set up the network), are not applicable
to choked flow across throttling devices such as nozzles and valves, and
are not readily applicable to transient flow where a combination of
conditions (e.g., constant pressure, constant volume, or variable volume)
may exist at supply and receiving tanks.

The simulation described in this report overcomes these drawbacks.

In addition, it provides a new technique for rapid convergence of mass
flow rate distributions at junctions common to two or more pipe lines,
can account for any ideil gas, and employs numerical techniques for all
calculations requiring iteration or integration. It can be expanded to
accommodate practically any piping network.

%i

M
T " ' ' 2.. " " I ,
•. '. ,
-' P, J' ' .~.2r W' , ... . -. .'.'_ .' .= " " "-"
The simulation assumes adiabatic flow in each line. It accounts for
irreversibilities encountered in throttling devices and for frictional
effects in the pipe and fittings by treating them as equivalent pipe
lengths. Figure I shows the network used in the simulation. Circled
numbers in the figure indicate tanks, uncircled numbers indicate pipe
sections, and circled letters indicate junctions. Pipe sections connect
two junctions or a tank and a junction. Each pipe section can consist of
up to 10 pipe lengths arranged in series. Throttling devices can be
substituted for pipe lengths, but cannot be the first or last lengths in
a section.

01
11 ,
PPol

DIRECTION 1
OF
FLOWP01

10
79

a S

i 13
Figure 1. Compressible Flow Network

Figur .k Compressible Fow Network

I-

DISCUSSION

THEORY

The general approach used here in solving the transient compressible 0

flow problem in a piping network is similar to that used in reference 1.
Namely, adiabatic frictional flow is assumed to exist in each line of the
network and the system's mass flow rate is determined by balancing mass
flow rate and stagnation pressure at a single point (called the common 0
point), which is subject to the total system mass flow rate. For the N
initial time increment in this transient flow problem, calculations first
proceed in a forward direction from the most remote supply tank to the
common point, and then in a backward direction from the most remote 0

receiving tank to the common point. Stagnation pressure values at the

common point as calculated by these forward and backward computations or
"passes" are then compared. If they are not approximately equal, the
mass flow rate calculations undergo repeated iterations until they
balance. Once convergence of stagnation pressures at the common point
occurs, the simulation employs a Runge-Kutta integration scheme to set up
tank conditions for the beginning of the next time increment, whereupon .

the calculation process for stagnation pressures at the common point is

repeated.

The equations used to calculate compressible adiabatic flow are taken

from reference 6 and are not repeated here. The convergence scheme,
however, was developed along with the simulation itself. The scheme is
programmed as one of the simulation's main subroutines (subroutine RATIO)
and is discussed along with other subroutines later.

3 S:~
ASSUMPTIONS

Assumptions made in conducting this analysis are as follows:

1. Flow is one-dimensional, quasi-steady, and adiabatic.

2. Ratio of specific heats is constant.
3. Ideal fluid compressibility factor, Z, is constant and equal to
unity.
4. Friction factor, f, for any length of pipe is equal to the

average of the inlet and outlet steady-state friction factors. N

5. Minor losses due to elbows, reducers, etc., are introduced as

equivalent pipe lengths.
6. No reverse flow exists in any line.
7. Isentropic stagnation pressures are balanced at each junction of
two or more pipe lines.
8. Isothermal thermodynamic processes exist in each tank.

STRATEGY
NV

The junctions between pipe sections are assumed to be isentropic,

since any frictional effects due to the junctions may be compensated for
by adding equivalent pipe lengths to the network. The isentropic
condition at junctions forces the conclusion that there must be a single
common stagnation pressure for all pipe lengths common to a junction.
This assumption is critical to the analysis because it leads to the
condition that the stagnation pressure is the property to be iterated
upon at each junction in order to determine mass flow rate distribution.

4.

.',
X0FX - .W A

One important exception to this condition is evident when throttling

or choking occurs in a pipe terminating at a junction. The stagnation
pressure calculated at the end of this pipe may be greater than the

junction stagnation pressure, with the pressure difference being lost in

the choked pipe as the flow enters the junction.

The overall approach to calculations in the forward direction is to

assume a mass flow rate in the most remote and always present supply tank.
This is tank 1 and mass flow rate ;l in figure 1. Throughout the
transient there must be some absolute flow from this tank since it
controls the overall calculation process. Calculations proceed to the

first junction (junction A in figure 1), where a stagnation pressure at

the end of pipe section 1 is calculated. The program next assumes a
mass flow rate from tank 2 into pipe section 2 (m 2 ), and calculates a
stagnation pressure at junction A from section 2. If this stagnation
pressure is within 1 psi of that calculated from section 1, calculations
proceed to section 4; if it is not, adjustment is made to m 2 , via the
convergence scheme (subroutine RATIO), while 1l is held at the initial
value. As soon as agreement is reached (usually 2 to 5 iterations on
m 2 ), an arithmetic average of stagnation pressures at junction A is
made and calculations continue.

Mass flow rate i 4 is then taken to be the sum of mi and m2 .

The starting pressure for section 4 is taken as the junction stagnation
pressure obtained from sections 1 and 2. A stagnation pressure at

junction B from section 4 is then calculated. The iterative procedure

used to calculate the tank 2 flow is next employed again for tank 3 to
obtain a value for m3 and a stagnation pressure at junction B. This
completes the forward pass, since junction B is the common point where
the forward and backward calculations meet.

5i
The backwrard pass begins at the receiving tanks and proceeds toward
the common point, junction B. It is assumed that whatever mass flowed
into the system from the supply tanks will flow out of the system into
the receiving tanks during any particular time step. To ensure that this
condition is met, the program takes a percentage of the mass flowing into
junction B from the forward calculations, and uses that as the mass flow

in section 9. The remainder goes to section 10. The flow in section 9

is then divided between sections 11 and 12 (initially 50 percent in

each). Calculations proceed to determine the stagnation pressure at

junction F that would be required to produce that estimate of mll .

Stagnation pressure at junction F from section 12 is calculated next and

compared to the pressure from section I1. If the two do not agree, a
mass flow rate adjustment is again made via subroutine RATIO.

Once agreement is reached, m 9 is used to back-calculate the

stagnation pressure required at junction E for a flow rate of m9 . The

program proceeds identically in sections 13 and 14 to junction E. If the
stagnation pressure at junction E from section 9 does not match the one
from section 10 within limits set into the program (presently 1 psi), then
subroutine RATIO adjusts m 9 and mlO until agreement is reached. The
correct ; 9 and mlO values are rapidly obtained, resulting in a
calculated stagnation pressure at junction E from the backward
calculation. Note that if m 9 and ma1 are adjusted, then ; 1 1
through m14 must also be adjusted and all backward pass calculations
repeated from the tanks.

Calculations proceed to the beginning of section 8 where stagnation

pressure at junction D is obtained. The flow then splits into the
parzillel branches of sections 6 and 7. Under the initial assumption that
the flow is equally split between sections 6 and 7, the pressure at
junction C based on each section's contribution is then determined.

6 S
If agreement is achieved within the set limits of 1 psi, the parallel

branch calculations are complete and the simulation calculates the

conditions for section 5. If there is nG agreement, subroutine RATIO .4

adjusts the percentage of flow going to section 6 (which then gives a new
flow rate in section 7 since m7 = m5 - m 6 ) until convergence is
achieved. Stagnation pressure at junction B is then calculated. This
completes the backward pass, since junction B is the common point.

The overall convergence scheme of subroutine RATIO is once again

accessed. It now compares the stagnation pressure at junction B from the
forward pass with stagnation pressure from the backward pass. If they
match (within 1 psi) the network is solved for that time step. If they
do not match, an adjustment to ml is made via subroutine RATIO and the
forward and backward pass calculations are again repeated with the new
value for m1 . Iterations on mi continue until stagnation pressure
values at common point B as calculated from the forward and backward -.

passes agree. At that time, all distributed mass flow rates and junction
pressure values are correct, values for that particular time increment -

are output to the output file, and calculations for subsequent time
increments begir!. Initial conditions used for the subsequent time
increment are determined by a Runge-Kutta integration. P

RUNGE-KUTTA INTEGRATION SCHEME

The accuracy of the program in predicting the transient behavior of

the network depends, in part, upon the technique selected to update the
conditions within the supply and receiving tanks. The experience of the
authors with different numerical integration schemes led to the selection
of a fourth-order, Runge-Kutta technique, which is ideally suited to this
particular application.
.. ,

-S:
The technique is straightforward. The initial condition of each tank
with respect to volume, temperature, pressure, and type of gas being used

is known. The initial mass of gas contained within each tank is first
calculated by the perfect gas law, the forward and backward passes are

run, and the overall convergence criteria are applied at the common
point. Once convergence at the common point is achieved, the mass flow

rate from (or into) each tank, as well as that through all the pipes, and
across each Junction is known. With the initial mass within a tank and
mass flow rate known, it is possible to calculate a new mass within the
tank by means of a Runge-Kutta integration. Figure 2 illustrates this
concept. Once the new mass is found, the perfect gas law is once again
used to obtain new pressures under the assumption of isothermal tank
processes. Having new pressures at the beginning of the next time
increment allows the calculation process to begin again.

KNOWN ORIGINAL
/ MASS

KNOWN SLOPE, m,

NEW PREDICTED
MASS

MASS
(Ibm) I

TIME (s)
ONE TIME
STEP

Figure 2. Runge-Kutta Zntegration Scheme

8%
In the event that the tank has a variable volume, the user may employ

a differential equation to describe the volume rate of change of the tank

with respect to time. The differential equation can be non-linear, but
must be ordinary. In running the simulation, the user inserts the

equations into the FORTRAN main program in a section clearly labeled with

comment statements. An example would be the following equation:

dV2/dt = 2 x Time x Time + 32.0

which would be inserted into the program as

F(2) = 2 * TIME * TIME + 32.0

in the section marked for tank 2. It would also be necessary to°

initialize the volume of tank 2 in this example at time = 0. The FORTRAN

statement is Y(2) = 32.0, which would indicate an initial volume of 32,'
cubic feet.

CONVERGENCE TECHNIQUE

The assumptions that the flow is adiabatic and frictional provide no .

explicit equation for pressure drop in terms of a mass flow rate through
" ,X ' "za,"
. "', "'- J -X:
c"a-.-"". "" ,: ' -;-:
' ",'; '::":;:'; :'' <:-'.':<'-, '-': ;',• : ::< :,:., pi"
a pipe, as would be the case for isothermal pipe flow. Instead, the ".
polminvolves a set of very non-linear algebraic equations that must
be solved simultaneously. The problem experienced when the wrong massl
flow rate is initially selected is that the stagnation pressures do not
match at the junctions. For example, consider section 2, with 2
flowing to junction A. At junction A a value for stagnation pressure

from sectioni has been obtained andt is now necessary to arrive at the
correct value ofm 2 In order to have the junction A stagnation pressure
from section 2 match the pressure from section 1. In short 2 must be

9|
made to converge to the correct value. Since the equations being used to
calculate junction pressures are non-linear, a technique had to be
devised to cause m 2 to converge rapidly. Indeed, one of the main
contributions of this study is the method that was developed to cause the
mass flow rates to converge. The method is extremely fast and general.
For example, with an initial estimate of m2 of 0.1 lbm/s when the
actual converged value should be 100.0 lbm/s, the method used here will
yield the correct value in only six iterations.

The method is as follows: an initial value of ;2 is estimated from

which a stagnation pressure at junction A is calculated. If the estimate
does not result in convergence it must be adjusted and the procedure
repeated. If the stagnation pressures from sections 1 and 2 are far
apart at junction A, then a relatively large adjustment must be made, and
the direction in which the adjustment must be made is known. If the
value calculated for section 2 was too low at the junction, then m 2
must be reduced so there is less pressure drop through section 2, thus
increasing the section 2 stagnation pressure value at the junction. If

the calculation for section 2 was too high at the junction,

then ; 2
must be increased in order to lower the section 2 pressure value. It is
not desirable to simply add or subtract a fixed amount to rn2 because it
is not possible to know ahead of time how large m 2 should be.

For efficient convergence, a method was developed that is sensitive

not only to the size (and sign) of the pressure difference at junction

A, but to the size of ;m2 as well. The method uses the exponential %
function:

Ynew = (Yold) exp (Constant * x)

where y 2 and x = junction A stagnation pressure from section 2

minus the stagnation pressure from section 1.

10
%!
vrvv4z'%i ...

The choice of constant will depend on how quickly convergence is

desired and how sensitive x is to small changes in y. The constant value I'

in'the simulation as it currently exists is good for most situations, but

may have to be changed for very small pipe diameters. Tables la and lb
illustrate convergence for two very different size y-values when the
constant equals 0.00005.

Table la. Small y-Value Convergence

I
Yold Ynew

1000.0
1 0.1 0.1051271
I
+ 50.0 0.1051271 0.1053903
- 10.0 0.1053903 0.1053376
+ 1.5 0.1053376 0.1053455
- 0.2 0.1053455 Convergence 0.1053444

Table lb. Large y-Value Convergence

x Yold Ynew

1005.0 85.0 89.380385

+ 500.0 89.380385 91.643061
+ 300.0 91.643061 93.028068
+ 60.0 93.028068 93.307571
- 5.0 93.307571 93.284247 0
+ 0.2 93.284247 Convergence 93.285180

This convergence method is also used to adjust the ratio of mass flow
rates at a junction having two pipe sections leaving or entering from it.

11
At junction C, for example, RATIO(6) can be defined as ;6/mtotal. I
Then, by adjusting RATIO(6) using two equations (see figure 3) it is
possible to split the flow by the appropriate amount:
S.

a = RATIO(6)ne w

* RATIO(6)old exp(Clx) (for C1 > 0 and x > 0)

* 1 - exp(-C2x ) RATIO(6)old * exp(-C 2 x) (for C2 > 0 and x < 0).

'p

RATIO 1.0
(6

RATIO1 (1 - a) "-P(C2x) / RATIO a exp (-C1 x) ,O

4- NEG x x POS

Figure 3. Ratios of Mass Flow Rates F

The following summary shows where in the piping network the convergence
scheme of subroutine RATIO is applied, as well as the the order of use:

I. At junction A to adjust m

2. At junction B to adjust m 3
3. At junction F to adjust•mi3
3. At Junction G to adjust m.4

5. At junction E to adjust r9

6. At junction C to adjust m6
7. At junction B for overall convergence, to adjust n, I

Uses I and 2, above, are for the forward pass; uses 3 through 6 are for
the backward pass.

12
"
''.5
TV MTV N*

THROTTLING DEVICE; FORWARD PASS

Throttling devices control, monitor, or limit flow rate in a pipe.

Examples are valves, orifices, and nozzles. The effectiveness and/or
discharge coefficient of such devices may be a constant or may be a
function of time, Reynolds number, or some other parameter. Throttling
devices are represented schematically as shown in figure 4.

UPSTREAM DOWNSTREAM

Figure 4. Throttling Device in Pipe Section

The equations used in the program to define the mass flow rate
through a throttling device are taken from reference 7. These equations,
which were derived for orifices with pipe taps, are as follows:

Y 1 -0.333 + 1.145 (2 + 0.78

2 * 12B3)]
3 P 1 7P
" P2 () "

1/
2
m CDYA(2gP I(P 1 - P) (unchoked condition) (2)

mrit : CDAP -TY 1 2 (choked condition) (3)

L3 L

13

-ft. ft
P2 crit P1 crit \ y +1(4)

Po2 - Po0 + P2 - 1 (5)

where:

A = throttling device reference area (ft2 )

CD = discharge coefficient I
2
g = acceleration of gravity (ft/s )

= mass flow rate (lbm/s)

2
P1 = upstream pressure (lb/ft )
2
2 = downstream pressure (lb/ft )
R = gas constant [(ft-lbf)/(lbm-OR)] (53.3 for air)

T =temperature (OR)
8 = diameter ratio (throttling device/upstream pipe)
p= weight density (lb/ft 3 )
Y = specific heat ratio

and the subscript "o" specifies stagnation properties and "crit" specifies
properties at critical (throttled or choked) flow. I

The conditions upstream of the throttling device are known prior to the
point in the computer program at which the throttling device calculation
takes place. These known conditions include values for m, M1 , Pol, and I
To,. The calculation proceeds as follows:

1. The critical mass flow rate (mcrit) is calculated via equation (3). N
S

2. If m is greater than mcrit' then the estimate of flow rate from

tank 1 is too high and is reduced by mcrit/ . Calculations
begin again from tank 1.

~. ~ ~ ~ ~w.~ga *... '' ~ ~ w~ .:.-

14
3. Ifin in

.
If -- erit , then the flow is choked and will pass m. For this
crit' oh and (5).
case downstream properties are calculated via equations (4)

4. If is less than crit the flow is not choked and will pass M.
in

Downstream properties are calculated next via equations (1), (2),

and (3).

THROTTLING DEVICE; BACKWARD PASS

When a throttling device is encountered in the backward pass

calculations, downstream conditions (point 2 in figure 4) will be known
(i.e., ;, M 2 , Po2, and T. 2 ). The main difference between the forward
and backward pass throttling device calculations is that the mass flow rate
is not iterated in the backward pass when ; is greater than mc. Instead,
the upstream pressure is adjusted to accommodate the required mass flow
rate. The specific procedure is as follows:

1. Downstream static pressure, P2 , is calculated.

2. Critical upstream static pressure, Pl crit, is calculated via

equation (4).

3. Unchoked upstream static pressure, P1 unchoked, is calculated via

equations (1) and (2). Note that this is the upstream pressure
required to produce a mass flow rate across the device equal to M.

4. If P1 unchoked is less than or equal to P1 crit, then

I unchoked is the upstream static pressure, PI, and the I

upstream stagnation pressure is calculated via equation (5). '-

1'

15
5. If P 1 unchoked is greater than P1 crit, then the flow across I
the device in choked and upstream stagnation pressure, Polo is
.9.
calculated via equation (4). In this case, when total convergence
for the entire network is satisfied all the properties downstream
of the throttling device will have to be recalculated in a forward
pass fashion. The reason for this is that the mass flow rate is so
high in this line that there is a loss of stagnation pressure from
the end of the line into the tank.

This concludes the calculations across a throttling device. Although

most of the calculations are straightforward, calculating the unchoked

upstream static pressure, P1 unchoked' of step 3 can be difficult. Here,

this was overcome by means of a Newton-Raphson technique (reference 8).
This technique was implemented by obtaining a single algebraic equation at

point 1 for p1 as a function of m, effective area, discharge coefficient,.'S

specific heat ratio, and the static pressure at point 2. This equation,
which was solved by the Newton-Raphson technique for pl, is as follows:

PA
m2) 2 /
1* 0 -IP1 1 0i
0' E2g(p2 RT°- G R) i

where

C* is the product of discharge coefficient (CD) and the effective

area of the throttling device,

and% GD - m -

2
A2 gYRT

16
Next, an equation was derived to calculate the static pressure at

point 1 in terms of i

2p

2
P, RT Y ~
" 1 o 2gY - P1

which is the value for unchoked static pressure at point 1.

PROGRAM STRUCTURE

MAIN PROGRAM

The main program (program MAIN) initializes all constants and reads

all input data from an input file. Next, the initial estimate for the
controlled flo,, rate from tank 1 is specified via the terminal. Once all
the data are fed into the computer, program MAIN calls various
subroutines to calculate the pressures and flows at the pipe junctions
throughout the system.

Program MAIN calls subroutine CTANK to calculate pressures in all

pipes leading away from the input tanks. Subroutine REVERSE is called
upon to calculate the pressures and flow rates from the receiving tanks

back to junction D. Once the flow is known at junction D, subroutine

REVERSE calculates all necessary pressures and flow rates in the parallel

branch. REVERSE then evaluates the pressure and flow rate in the last
pipe (section 5) of the backward pass to the common point, junction B.

Program MAIN then checks the difference in stagnation pressure at the ft.--

common point. If the difference is large, MAIN will change the control .
flow rate via subroutine RATIO, and will repeat the entire calculation

' 1C
.~

17
ZLC
- ~~ I - -7777771s

I.

.2

procedure until agreement is reached and state conditions at all

junctions of the flow network are satisfied. MAIN will print the
results, increment the time step, and repeat the procedure.

Appendix A presents flowcharts of the main program and major

subroutines. A brief description of the subroutines follows.

MAIN SUBROUTINES

Subroutine PRES calculates the end pressure for a pipe length in the

forward pass under the assumption that specific heat ratio, gas constant,
mass flow rate, stagnation temperature, stagnation pressure, pipe
dimension, and pipe length are known.

Subroutine TANKPR calculates tank stagnation pressure at the end of

each time step using a Runge-Kutta integration.

Subroutine BKVALVE performs throttling device and pipe flow

calculations for the backward pass. A.

Subroutine FVALVE performs throttling device calculations for the

forward pass.

Subroutine CTANK calculates end pressures and flow rates of pipes

leaving the tanks in the forward pass network and in pipe section 14. %

Subroutine RATIO, the convergence subroutine, adjusts the flow rate,

using an exponential function, until the correct value for the flow rate

is found.

'A
~-~
'JW~' .~jA ~~ -~ 5 . 'A A \ A ~ Is~~
5
A 'A~ A *\* * A, ~ -. **%

18
Subroutine REVERSE controls the backward direction calculations,
which calculate end pressure and flow rate for all junctions between the
common point and the receiving tanks. It also calculates all pipe length
pressures and flow rates in the backward pass.

REAL FUNCTION SUBROUTINES

Real function FRIC determines the Darcy friction factor at a point in I

the pipe where the Mach number, stagnation temperature, stagnation
pressure, specific heat ratio, gas constant, mass flow rate, and pipe
diameter are known.

Real function DYNVIS uses Sutherland's equation (reference 9) to

calculate the dynamic viscosity.

Real function MALPHA calculates the Mach number at the end of the
pipe length (outlet condition), using fL*/D ratios:

fL*/D = (fL*/D)inlet - (fL*/D)outlet

where L* is the length of duct required to develop a flow from the Mach
number at the position under consideration to the sonic point, f is the
average friction factor, and D is the pipe diameter.
I

Real function BETA solves the equation that relates Mach number to
friction given a known Mach number, M, and specific heat ratio, Y:
I,

L'D -M (Y l)2
)M
2
7M2 2Y 2 + (Y - l)M

19

V- U.. A. -44
Real function FF is a digital equivalent of the Moody chart for
friction factor as a function of pipe roughness and Reynolds number.
(Units are in feet with roughness = 0.00015 for commercial steel pipe.)

Real function MARAT calculates the Mach number at the pipe inlet.

Real function TEMP calculates the temperature at any point given the
Mach number and specific heat ratio.

Real function CVEL calculates the speed of sound for a specified

temperature.

Real function VEL calculates the velocity of the fluid at any point
in the pipe given the Mach number.

Real function RHO calculates the density of the fluid at any point in
the pipe given the Mach number and velocity.

Real function REND calculates the Reynolds number at any point given
the Mach number, temperature, density, and velocity.

DEFINING THE NETWORK

ORGANIZATION OF INPUT DATA

Input data for the compressible pipe flow computer program is set by
the user via a computer file. Figure 1 is the model for the input data,
and networks up to the complexity shown in figure 1 can be simulated.
Tank 1 must always be present, but any of the ozher tanks can be deleted .

as long as at least one receiving tank is available. Pipe sections can

also be deleted as long as continuity in the flow path is maintained.
Each pipe section can be defined by up to 10 individual pipe lengths,

20

% -
.w

4-

each with its corresponding diameter. Throttling devices are substituted

for pipe lengths, but they cannot be located as the first or last lengths
of a pipe section. This is of no serious consequence, since if a
throttling device located at the end of a pipe section in a real network
must be simulated, a relatively short pipe length of large diameter
(which would have a negligible effect on the calculaticns) can be input
in the simulation after the throttling device. The organization of the *

input file (see appendix B for a sample input data file) is as follows:

Line 1 - Identification line for the problem as defined in columns 1

to 110 of the data file.

Line 2 - Stagnation temperature (OR) in the network. The FORTRAN S

format is F10.4.
k'S

Line 3 - Initial stagnation pressure (lb/ft 2 ) in each tank. The

format for this line is 7F10.4 and is in the following tank order: 1, 2,
3, 11, 12, 13, 14. If a tank is not present, stagnation pressure of 0.0

is input. v

Line 4 - Initial volume (ft 3 ) of each tank. The format for this S
line is 7F10.4. Tank order is the same as for line 3. If a tank is not
present, volume 0.0 is input.

Line 5 - Flags specifying which tanks are in the system. Foidlat is 0

714. Tank order is the same as for line 3. An integer value of 1 is .

used if a tank is in the system; 2 if a tank is not in the system.

Line 6 - Flags to specify the conditions in each tank. The format

for this line is 714 with the tank order the same as for line 3. An
integer value of 0 is used if a tank is not present; 1 for constant
pressure; 2 for constant volume; 3 for variable volume defined by
user-supplied differential equation. I

21

%%
Lines 7-8 - The number of different pipe lengths in each pipe
section. Maximum is 10. Format is 714. Information is given in the
following order:

Line 7: Pipe sections no. 1, 2, 3, 4, 5, 6, 7

Line 8: Pipe sections no. 8, 9, 10, 11, 12, 13, 14

If a pipe section is not present in the network an integer value of 0 is

input. Each throttling device counts as a single pipe length.

Lines 9- - Pipe length and diameter data (ft) for each pipe listed in
lines 7 and 8. For each pipe section present in the network as shown in
the order given in lines 7 and 8, two lines of information are required.
The first specifies the length of each pipe in that section in the order
of direction of flow; the format is IOFlO.4. The second line contains
the diameter data in the same order and format as the length data. Even
though the format of each line is lOFlO.4, only data for the number of
lengths in that section are required. If a section is not present no
lines of data are input for it. If a throttling device is to be
simulated, a value of -99.0 is input for length and diameter. (Area and
discharge coefficient data for the throttling device are inserted into
the throttling device subroutines as discussed in the next section.)

Last line - Contains the values of specific heat ratio and gas
constant (ft-lbf)/(lbm-°R) in that order and in format 2F10.4.

THROTTLING DEVICE AND TANK CONDITIONS

Throttling device and tank condition data are still required to

completely define the network. Since the respective subroutines are

based on equations that are peculiar to a particular type of throttling

22
device or tank condition, the user must edit the program to change
equations and constants for discharge coefficient, area, mass flow rate,
expansion factor, etc. The subroutine source listings contain comment

statements to help the user.

Throttling device calculations are performed by subroutine FVALVE for

pipe sections 1 through 4 and subroutine BKVALVE for pipe sections 5
through 14. In these subroutines there is a separate calculation for

each pipe section so that different throttling device characteristics can 0

be defined. The mass flow rate and expansion factor equations that are

presently used in the subroutines are shown in the "Theory" section of

this report. Should these not apply to a particular throttling device,%
then other equations must be substituted. In any case, discharge
coefficient and related area must be inserted in the form of a constant .

or a function of time (time is transmitted to the subroutine through its N

argument).

Subroutine CTANK calculates tank conditions at the end of each time

step using the Runge-Kutta integration scheme. When a variable volume
tank condition is specified (via an integer value of 3 on line 6), the %
main program must calculate the rate of change of volume with respect to
time via a user-supplied differential equation prior to calling the CTANK

subroutine. As explained earlier, the equation to define the slope for

the variable volume must be inserted into the main program. Similar to
the throttling device subroutines, these calculations use separate
calculation sections in the main program for each tank. The CTANK
subroutine assumes an isothermal thermodynamic process between time
steps, but can be changed by the user to include other processes as

required. I

This concludes all requirements for defining the network.

23
d. o.I 'r _ e

% % % %
~
*5'.%S.*S~* ~ _r.A %I%'S
:Y. - - -
I'.'%~ ~ S%~I~~%S a
PROGRAM EXECUTION

The compressible flow computer program is presently configured to run

on a Digital Equipment Corporation VAX 11/780 computer with the user
providing setup data via interactive terminal. Prior to any
calculations, the program asks the user for the input and output data
file names, the time increment interval in number of time steps at which
results are to be written to the output file, the integration time step -

size (s), the maximum run time (s), and the initial estimate of mass flow
rate (lbm/s) from tank 1. This is all the information required to
execute the program. Note that if the user sets integration time step
size equal to the total run time, then the program provides a
steady-state network analysis.

Program termination occurs when one of three criteria are met: when
the maximum specified run time has been reached, when the pressure in the
supply tank(s) is less than or equal to that in the receiving tank(s), or
when flow from tank 1 has stopped. As mentioned in the "Theory" section,
tank 1 is the main supply tank and controls the calculations; hence it
must always produce a flow in the downstream direction. ..

OUTPUT

Appendix C contains the output file which resulted from executing the 0
program using the sample input file of appendix B. Note that all input
data contained in the input data file are initially written to the output
file, permitting verification of input data. Results of program
calculations at each time step for which results were to be written are
listed next. At each output time step, the following data are written to
the output file:

I. Time averaged mass flow rate (lbm/s) through-each section

2L4
" "
g' 4.."
." .4"
" ." ,. .---. r ., . . ' ' ." ,,',_o
" - ' ' ' ," . ' " . ,'- " ° ." - ' '' " ' ' ' '', '' '4"' ''',m- SA-
I .- '
2. Time averaged stagnation pressure at the end of each pipe section
3. Time averaged Mach numbers at the end of each pipe section
4. TPnk pressures at the end of the time step

5. Tank volumes at the end of the time step.

SUMMARY AND CONCLUSIONS

A method of calculating the transient flow of compressible fluids

through complex piping netwrks has been developed, and a computer
program implementing this method has been written. This report has
documented the theory and logic upon which the solution method and

computer program are based. It also serves as an introduction to the

piogram, which was designed to be used by persons having minimum

knowledge of compressible flow theory.

I.r

The program is shown to be sufficiently general to handle a wide

variety of complex piping systems. Supply and receiving tanks can each
be defined as having either constant volume, constant pressure, or
variable volume; flow can be either steady or unsteady; flow-limiting
conditions resulting from wall friction, changes in pipe diameter, and
throttling devices can be simulated; and the number, length, diameter,
and configuration of pipe sections can be easily varied. Numerical
techniques such as Runge-Kutta integration and Newton-Raphson iteration
have been incorporated in the program wherever possible. The program,
besides being versatile and adaptable to a wide variety of piping
networks, provides rapid solution convergence.

A%L "a')

-1 J

/5*...F-,--. '-- ' ~ ~ %/N 5

.\% ~~d~
25
REFERENCES

1. T. G. Keith, Digital Computer Analysis of Submarine Emergency Debal-

lasting System, Piping Section, DTNSRDC Report 2569, David W. Taylor
Naval Ship Research and Development Center, Bethesda, MD, May 1968.

2. M. E. Wilton and J. P. Murtaugh, Compressible Flow Network Analysis,

COSMIC Report No. LEW-11859, Computer Software Management and
Information Center, Athens, GA, 1971.

3. S. Krause, On One-Dimensional Frictionless Compressible Gas Flow

with Heat Addition and Area Change Under Particular Consideration of
the Transonic Region, Report No. DLR-FB-67-02 (Royal Aircraft
Establishment Library Translation No. 1237, Farnborough, UK), 1967.

4. J. F. T. MacLaren, et al., Solving the Equations for Unsteady Gas

Flow by Two Methods, Simulation Councils, Inc., 1976.

5. P. T. Daniel, "The Analysis of Compressible and Incompressible Fluid

Networks," Transactions of the Institution of Chemical Engineers,
Vol 44, 1966.

6. A. H. Shapiro, The Dynamics and Thermodynamics of Compressible Fluid

Flow, Vol 1, Ronald Press Company, New York, 1953.

7. J. P. Holman, Experimental Methods for Engineers, McGraw-Hill Book

Company, New York, 1971.

8. B. Carnahan, H. A. Luther, J. 0. Wilkes, Applied Numerical Methods,

John Wiley and Sons, New York, 1969.

9. F. M. White, Viscous Fluid Flow, McGraw-Hill Book Company, New York,

1974.

2I

- . ~~ ~ •'~*j -~ m r*~---
.~(* a~ -' ' v v -: |
APPENDIX A

FT-O WCHA TS
-

%fAIN PROGRAM

0IMENSION ALL VARIABLES

READ INP.T DATA FILE

A
DIETERMINE INITIAL MASS IN EC

ESTIMATE FLOW11)

CALL SUBROUTINE CTANK TO FIND

AND FLOW
PRESSURES T RATES
! FROM
SUPPLYTNKS; CALCULA E FCESSURE AT
I
COMMON PINT JUNIN 3

''
FLOW RATES IN
PR C3URIS AN D)

RECEIVING TANKS ANO ALL PIPES SACK TO

JUNCTION 8
P
SCOMPARE PRESSURES AT JUNCTION a 1

_ _ FROM
_ SUBIROUTINES
_ _ _ _
_ CTAINK
_ _.,.'I
CA LL SUBIROUTINEI p T __ OBITAINED
_

ADJUST FLOWAI) AND REVERSE. (TMIS IS THE OVERALL

55
OVALOPSS IS COMPL-ET& CORRECT
FLOW(I) IS KNOWN: CORRECT MASS FLOW
RATE OISTRIBUTION AND JUNCTION
STAGNATION PRESSURES FOUNO

WRIT OUTRESUTSAS A FIRST ESTIMATE USE

F - UT REULOSDo
WRT' FLOW)I

A MAXIMUM TIME SIEENlEXCEErDED

7 t'NO CALL SUBROUTINE RUNGE rO
IPDATE TANK PRESII
U

AND VOLUMES

Ass

STOP

'I

% %.

Ntm

!w*
.
P~*~ . ~V P. ~ V %~VV' 16~~ V
SUBRCUTINE CTANK F DMN1MADIIILZ

CALCUL.ATE FROM TANK I TO JUNCTION A

VIA PIPE SECTION 1

. CALCULATE PROM TNK 2 TO JUNCION A

VIA PIPE SECTION 2

CHK
,I
STAGNATION PRIES.SURES FROM

PI ECT.IONS I AND 2 AT JUNCTION A

AOJUST FLOW 21 VIA NO

SUBROUTINE RAIO IO PRESSURES CONVERGE ?

YES

CALCULATE FROM JUNCTION A TO

JUNCTION 8 VIA PIPE SECTION 4 USING
FLOW (4?. FLOW (1). FLOW f2)

C
___fCALCULATE FROM TANK 3 TO JUNCTION 61
VIA PIPE SECTION 3

CHECK STAGNATION PRESSURES FROM

PIPE SECTIONS 3 ANO 4 AT JUNCTION I

ADJUSTI
LTO VA NO PRESSURES CONVERGE ? ,7

POCTANK * STAGNATION PRES3URE AT

JUNCTION a j

1401 THIS F.OWCHART REPRESENTS 11 1-1

MOST COMPLE CASE IN WHICH PIPE
SECTIONS I THROUGH 3 ARE ALL
PRLENT. CASES WHERE NOT ALL RETRN.
SECTIONS ARE PRESENT ARE
AUTOMATICALLY HANDLED BY THE
PROGRAM, WPICH COMPENSATES
POI THE MISSINO SECTIONS. AND
VmlIES FLOW CONTINUITY w
CHOKING LOGIC. WHICH IS OETAILO,
IN THE BODY OF THIS REPORT ANO
INTHE COMPUTER PROGRAM ITSELF. -

HAS SEEN SUPPIIES3ED 1N THIS

FLOWCHART.

A- 2

%. %* V

'i _ NA
I~ . fl I,
,. M '-J
l LI N,_1 , .V: ,: -X ,,:. . , . , :. , , .. ; , y v'-.y v:. . ., .v- w- " xvy ~*y

SUaROUTINE SKVALVE
"P aES PIPE L.NGTH REPRESENT A ;PP VALVE

OR A VALVE ?

JA' CALCULATE OUTLET STTI PRESSURE.

AT THE OUTLET CALCULATE RATIO

R~rC.LT. .0 ?CALCULATE VALVE AREA

2SOTMA 4JUS
-CALCULATE PCASSUME UNCHIOKED CONDITION:
INLET STATIC PRESSURE
(PSTATIUN)

i - CALCULATE MACH2 BASEDI CALCULATE iNLET STATIC PRESSURE FCRI

ON RATIO
OMAC2 CHOKED CASE (PSTATICH,

CALCULJATE A FqCICM FACTOR SAS ON T ?

I T
OUTLET CONOITIONS Y S

HODVALVE -7 l uMKE VALVE

ETMT AM STATI PSAIT7 1 ATI .PT~U

I I ,
USING MACHI CALCULATE INLET FRICTION CALCULATE INLET STAGNATION PRESSURE
T T
FACTOR P01 P02 (PSTAT1 PS A 21 :.e

USING MACHI CALCULATE PSTATI

'pI

USING MACHI AMC0 PSTATrI CAL.CUL.ATE AMAN

.5

A-3
. :ip
SUBROUTINE REVERSE 01MIpNSCON ASO INITIALIZE

~TI~g *VIA SUBROUTINE SAKVALVE

RAT4, I - FLOW(IILCWCTANK

FLO4121FLW(1-FLW )A CI14 OMA 01411T!.0W(101 -

VIAIOS
11 S N

t!
t'It IICU I
VIA
PlPI SECTION 12

FROM&K1 JUNCTION
LW )7TI CHECK STAGNATION PRESSURES FROM
, PIPE SECTIONS 11 AN 12 AT JUNCTION F

YES

CALCULATE FR4OMJUNCTION F TO
JUNCTION E VIA PIPE SECTION 9

CALCULATE P
TANK To
UROM
JUAOTCTAN
O4
PIPE SI IOECTUON 1,14

?+ DO PRESSURES CONVERGE
W00H SAOUTINE RATIO

JU4TO Q
CACLT
FRGOM ',

JUNCT'ION I VIA PIPE[ SECT ION 10 (IL.Wfll - A"Ties- FLOWCTA S- I

PIPIE SECTIOQNSI AND 10 AT JUNCTIC14 E WITH SUBIROUTINE %IArlO +

YES

00 PRES[.SURE.S. CONV!IRGE[

J,-

Il

N N
'.

SUBROUTINE REVERSE (CONT'O)

CALCULATE FROM JUNCTION j TO

JUNCTION 0 VIA PlIt SECTION 8

RATIOS . FLOWIIPLOWCTANK

CAL UATE FROM JUNCTION 0 TO

JUNCTION C VIA PIPE SECTION S

FLW04LWCTAmx1 CALCULATE FROM JUNCTION 0 To

(I JUNCTION C VIA PIPE SECTION 7

RATIO FLOWCTANKI
FSLOW

AOJUST RATIOS CHECX STAGNATION PRESSURES FROM

WITH SUBROUTINE RTO PIPE SECTIONS 6 ANO 7 AT JUNCT:ON C

oo PESURE CON V.ERGE'

NO2I

CALCULATE FROM JUNCTION C TO

JUNCION 8 VIA PIPE SECTION 5

POIEVERSE * STAGNATION RRESSURE

OF Pipe SECTION S AT JUNCTION 8

7'p
RETURN p,

NOTE; THIS FLOWCHART REPRESENTS THE %

MOST COMPLEX CASE IN WHIC"
PIPE SECTIONS 5 THROUGH is ARE
ALL. PRESENT. CASES WHERE NOT
Al PI SECTIONS ARE PRESENT ARE
AUTOMATICALLY NANOLED BY THE
PROGRAM. WHICH COMPENSATES
FOR THE MISING SECTIONS ANO
VEIFaIlS FLOW CONTINUITY.

A-5 1
SUBROUTINE PAES

CALCULATE RATIO - WA
BASED ON MOOT

No RATIO. LT.1.0 N

N
A.
. XU
B ASED ON RAIO
CALLR.ATE FLOWM,

iSUPEAC)4KED INE 1
RAI L.A~
1.0,
PLMACAX
MO1.T
o

1 J /rEVEN IF=RAI"IO.!EQ.1"0 AT THE INL'ET

UNCHOKEO INLET; FRIcTION WIL' CAUSE SUPER-CHQKED

1LO . MOOT CONOITION: SWITCH = 2: FLO . FLOWMAX

.C MACHI SASID
LCULAT ESTrMATE FRICTION FACTOR

CA RATIO,1411 FROM MACH1 AND FLIO

I tALCUL.ATE AN (Nt.Er MACH NUMBdER

(MACHIMA4 THAT WILL. JUST CAtUSE THIE

IPE TO0 CHOKE AiT THE EXI

.............. Z...FROM
tSE MACHIMAX TO CALCULAITE ART
U 1
," ' C RRSMNOING FRICTION FACTOR

YES
-'5

,SP' ' 'S9 *..S 1 -' *

5.-6
Yes.

i .%

SUBROUTINE PRES (CONT"O

,%,'

CACULATE STATIC PROPRTIES AT INLET

BASED ON MACHIMAX"

CALCULATE THE MAXIMUM FLOW THE PiPE

CAN HANOLI (FLOWMAXI

L [iO.MOOT. T LWMAX" SWITCH -2

YESYE
No F) MACH2 . 1.0

UNCOKEMOOT. .FWMAX C FLO . FLOWMAX

CALCULATE MACHZ2BASED CHOKED EXIT; SWITC M I CALCULATE EXIT STATIC

ON MACHI AO PIN LENGTH
ONMC1ADPP MACH2 - 1.0 EGHANO BASED ON FLO
PROPERTIES MACM2"2 -

'5

FLO MOOT CAXLCULATE P02

I--

%
J'

I
,.

* J.

.- !S
55
I.r
Ky

K
SUIOUTINE FVALVE

F CALCULATE
VALVE AREA

CALZULAT1E INLAT STATIC 7%opEvvnEssdl

C~OKEDO'.W III
CACLT THEl.rl :

MOOIM10.
FIOWMMAX
DT.1T.OWMAX7hQTL-F.7-AX

SuPCOKD SITUATION C CON UNCHONOION

ITERATE VALVE IOUAMONS

SWITCH 2SWITCTO CALCULATE PSTAT2

CALCULATE PSTAT2 "ROM

CHOKING PRESSURE RAflO
,I

MACM2 1.0 F 2 0P0A. PSTAT1-PSTATZ)

RA PS
P l."PTAT1-VA"T2) I i PTO . P2 A T AT"

IASEI
CALCULATE MACHM
ON PRATIO

A-

RrI'UR 1 .5

A- - ,/

I,:
:'..

%
APPENDIX 3
SAMPLE INPUT DATA FILE .1J

SA4.-5
.40.0"" N75TmCR
2O3. .0,2 3 83 .*0,2C 3*C O.0,, 1 .0.0 , 10.0, I40."
320.,320.,320. ,200., 00..200.,200.
1,1,1,1,1,1,1
2,2,2, 2,2, ,2-'
3,3,3,2,2,2,2
2 , Z22 2, 29292
12. 13 . 14 .
.111,. J8,.070
13 ,4 1 Z12
.035,.081,.077
14., 16., 15.
.075,.081,.077
13., -39
.033, .111
15. 16.
.125,. 105
14.,16.
.0999.0C85
13. ,15.
• I 5,.166
15. "12
.175,.185
10. ,13.
.095 , . 098
12z. , 3.
.079, .083
13., -9 3.
.071,.083
11., 1.
.09, 083
12. .13.
•11 ,.125
12. . 3.
.102, . 101
I.4, 53.3

B-1/B-2
Reverse Blank %
APPENDIX C
SAIMFL.E OUTPUT DAT7A FPILE1

L%,U T 1ATA FIltENAME * INPUT

1
STAGNATION ?!0,!qirUE*f 540.0000 Of.305!s uutmN!
SD!CItIt HOAT WA?3 A 1.4000
;A 'C STA'T*0 I 330 FTOGAE AmxIN 5
7A14 04TA

24!SSURE 2 750.00 LS/INSS*Z VOLUP! 3Z.00 FT*03

TANK~ I I STA41ATICN

PRESSURE! 2750.00 L!/14**Z v3LURE * 32.10 9TO83

%2* 2STA.AT3:4

1 750.00 L!/?NW2 VLUFS 32.20 9 **3

'44K a 3 STA;AT!0Q4 PIESSURE

1AI~~V. I ~STAGNATION PRISSURS 1 25.00 LB/IN**z VLM 00 ~'

AGNAT:.2N '05551,9! 140.01 L5/IN**2 VJ LU!* 3*S
'A4K a IzZ
I 225.12 LS/:NjsaZ VM 22.00 1*
OCT
3
1*l 576 daTION 00!SSLZE

STAGNAT:O24 PQ!SSURE * 140.00 LS/I'4sa2 Y3LUME 22.30 FT3

EllA a 1'

ol A! S!CTIh I I. PIPE a I L!NGTH a z2.0 003. AT IAMETE4 1 0.1110 AT

2:3!S IOCTI Z? a 1 pipe 0 2 L !N 'TI " 13.0000 AT CIAaE'0 0.0900 AT%
P!P! SECT: ;b L PP 0 3 L ." T?* 14.0000 AT OIAMFTEt 0.0700 IT
PIPE 0 LZN&,TII 13.0000 IT C!A4ETEA I 0.0050 IT
2:9i SC?::F4 a z
210! S CT::N I z 313! a 2 LENGTM 14.0000 AT CrAm!tfl 0 .0!10 IT'
'P2f S1CT::1 * 2 PIPE0! 3- z,;TN a Z.0300 27 A4!E 0.070
SIPE S!CTI~h a 3 PIPE # I. L!NGTM 14.0000 £T OI14tTE.4 1 0.0750 AT
S!CtIM, a 3 pipe 0 2 LINGTH 16.0000 Fr7 0I A04eT Et 0.0310 a'
'1P f 0.077' AT
3 PIP! A 3 LE4GTma 15.0000 FT OI A tTEA *
0!01 S!C TI IN w
I. LZNSTm 13.0000 FT CIAMETE4 0.054;
2i's Ss!::h a 4 PIP! a
PIPE S!CT oh a 4 aloe a z LuNGTH * 10.0300 '7 0 IA -E T 4 a 0.1110 '

1 L;'-3T 4 a 15.U000 2r 02A Mf E a 0.1250 AT

012! sc T :14 9 5 0PP *
P~ae S!r:*a 5 pip! a Z IS A I0'4TACL OeVICt
3 L!MG7TM 0 15.0000 AT ClAmETi.~ 0.0)S0 '1T
SLOE sac'zpq 1 5 PIP! 0
0IAA E TE~ A 0.0392 I
212s ssCTZ"-% a 1 PtPE 6 t. LIN;T4 1 14.0000 1T7
Z LsNGT4 I t6.0a00I OIA1E
C1 T Ex 0.0450 9T
It'! 551712'. a s PIPE A
0.0135
LzN*;T'i 13.0000 FT IT
210!
a
S1'!1CT1'.h 7
r~
pS
A't~t I
!~Ma
a
1.00'' ~ IIA-ETE~A
.. TV ~ 016
SET' 10'
I! 2 L 54NST4 a 15.0000 4 T I,A E.! 0.1750 IT
I PIPE I I L!N;TM a 15.0000 AT C1A"V1.1 0.2.050 '7
3:91 5OCTt2'. a Ot~A S4T EA a 0. 150 I7
311!I IP 0 1 L NGT'. - 10.0000 7
SI"1? 9;ST PIPE
01!4 1 13.0000 £7 O:Aa575d I 0.09!3 £T
1 prof 0 1 L MIT'4 - 12.0000 FT 2tlAmeE. 0 0.0140 *
2101 -2
3073'

it PIPE 0 1 L!.1P-4 0 13.0000 17 0~!! 0.0710 IT

2:2! sict::. 6
'1'! a 417' L=IN'Tw 19.0000 A CtAme7!d 0.331 I
51::'. 0
SECT 12 P:Pj 0 1 L!Wm 1 11.0:00 F1' 2:4-fT!51 I a.2450 £T

'po -CT:Zl a 21P! 0 2 Llm- TN4 - 11.0300 '7CIT'! 0.0763 'T

'p~a 1
3ic7::N '~z4a 1*- * 1.0000 '7T A!! .. i
':'0~~~~ sO:Na- 35a 2L.,P,4 Z .0002 C7 A 0.1Z51 '7
':'! ~ a2'"! ~ ~a 2
~ 0C'. L '.Na
L31T'l
1210 7
9.:003 IT
I'!
::A.!,
I
.
.122 -r
j'.'
.:12~

C-1~

A. IN lw_ %-
APPENDIX C
SAMPLE OUTPUT DATA F:LE:

* SAOL! NA?.C-7K

'NOU' :4?1 LIST%4G

C. RU1r 4?A FtLfhAmi INPUT

STAGNATION T!PPER1?U4e - 540.0000 OE RfIS RANKINE

SD!CI'FIC N;Ar 35thr 1 1.4000
;AS 03 4ST1;T .13.3000 FT/0EGR4!!S RANKIN!

I 2710.00 LB/IN**Z VOLUpe 32.00 FT*w3

IA~ SA;NAIC4E 20SSE
!7!0.00 L!/INSZ * 32.30 Ir**3
* 2 SA5%A?::4 PsOSSU2E

1 750.Q0 L!/14*Z V2U!* 3!10 I*Wa3

'A4K a 3 STANAT'0 PRESSURE

StAAGNtrIO01 ;R!SSUR! 22.00o La/IN4*.z YCLUME 21300 FT*a3

fAN * 11

I 225.CO L9/:N*-Z VOLUMS 200 AT,-3

,4x 0 11 STA;NAT:.4 OCISSLIE

1 40.00 LS/14**2 Y3L'J0E * 230 I*w3

flIK a 14 SA3-NAT7.- PQOSSU2E

a %. t.NG' 1 12.0003 IT DjAmETEA 1 0.1110 IT

Plot sIcTIZPS I P ips 5
I z LgNGTI4 - L3.3000 AT ClA E'EA 0.0900 AT
2:2 SiCTI,% a P1!PE
PP 0 3 L NGTH4 A 14.0000 AT 2 1AM E T!.4 0.0700*T
A!P! 500Z4 * 8
A 0.GaS00 T
%~~ PPps
PI 1 LSN',T4 - 13.0000 AT ttA15430
S.'
'0!1CI I!U 2L6NGT4 A 14.CC00 'T ~IA-!TE( A 0.3!10 IT
':iS~LI I0a 3L1N;4 I 12.01cQ X? al~veTf.( 1 0.0770 91
*IF C~.* 3 pipe a I L;!NGTM 14.0000 A7 osAAElE.4 1 0.0110 AT
22 7.CTIIK a 3 *Pp 4 -1 LENGTH 16.0000 AT O,4!T!, A 0.0310 AT
3 LENGTH4 * S.0000 FT CAE! .73F
0:01 S;CTI.M 0 3 PIPE a
Pips a 1 LENGTH 13.0000 PT CIA-ET;!4 * .0
21F! S;CTI-' a 4
PI' g DIOS:: 4~ LOIdGTo 13.0300 'T 0IA -ETfi* 0113
'3'!~~.... ao 07C1 J4 IAF T ex
20SCt^,kf a S DISa -- IO.1V~
a 3L5NGT7.4 * 15.0000 FT C3A~37ij * 0.0)80 17
'lot SeOT:ha 9 3 PIP
aT
E5~~
DI PI' !jNIF .4 S 14.0000 AT I1A pe E4 O .0990 IT
a
S, !C"!F
-:3T a600
a a., '7jAE! 0.0050 --
LSN-l4 16.0000 FT CIh!Td 0.35 *
212E StC71k a 6 PIPE 'A 1
iLz.437T-* 13.0000 'T ZIA-eTSA - 0-550 IT
2I2! SICl.2, A I or;'! a
I L!N;fl4 15.0000 IT 0!Aw!,E4( - 0.1750 AT
3:2f SECT::h 0 3 DIP! 0
AIDE a L ;45T 4 1 12.0000 FT C IAF1TE.4 1 0.1!50 cy
2!-I'! 5503114 I 0
;457-4 4 10.0000 '7 O1A'!7!' I C.0150 AT
2t3! SiCT I:N a 9 PIPE 0 1
ISa 2L;4744 13.2100 AT *:ME 0.0490 AT
9
0C3'a
I L 5NGT. tZ1.0000 FT I~~! 3.0790 '
* '*S SSCT*.]N a II PIPE a
I LSN4TH 1.3000 zT 3~!! 0.03!0 '
210! SIC~TI; a 13 ;Plot a
it PIPE a I. L *T1 13.2000 AT OZAw!3!4 1 0.0710 27
I 2
'I'S $1Ct: 6
0 .3131 IT
15 a 1?
3.34ps 0 * 11.0000 AT CZAET!4
'0'! 7
SC3!M
Pips t I .9434 * 1.o300 27 3:1~s .3653
1:2 SECT:I4 0 12
~:~ sot::
a IL I'!* L1'N.?N
;T
'N7:
* 11.000 :7
12.00 AT
::F!.
15S a
0.0781 AT
.~
': C 1! '.!a.47
a 2 13.30001AT ::A." I .:5
' S' 07 .I.
L * Z..12co* r~!4* 312
':' :0::a*'''
a~~ 3*:s01 XT''

Lo5
ZESULTS AT !NO CP TIME STEP 0 a
%M 0.00000 SEC

qESULTS AT ENO CO 0105 SECTIONS I TO 1'

*roe if I MAC a 0.0338 PLO0W KATE 7.4,336 LEIR/SEC STAG%6TIU4 POISSURE 2522.7 LSF/IN**Z
PIPE a 2 MAC" a 0.1211 PLOW RAT! 1.03S7 LaiP/SEC STAGNAT.JI, PRESSURE 2523.0 LSF/1N32*
a:'Ia 3 14ACIN a 0.1631 SLO3W RATE 1.6173 LSN,'SEC S716AATI.it. PRESSURE 2Z56.) L!6/:4**2
vIPE1a 4 MA6CH a 0.1257 PLOW RATE IS.5.13O L14/ SEIC STA G'4AT I i. PRtESSURc Z252.2 LPF/zJ''2 .
31'!9 a 5 4ACM 6 0.3894 CLOW RATE 25.3663 L 9 4zs EC ST G'4AT 1,. PRESSURS 1439.3 1.E':NSM2 %6
SIDE 6 MACH a 0. 3910 PLOW RAT! 4.5945 Le N/SEC ST AGNA TI-l. PRESSUR! 154S.3 L!P/INSMZ 1
PIPE 4 ? MACH a 0.t035 PLOW RITE 20.49ZI LGAMVSEC STIGNAtlt. PRESSURE 15695 L!FPINO'Z
&IF! a a MACH 0 0.106? PLOW RATE 25.0865 LBO/SEC STAGNATI.~ PRESSURE 1544.3 L!F/140MZ
SIDE a I MACH 0.3633 PLIW KAI 1Z.3246 LIN/SEC STAGNAT.31. PA IS SURE 1038.2 L 2/IN'1**2
012! *to MACH A .5600 GLOW RATE 12.1620 LON/SEC STAGNATji, PRESSURE 769.1 1.SP/lN*S2
0191 all "Arm 4 1.3000 a1.3W ITEI 6.5106 LEN/SEC ST IGN T~td. PRESSJRE 225.0 LSP/14**Z
PIPE 012 MAl)4 a 1.1000 FLOW RITE 5.3140 LBE'/SEC STAGNAXiji. PRESSURE_ 1413.a L!0/IN*0Z
'12! 413 '96CM a D.2169 GLOW RATE 1.3251 LIM/S2C STAGNATIJt. PRESSURE Z25.0 LSP/XN*
PIPE *14 MACH a 1.0000 FLO3W RAT! 10.93149 L3MN/SIC STAGNdITl. PRESSURE 140.0 L!F/I.*Z

ISSULTS AT E4C 09 TI-E STEP 0

rre - O.S0G30 SEC

!A4NC RESULTS

rAmK a I STA;NATION PRESSURE * 2726.76 CS/?MS*Z VOLUME - 32.0000 PT**3

TANK 0 2 STA~h4ION PRESSURE . 2124.90 LS/IN*82 VOLUME - 32.0000 07ss3
TANK. a 3 STAGNATION PRESSURE s Z?19.96 LS/1NSO2 VOLUME - 32.0000 PT**l
TANK 2 11 STAZNATION PRESSURE . 257.53 L3/XM.** VOLUME - 21.0000 F?*03
r~l 0 12 STA NAT1'9 PRSSSURE . 163.05 L3/tI&**2 VOLUME - 20.0000 FT**3
T a.11 4 1! STA MATI3N PRESSURE a 234.12 LS/tN**Z VOLUle - 20.Q000 FT**3
7AK 14 STAGNATZCON PRESSURE - 194.65 LB/IN8S2 VOLUN4e * 000 7*

RESULTS AT ViC0!O9 PIPE SECTIONS I To 14

PIPE I 1 MACH C 0.1360 PLOW RATE 7.3712 L3M/SEC STAGNAT1jh PRESSURE 2500.9 LSF,'1N**Z
,tee * 2 MACH a 0.1213 PLOW CATE 7.;759 Lf'9/SIC STIGN*TTUjN PRESSURE 2499.1 L3F//N*02
ItR 8 1 MACH a 0.16410 PLOW RATE 9.9456 LSM/SE STIGNAT~j" PRESSURE Z21.3 LSP/1%N*Z
'I'E 1 11M
.1260 GLOW C ITS I S.31'61 LEN/SEC ST4GNA T.J,. PRE SSUR 1 2231.9 L!F/,4wS
'0'! a S MACH a 0.!394 FL:W CI?! 24.3986 LBM/SEC STIGNATIJ,. PRESSURE liZ7.1 L!F/I4saw2le
'EPEI 1 5 MACM a 0.0900 FLOW RATE 4.3603 LOM/SEC STAGNAI... PRESSUCE 1573.7 L!/r.'d*2
'IRE 0 7 MACH a 0.1055 FLOW CITE 20.3395 LBM/S EC ST IGNA TIJI. PRESSURE 1573.7 L!F/'N0*Z
'1'a -ACM a 0.1067 GLOW RIT! 24.J936 L314/SEC STAGNAT.jI. PRESSURE 15!3.3 L'9/'4I'2
sI'I9 ACM a 0.3633 &L:W RAT! 12.232! LS'/SEC ST 4GNATI.;. PRESSURE 1050.4 L Vv/ t'1*a
3:2! *to WacoN a 0.3400 FLOW1 RATE 12.i67 LB"/S!C STIGNATI.4. PRESSURE 7S3.3 L'IINSW2Z
21e 011 M6CM A 1.0000 GLOW RATE 6.4621 L!M/SEC ST A444U TI;. RESSUR% 257.3 L''IINIWgZ
AID! 612 MAC). 6 1.0010 4L.3w DATE 5.1707 LEN/SEC ST&GN6TI1jf. 'RESSURS 161.1 L 3F /1'4*02
2 :P 10a MACH S 0.1115 PLCW RATE 1.3113 LE1N/ SE C SThCMAT.. -72ZSSUC! 234.L Ll N/ 14*2
'1! 614 MACM a 1. ac 0 'L3%l RAT! 10. 1555 LEMl/S !C STIGNAT.Jl. *RESSURi 194.7 L!F/t'4a'2

610,
4!SULTS AT EN CF T"MS STEP 2

TA'4K I I $TAGNIT!OR =SISLqE - 2'O3. 74 LS/INS*Z V:2. . 32000C

12 01TW 3
TANK 6 2 SA-INATICd 3q!SSURE - 2 703,00 LI/IN** VCLI!" - 32.3C00 FT*A.?
tA'iA a 3 STA;,MA?:C4 09 3SSUlf 2690.15 L3/:)N3*2 YCI.A- 32. 000 FT*13%
r& K AI sA;NAT
S Ti :C4 §RSSURE - 23 9. S2 L S/ tm.. VCLU'46 - .0. 10Z2C '7883
T A-iK a 12 srA;NATI.3 002SSU S . 197.99 Ll/ IN**! Vci~u'q - 3. JcoC O7a83
TANK( a 13 STAhMAT3N -~3IQ
-S2E - 2 43. 17 L!/:Nv*Z VcLu4S - :0.;000 FT*-!
TANK A 14 STA;N&7'tCN ORiOSSURE - 24i.49 L3/ IN*sZ VOLumE - 23 .30^3 FT*13

VISULTS AT E'0 GO SI S!C?3:%S I tC 14

3P ~ i 1.136a 7L~8T
.35.9 La-/SIC 3A:.P1su0 27. LF/Z2
3!zz 0 2 "AC~ a . 12.2 cL3W 2 VS 7. 1475 UM/SEC S13IN 7: , 3.6
; RESSU1. S247 L!'/:sZ
*A' 0* .16 31
J.i~ *LO- 2371 3. 40 5 L~m/ SEC STAG5'4 ESUR 203.1 L if/' N*-
a .cMq
V!~ a 0.1270 FLI3- ATE 15.ZS34 L 3 4/S 5C ST AG A'LAE. ~Sjq
2235. 3 L ! /CPl*S2Z.
:'
*S ACH P .32)4 GLZ3 RAT: 24..7048 L/SC ST3GNAT I..l.OR 5SURI5 1. L!F/1H*92
3:Ea 6 AC94a .3.3910 FL:Z' DATE 5 I?
'.31 S/S!C STI7.Nh Trr. 04ESSURE 15591..4 L97/1NO*,
Ps! a T -8C1M a 3.loss FL~ LB7v0..,
SE/3C Sl AG'i*Tl~i. P#3 SS j; 1561.4. I/14**2
a AM* . 1367
I:~ cLZW RAT; Z4.704(8 L!-/SEC ST AG 4AT :, Aq;SSu#1 152 1.3 Lic/:4-Z
** - "A-). i 0. !63! GLZ'. CAT3 12: Ll-/SEC STAG,AT., OR;j9S SF/hJQ
E' t0 0I 1.56, 2* OAT I7 . 547
1Z L!-/SEC ST A'1 A I.;. *223Sups ??5.7 Li/!4*02
A:2
S kI -AI: a 1 . 000 FLIW a ATE S.-Z37 L!-/SEC STA' A', 02ES.9! 2 1 %*
't'e2
E2'
-12 "AZM a 1. 3000 .FL'. RATS 5.7364 L! m., c ST &GNA - '23 R S. 197. 3 L ''Nsaz
A13
S12 ma.: it 0. 1362 FL:.1 4 %T 1.7940 L3 / S EC STAG NA T j. 2 q .5313 14 3.2 L6 a2
''1 814 a~~,. 1.30 FL:W #37 10.7!17 L B64/ 3EC S~GAI, QSU0 24 3.9 '/~

C-3/C- 4
Reverse 3Jlank