**AIAA 98-3682**

**A GENERALIZED FLUID SYSTEM SIMULATION PROGRAM TO MODEL FLOW
DISTRIBUTION IN FLUID NETWORKS**

**A.K. Majumdar**

**J.W. Bailey**

**P.A. Schallhorn**

T.  **Steadman**

**Sverdrup Technology, Inc.**

**Huntsville, AL 35806**

**A GENERALIZED FLUID SYSTEM SIMULATION PROGRAM TO MODEL FLOW
DISTRIBUTION IN FLUID NETWORKS**[^1]

Alok Majumdar, John W. Bailey, Paul Schallhorn and Todd Steadman

Sverdrup Technology

Huntsville, Alabama - 35806

***ABSTRACT***

*This paper describes a general purpose computer program for analyzing
steady state and transient flow in a complex network. The program is
capable of modeling phase changes, compressibility, mixture
thermodynamics and external body forces such as gravity and centrifugal.
The program's preprocessor allows the user to interactively develop a
fluid network simulation consisting of nodes and branches. Mass, energy
and specie conservation equations are solved at the nodes; the momentum
conservation equations are solved in the branches. The program contains
subroutines for computing "real fluid" thermodynamic and thermophysical
properties for 33 fluids. Seventeen different resistance/source options
are provided for modeling momentum sources or sinks in the branches. The
system of equations describing the fluid network is solved by a hybrid
numerical method that is a combination of the Newton-Raphson and
successive substitution methods. This paper also illustrates the
application and verification of the code by comparison with the Hardy
Cross method for steady state flow and an analytical solution for
unsteady flow.*

# INTRODUCTION

A fluid flow network consists of a group of flow branches, such as pipes
and ducts, that are joined together at a number of nodes. They can range
from simple systems consisting of a few nodes and branches to very
complex networks containing many flow branches simulating valves,
orifices, bends, pumps and turbines. In the analysis of existing or
proposed networks, some node pressures and temperatures are specified or
known. The problem is to determine all unknown nodal pressures,
temperatures and branch flow rates.

The oldest method for systematically solving a problem consisting of
steady flow in a pipe network is the Hardy Cross method \[1\]. Not only
is this method suited for solutions generated by hand, but it has also
been widely employed for use in computer generated solutions. But as
computers allowed much larger networks to be analyzed, it became
apparent that the convergence of the Hardy Cross method might be very
slow or even fail to provide a solution in some cases. The main reason
for this numerical difficulty is that the Hardy Cross method does not
solve the system of equations simultaneously. It considers a portion of
the flow network to determine the continuity and momentum errors. The
head loss and the flow rates are corrected and then it proceeds to an
adjacent portion of the circuit. This process is continued until the
whole circuit is completed. This sequence of operations is repeated
until the continuity and momentum errors are minimized. It is evident
that the Hardy Cross method belongs in the category of successive
substitution methods and it is likely that it may encounter convergence
difficulties for large circuits. In later years, the Newton-Raphson
method has been utilized \[2\] to solve large networks, and with
improvements in algorithms based on the Newton-Raphson method, computer
storage requirements are not much larger than those needed by the Hardy
Cross method.

A Generalized Fluid System Simulation Program (GFSSP) for modeling
cryogenic fluids in a complex flow circuit was developed \[3,4\]. The
objective of the present paper is to describe further development and
validation of GFSSP to model unsteady flow.

#  MATHEMATICAL FORMULATION

GFSSP assumes a Newtonian, steady, non-reacting and one dimensional flow
in the flow circuit. The flow can be either laminar or turbulent,
incompressible or compressible, with or without heat transfer, phase
change and mixing. The analysis of the flow and pressure distribution in
a complex fluid flow network requires resolution of the system into
nodes and branches. Nodes can be either boundary nodes or internal
nodes. Pressures, temperatures, and concentrations of fluid species are
specified at the boundary nodes. At each internal node, scalar
properties such as pressures, temperatures, enthalpies, and mixture
concentrations are computed. The flow rates (vector properties) are
computed at the branches.

##  GOVERNING EQUATIONS

Figure 1 displays a schematic showing adjacent nodes, their connecting
branches, and the indexing system used by GFSSP. In order to solve for
the unknown variables, mass, energy and fluid specie conservation
equations are written for each internal node and flow rate equations are
written for each branch.

### Mass Conservation Equation 

![](media/image1.emf) (1)

Equation 1 requires that, for the transient formulation, the net mass
flow from a given node must equate to rate of change of mass in the
control volume. In the steady state formulation, the left hand side of
the equation is zero. This implies that the total mass flow rate into a
node is equal to the total mass flow rate out of the node.

### Momentum Conservation Equation 

The flow rate in a branch is calculated from the momentum conservation
equation (Equation 2) which represents the balance of fluid forces
acting on a given branch. A typical branch configuration is shown in
Figure 2. Inertia, pressure, gravity, friction and centrifugal forces
are considered in the conservation equation. In addition to these five
forces, a source term S has been provided in the equation to input pump
characteristics or to input power to a pump in a given branch. If a pump
is located in a given branch, all other forces except pressure are set
to zero. The source term, S, is set to zero in all branches without a
pump.

![](media/image2.emf) ![](media/image3.emf){width="3.2618055555555556in"
height="3.4375in"}

**Figure 1** **- Schematic of GFSSP** **Nodes, Branches and Indexing
Practice**

![](media/image4.emf){width="2.9923611111111112in"
height="3.542361111111111in"}

**Figure** **2** **- Schematic of a Branch Showing Gravity and
Rotation**

$\begin{array}{r}
\frac{\left( mu_{t + Dt} - mu_{t} \right)}{g_{c}Dt} + \frac{\overset{\text{.}}{m_{\text{ij}}}}{g_{c}}\left( u_{\text{ij}} - u_{u} \right) = \\
 \\
\left( p_{i} - p_{j} \right)A + \frac{\rho\text{gV}\text{cos}q}{g_{c}} - K_{f}\overset{\text{.}}{m_{\text{ij}}}\left| \overset{\text{.}}{m_{\text{ij}}} \right|A + \frac{\rho K_{\text{rot}}^{2}w^{2}A}{2g_{c}}\left( r_{j}^{2} - r_{i}^{2} \right) + S
\end{array}$ ( 2)

![](media/image5.emf)

The two terms in the left hand side of the momentum equation represent
the inertia of the fluid. The first one is the time dependent term and
must be considered for unsteady calculations. The second term is
significant when there is a large change in area or density from branch
to branch. The first term in the right hand side of the momentum
equation represents the pressure gradient in the branch. The pressures
are located at the upstream and downstream face of a branch. The second
term represents the effect of gravity. The gravity vector makes an angle
(θ) with the flow direction vector. The third term represents the
frictional effect. Friction was modeled as a product of K~f~ and the
square of the flow rate and area. K~f~ is a function of the fluid
density in the branch and the nature of the flow passage being modeled
by the branch. The calculation of K~f~ for different types of flow
passages has been described in detail later within this report. The
fourth term in the momentum equation represents the effect of the
centrifugal force. This term will be present only when the branch is
rotating as shown in Figure 2. K~rot~ is the factor representing the
fluid rotation. K~rot~ is unity when the fluid and the surrounding solid
surface rotates with the same speed. This term also requires a knowledge
of the distances between the upstream and downstream faces of the branch
from the axis of rotation.

### Energy Conservation Equation 

The energy conservation equation for node i, shown in Figure 1, can be
expressed mathematically as shown in Equation 3.

![](media/image6.emf) (3)

![](media/image7.emf) (4)

Equation 3 shows that for transient flow, the rate of increase of
internal energy in the control volume is equal to the rate of energy
transport into the control volume minus the rate of energy transport
from the control volume plus the rate of work done on the fluid by the
pressure force plus the rate of work done on the fluid by the viscous
force plus the rate of heat transfer into the control volume.

For a steady state situation, the energy conservation equation, Equation
3, states that the net energy flow from a given node must equate to
zero. In other words, the total energy leaving a node is equal to the
total energy coming into the node from neighboring nodes and from any
external heat sources (Q~i~) coming into the node and work done on the
fluid by pressure and viscous forces. The MAX operator used in Equation
3 is known as an upwind differencing scheme which has been extensively
employed in the numerical solution of Navier-Stokes equations in
convective heat transfer and fluid flow applications \[5\]. When the
flow direction is not known, this operator allows the transport of
energy only from its upstream neighbor. In other words, the upstream
neighbor influences its downstream neighbor but not vice versa. The
second term in the right hand side represents the work done on the fluid
by the pressure and viscous force. The difference between the steady and
unsteady formulation lies in the left hand side of the equation. For a
steady state situation, the left hand side of Equation 3 is zero, where
as in unsteady cases the left hand side of the equation must be
evaluated.

### Fluid Specie Conservation Equation 

The flow network shown in Figure 1 has a fluid mixture flowing in most
of the branches. In order to calculate the density of the mixture, the
concentration of the individual fluid species within the branch must be
determined. The concentration for the k^th^ specie can be written as
shown in Equation 4.

For a transient flow, Equation 4 states that the rate of increase of the
concentration of the k^th^ specie in the control volume equals the rate
of transport of the k^th^ specie into the control volume minus the rate
of transport of the k^th^ specie out of the control volume.

Like Equation 3, for steady state conditions, Equation 4 requires that
the net mass flow of the k^th^ specie from a given node must equate to
zero. In other words, the total mass flow rate of the given specie into
a node is equal to the total mass flow rate of the same specie out of
that node. For steady state, the left hand side of Equation 4 is zero.
For the unsteady formulation, the resident mass in the control volume is
changing and therefore, the left hand side must be computed.

**Thermodynamic and Thermophysical Properties**

The momentum conservation equation, Equation 2, requires knowledge of
the density and the viscosity of the fluid within the branch. These
properties are functions of the temperatures, pressures and
concentrations of fluid species for a mixture. Three thermodynamic
property routines have been integrated into the program to provide the
required fluid property data. GASP \[6\] provides the thermodynamic and
transport properties for ten fluids. These fluids include Hydrogen,
Oxygen, Helium, Nitrogen, Methane, Carbon Dioxide, Carbon Monoxide,
Argon, Neon and Fluorine. WASP \[7\] provides the thermodynamic and
transport properties for water and steam. GASPAK \[8\] provides
thermodynamic properties for helium, methane, neon, nitrogen, carbon
monoxide, oxygen, argon, carbon dioxide, hydrogen, parahydrogen, water,
isobutane, butane, deuterium, ethane, ethylene, hydrogen sulfide,
krypton, propane, xenon, R-11, R-12, R-22, R-32, R-123, R-124, R-125,
R-134A, R-152A, nitrogen trifluoride and ammonia. For RP-1 fuel, a look
up table of properties has been generated by a modified version of GASP.
An interpolation routine has been developed to extract the required
properties from the tabulated data.

### Friction Calculations 

### 

It was mentioned earlier in this paper that the friction term in the
momentum equation is expressed as a product of K~f~, the square of the
flow rate and the flow area. Empirical information is necessary to
estimate K~f~ . Several options for flow passage resistance are as
follows: pipe flow, flow through a restriction, non-circular duct, pipe
flow with entrance and/or exit losses, thin sharp orifice, thick
orifice, square edge reduction, square edge expansion, rotating annular
duct, rotating radial duct, labyrinth seal, parallel plates, common
fittings and valves, pump characteristics, pump power, valve with a
given loss coefficient, and a Joule-Thompson device

COMPUTER PROGRAM

GFSSP was developed on an IBM compatible PC using the Lahey EM32 FORTRAN
compiler. The same source code also runs on a Silicon Graphics
workstation. The code was developed with a modular structure to achieve
two objectives. First, the code's solver module was separated from the
preprocessor module such that users are not required to write any code
to develop their models. Secondly, the code can easily be extended to
enhance its capability by adding new modules. The main routine controls
all of the program operations and makes the decisions whether to
continue or to stop the calculations.

The computer program has three major parts. The first part consists of
the subroutines for the preprocessor. The preprocessor allows the user
to interactively create the flow network model consisting of nodes and
branches. All of the input specifications, including the boundary
conditions, are specified through the preprocessor. The second major
part of the program consists of the subroutines that provide the initial
conditions and then develop and solve all of the conservation equations
in the flow network. The third part of the program consists of the
thermodynamic property programs, GASP, WASP, and GASPAK that provide the
necessary thermodynamic and thermophysical property data required to
solve the resulting system of equations.

Figure 3 shows GFSSP's process flow diagram. The user runs the
interactive preprocessor to generate the input data file. The input data
file contains all the information necessary for the model. The solver
module reads the input data file and produces the solution in
conjunction with the thermodynamic property programs. It should be noted
that the user interfaces with the program only through the input data
file. The user is never required to modify the source program to develop
the model. The user is also not required to understand the structure of
the solver module in order to develop their model. However, for the
completeness of the documentation, a flow chart showing the major
activities of the code is described in Figure 4.

![](media/image8.emf){width="3.3027777777777776in"
height="1.832638888888889in"}

**Figure 3 -** **GFSSP** **Process Flow Diagram**

![](media/image9.emf){width="4.663194444444445in"
height="5.829166666666667in"}

![](media/image10.emf){width="6.483333333333333in"
height="8.400694444444444in"} ![](media/image11.emf) **Figure 4 -
GFSSP** **Flowchart of Major Subroutines**

**RESULTS**

The paper also describes comparison of steady state flow network
simulation with the predictions of the Hardy Cross method and unsteady
flow simulation of a tank blow down with an analytical solution. A
typical water distribution network is shown in Figure 5 and Table 1.

Figure 6 shows a comparison between GFSSP and Hardy Cross predicted
flowrates. The comparison appears reasonable considering the fact that
Hardy Cross method assumes a constant friction factor in the branch
while GFSSP computes the friction factor for each branch during every
iteration. Therefore, as the flowrates change the friction factor also
changes.

![](media/image12.emf){width="2.4319444444444445in"
height="1.7986111111111112in"}

**Figure 5 - Water Distribution Network Schematic**

**Table 1 - Water Distribution Network Branch Data**

+----------+-------------------+----------------------+--------------+
| Branch   | Length(inches)    | Diameter(inches)     | Roughness    |
|          |                   |                      |              |
|          |                   |                      | Factor       |
+----------+-------------------+----------------------+--------------+
| 12       | 120               | 6                    | 0.0018       |
+----------+-------------------+----------------------+--------------+
| 25       | 2400              | 6                    | 0.0018       |
+----------+-------------------+----------------------+--------------+
| 27       | 2400              | 5                    | 0.0018       |
+----------+-------------------+----------------------+--------------+
| 57       | 1440              | 4                    | 0.0018       |
+----------+-------------------+----------------------+--------------+
| 53       | 120               | 5                    | 0.0018       |
+----------+-------------------+----------------------+--------------+
| 56       | 2400              | 4                    | 0.0018       |
+----------+-------------------+----------------------+--------------+
| 64       | 120               | 4                    | 0.0018       |
+----------+-------------------+----------------------+--------------+
| 68       | 1440              | 4                    | 0.0018       |
+----------+-------------------+----------------------+--------------+
| 78       | 2400              | 4                    | 0.0018       |
+----------+-------------------+----------------------+--------------+
| 89       | 120               | 5                    | 0.0018       |
+----------+-------------------+----------------------+--------------+

![](media/image13.emf){width="6.492361111111111in"
height="1.8736111111111111in"}

**Figure 6 - Predicted Flow Rate Comparison Between GFSSP**

**and the Hardy Cross Method**

![](media/image14.emf)

> **(a) (b)**

**Figure 7 - Physical Schematic (a) and GFSSP Model (b) for Venting
Nitrogen Tank**

[Blowdown of a Tank]{.underline}

Figure 7 shows a tank with an internal volume of 10 ft3 containing
nitrogen gas at a pressure and temperature of 100 psia and 80 ^°^F
respectively. The nitrogen is discharged into the atmosphere through an
orifice with a 0.1 inch diameter until the pressure in the tank becomes
50 psia.

[Analytical Solution]{.underline}:

The differential equation governing an isentropic blow down process can
be written as:

$\left( \frac{p}{p_{i}} \right)^{(1 - 3g\frac{)}{2}g}\frac{d\left( \frac{p}{p_{i}} \right)}{dt}\text{ = }\frac{gA}{\rho_{i}V}\sqrt{gg_{c}p_{i}\rho_{i}}\left( \frac{2}{g + 1} \right)^{(g + 1\frac{)}{2}(g - 1)}$
(5)

This is an initial value problem and the initial conditions are:

$$t\text{ = }0,\frac{p}{p_{i}} = 1$$

The analytical solution for $\frac{p}{p_{i}}$ is given by Moody \[9\]
as:

$\frac{p}{p_{i}} = \left\lbrack 1 + \left( \frac{\lambda - 1}{2} \right)\left( \frac{2}{g + 1} \right)^{(g + 1\frac{)}{2}(g - 1)}\sqrt{\frac{gg_{c}p_{i}}{\rho_{i}}}\frac{At}{V} \right\rbrack^{- 2\frac{g}{(}g - 1)}$
(6)

The analytical and GFSSP solutions are compared in Figure 8. Excellent
agreement was observed between the two solutions.

![](media/image15.emf){width="3.5930555555555554in"
height="2.297222222222222in"}

**Figure 8 - Comparison of the Predicted Pressure History by GFSSP and
the Analytical Solution**

**CONCLUSIONS**

The Generalized Fluid System Simulation Program was further extended to
handle unsteady compressible flow in a complex flow circuit. An
additional thermodynamic property code was integrated to increase the
fluid data base from 12 to 33. The code also provides seventeen
resistance options for modeling momentum sources or sinks in a network.
The code includes a command line as well as a visual preprocessor to
facilitate the creation of the input data file. The predicted steady
state flow and pressure distribution in a network was verified by
comparing with the predictions with the Hardy Cross method. Unsteady
predictions for a blow down of a pressurized tank were compared with
analytical predictions. Both comparisons were satisfactory.

**REFERENCES**

1\. Hardy Cross, "Analysis of Flow in Networks of Conduits or
Conductors", Univ. Ill. Bull. 286, November 1936

2.  Jeppson, Ronald W., "Analysis of Flow in Pipe Networks", Ann Arbor
    Science

3.  Majumdar, A. K. and VanHooser, K. P., "A General Fluid System
    Simulation Program to Model Secondary Flows in Turbomachinery",
    Paper No. AIAA 95-2969, 31^st^ AIAA /ASME /SAE /ASEE Joint
    Propulsion Conference, July 10-12, 1995, San Diego, California.

4.  Majumdar, A. K., Bailey, J. W. and Sarkar, B., "A Generalized Fluid
    System Simulation Program to Model Flow Distribution in Fluid
    Networks", Paper No. AIAA 97-3225, 33^rd^ AIAA /ASME /SAE /ASEE
    Joint Propulsion Conference, July 6- 9, 1997, Seattle, WA.

5.  Patankar, S. V., "Numerical Heat Transfer and Fluid Flow",
    Hemisphere Publishing Corp., Washington, D. C., 1980.

6\. Hendricks, R. C., Baron, A. K., and Peller, I. C., "GASP - A
Computer Code for Calculating the Thermodynamic and Transport Properties
for Ten Fluids: Parahydrogen, Helium, Neon, Methane, Nitrogen, Carbon
Monoxide, Oxygen, Fluorine, Argon, and Carbon Dioxide", NASA TN D-7808,
February, 1975.

7\. Hendricks, R. C., Peller, I. C., and Baron, A. K., "WASP - A
Flexible Fortran IV Computer Code for Calculating Water and Steam
Properties", NASA TN D-7391, November, 1973.

8\. Cryodata Inc., "User\'s Guide to GASPAK, Version 3.20", November,
1994.

9\. Moody, F. J., "Introduction to

Unsteady Thermofluid Mechanics",

John Wiley, 1990.

[^1]: ^ã^Copyright ^^ 1998 by the American Institute of

    Aeronautics and Astronautics, Inc.

    All rights reserved.