National Aeronautics and

Space Administration

IS20

George C. Marshall Space Flight Center

Huntsville, Alabama 35812
The NASA STI Program...in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientific and Technical Information (STI) 

Program Office plays a key part in helping NASA 

maintain this important role.

The NASA STI Program Office is operated by 

Langley Research Center, the lead center for 

NASA’s scientific and technical information. The 

NASA STI Program Office provides access to 

the NASA STI Database, the largest collection of 

aeronautical and space science STI in the world. 

The Program Office is also NASA’s institutional 

mechanism for disseminating the results of its 

research and development activities. These results 

are published by NASA in the NASA STI Report 

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant 

phase of research that present the results of 

NASA programs and include extensive data 

or theoretical analysis. Includes compilations 

of significant scientific and technical data 

and information deemed to be of continuing 

reference value. NASA’s counterpart of peerreviewed formal professional papers but has less 

stringent limitations on manuscript length and 

extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis.

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical conferences, 

symposia, seminars, or other meetings sponsored 

or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientific, technical, 

or historical information from NASA programs, 

projects, and mission, often concerned with 

subjects having substantial public interest.

• TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission.

Specialized services that complement the STI 

Program Office’s diverse offerings include creating 

custom thesauri, building customized databases, 

organizing and publishing research results...even 

providing videos.

For more information about the NASA STI Program 

Office, see the following:

• Access the NASA STI program home page at 

<http://www.sti.nasa.gov>

• E-mail your question via the Internet to 

<help@sti.nasa.gov>

• Fax your question to the NASA STI Help Desk 

at 443 –757–5803

• Phone the NASA STI Help Desk at 

443 –757–5802

• Write to:

NASA STI Help Desk

NASA Center for AeroSpace Information

7115 Standard Drive

Hanover, MD 21076–1320
i

NASA/TM—2011–216470

Generalized Fluid System Simulation Program,

Version 5.0—Educational

A.K. Majumdar

Marshall Space Flight Center, Huntsville, Alabama

August 2011

National Aeronautics and

Space Administration

Marshall Space Flight Center • Huntsville, Alabama 35812
ii

Available from:

NASA Center for AeroSpace Information

7115 Standard Drive

Hanover, MD 21076 –1320

443–757–5802

This report is also available in electronic form at

<https://www2.sti.nasa.gov/login/wt/>

Acknowledgments

Generalized Fluid System Simulation Program (GFSSP) development has been supported by several Marshall 

Space Flight Center projects that include Fastrac Engine, Space Station, Space Shuttle, Ares I, and Nuclear Propulsion. The contributions of John Bailey, Andre LeClair, Ric Moore, Paul Schallhorn, Todd Steadman, and Saif Warsi 

are gratefully acknowledged. The author would also like to thank Melissa Vandyke, Stan Tieman, and Sammy Nabors 

for their constant encouragement and support for continuing GFSSP development. 
iii

PREFACE

The motivation to develop a general purpose computer program to compute pressure and 

flow distribution in a complex fluid network came from the need to calculate the axial load on the 

bearings in a turbopump. During the early years of space shuttle main engine development, several 

specific purpose codes were developed to model the turbopumps. However, it was difficult to use 

those codes for a new design without making extensive changes in the original code. Such efforts 

often turn out to be time consuming and inefficient. To satisfy the need to model these turbopumps 

in an efficient and timely manner, development of the Generalized Fluid System Simulation Program (GFSSP) was started at NASA Marshall Space Flight Center (MSFC) in March 1994. The 

objective was to develop a general fluid flow system solver capable of handling phase change, compressibility, and mixture thermodynamics. Emphasis was given to construct a user-friendly program 

using a modular structured code. The intent of this effort was that an engineer with an undergraduate background in fluid mechanics and thermodynamics should be able to rapidly develop a reliable 

model. The interest in modular code development was intended to facilitate future modifications 

to the program. 

The code development was carried out in several phases. At the end of each phase, a workshop was held where the latest version of the code was released to MSFC engineers for testing, 

verification, and feedback. The steady state version of GFSSP (version 1.4) was first released in 

October 1996. This version is also commercially available through the Open Channel Foundation. The unsteady version was released in October 1997 (version 2.0). A graphical user interface 

(GUI) for GFSSP was developed and was part of Version 3.0 that was released in November 1999. 

GFSSP (version 3.0) won the NASA Software of the Year award in 2001. Fluid transient (water 

hammer) capability was added in version 4.0 that was released in March 2003. The main highlight 

of this version (version 5.0) is its capability to handle conjugate heat transfer.

This Technical Memorandum (TM) provides a detailed discussion of the mathematical 

formulation, computer program, GUI, and includes a number of example problems. Section 1 

provides an introduction and overview of the code. The mathematical formulation, which includes 

the description of governing equations and the solution procedure to solve these equations, is 

described in section 2. The program structure is discussed in section 3. Section 4 describes GFSSP’s 

GUI, which is called VTASC (visual thermofluid dynamic analyzer for systems and components). 

Several example problems are described in section 5. The new user may skip sections 2 and 3 initially, but will benefit from these sections after gaining some experience with the code.

This TM has been prepared for the educational version of GFSSP for distribution to 

accredited United States universities who have signed an agreement with MSFC for use in their 

educational curriculum. This version has been limited to developing models using only 30 or fewer 

nodes without any provision of user subroutines.
iv
v

TABLE OF CONTENTS

1. INTRODUCTION ........................................................................................................... 1

1.1 Network Flow Analysis Methods ............................................................................... 2

1.1.1 Network Definitions ........................................................................................ 3

1.2 Units and Sign Conventions ....................................................................................... 3

1.3 Data Structure ........................................................................................................... 4

1.4 Mathematical Formulation ......................................................................................... 5

1.5 Fluid Properties .......................................................................................................... 7

1.6 Flow Resistances ......................................................................................................... 8

1.7 Program Structure ...................................................................................................... 9

1.8 Graphical User Interface ............................................................................................ 10

1.9 Example Problems ...................................................................................................... 11

2. MATHEMATICAL FORMULATION ............................................................................ 12

2.1 Governing Equations .................................................................................................. 12

2.1.1 Mass Conservation Equation .......................................................................... 13

2.1.2 Momentum Conservation Equation ................................................................ 13

2.1.3 Energy Conservation Equations ...................................................................... 16

2.1.3.1 Energy Conservation Equation of Fluid ............................................ 16

2.1.3.2 Energy Conservation Equation of Solid ............................................ 17

2.1.4 Fluid Specie Conservation Equation ............................................................... 20

2.1.5 Thermodynamic and Thermophysical Properties ............................................ 20

2.1.5.1 Equation of State for a Real Fluid ..................................................... 21

2.1.6 Mixture Property Calculations ........................................................................ 21

2.1.7 Friction Calculations ....................................................................................... 23

2.1.7.1 Pipe Flow (Branch Option 1) ............................................................. 23

2.1.7.2 Flow Through a Restriction (Branch Option 2) ................................. 23

2.1.7.3 Noncircular Duct (Branch Option 3) ................................................. 24

2.1.7.4 Pipe With Entrance and Exit Loss (Branch Option 4) ....................... 25

2.1.7.5 Thin, Sharp Orifice (Branch Option 5) .............................................. 26

2.1.7.6 Thick Orifice (Branch Option 6) ........................................................ 27

2.1.7.7 Square Reduction (Branch Option 7) ................................................. 28

2.1.7.8 Square Expansion (Branch Option 8) ................................................ 29

2.1.7.9 Rotating Annular Duct (Branch Option 9) ........................................ 29

2.1.7.10 Rotating Radial Duct (Branch Option 10) ......................................... 30

2.1.7.11 Labyrinth Seal (Branch Option 11) .................................................... 31

2.1.7.12 Flow Between Parallel Plates (Branch Option 12) .............................. 32

2.1.7.13 Common Fittings and Valves (Branch Option 13) ............................. 33

2.1.7.14 Pump Characteristics (Branch Option 14) ......................................... 34

2.1.7.15 Pump Horsepower (Branch Option 15) .............................................. 35
vi

TABLE OF CONTENTS (Continued)

2.1.7.16 Valve With a Given Loss Coefficient (Branch Option 16) .................. 35

2.1.7.17 Joule-Thompson Device (Branch Option 17) ..................................... 35

2.1.7.18 Control Valve (Branch Option 18) ..................................................... 35

2.1.7.19 User-Defined Resistance (Branch Option 19) ..................................... 36

2.1.7.20 Heat Exchanger Core (Branch Option 20) ......................................... 36

2.1.7.21 Parallel Tube (Branch Option 21) ...................................................... 37

2.1.7.22 Compressible Orifice (Branch Option 22) .......................................... 37

2.2 Solution Procedure ..................................................................................................... 38

3. COMPUTER PROGRAM ................................................................................................ 42

3.1 Process Flow Diagram ................................................................................................ 43

3.2 Solver and Property Module ....................................................................................... 43

3.2.1 Nonsimultaneous Solution Scheme ................................................................. 46

3.2.2 Simultaneous Solution Scheme ........................................................................ 46

3.2.3 Conjugate Heat Transfer ................................................................................. 46

3.2.4 Thermodynamic Property Package .................................................................. 49

4. GRAPHICAL USER INTERFACE ................................................................................. 50

4.1 Menus ......................................................................................................................... 51

4.1.1 File Menu ........................................................................................................ 51

4.1.2 Edit Menu ....................................................................................................... 52

4.1.3 Advanced Menu .............................................................................................. 52

4.1.4 Run and Module Menus .................................................................................. 53

4.1.5 Display, Canvas, and Help Menus ................................................................... 53

4.2 Global Options ........................................................................................................... 53

4.2.1 General Information ....................................................................................... 53

4.2.2 Circuit Options ................................................................................................ 56

4.2.3 Unsteady Options ............................................................................................ 58

4.2.4 Fluid Options .................................................................................................. 58

4.3 Fluid Circuit Design ................................................................................................... 60

4.3.1 Boundary and Internal Node Properties ......................................................... 60

4.3.2 Branch Properties ............................................................................................ 61

4.3.3 Conjugate Heat Transfer ................................................................................. 67

4.4 Advanced Options ...................................................................................................... 70

4.4.1 Transient Heat ................................................................................................. 70

4.4.2 Heat Exchanger ............................................................................................... 70

4.4.3 Tank Pressurization ......................................................................................... 71

4.4.4 Turbopump ..................................................................................................... 72

4.4.5 Valve Open/Close ............................................................................................ 72

4.4.6 Fluid Conduction ............................................................................................ 73
vii

TABLE OF CONTENTS (Continued)

4.5 GFSSP Input File ....................................................................................................... 73

4.5.1 Title Information ............................................................................................. 73

4.5.2 Logical Variables ............................................................................................. 74

4.5.3 Node, Branch, and Fluid Information ............................................................. 74

4.5.4 Solution Control Variable ................................................................................ 74

4.5.5 Time Control Variables ................................................................................... 75

4.5.6 Fluid Designation ............................................................................................ 75

4.5.7 Node Numbering and Designation .................................................................. 76

4.5.8 Node Variables ................................................................................................ 76

4.5.9 Transient Heat/Variable Geometry Information .............................................. 77

4.5.10 Node-Branch Connections .............................................................................. 77

4.5.11 Branch Flow Designation and Resistance Options .......................................... 77

4.5.12 Unsteady Information ..................................................................................... 78

4.5.13 Inertia Information ......................................................................................... 78

4.5.14 Fluid Conduction Information ........................................................................ 79

4.5.15 Rotation Information ...................................................................................... 80

4.5.16 Valve Open/Close Information ........................................................................ 80

4.5.17 Momentum Source Information ...................................................................... 81

4.5.18 Heat Exchanger Information ........................................................................... 81

4.5.19 Moving Boundary Information ....................................................................... 81

4.5.20 Turbopump Information ................................................................................. 81

4.5.21 Tank Pressurization Information ..................................................................... 82

4.5.22 Variable Rotation Information ........................................................................ 82

4.5.23 Restart File ...................................................................................................... 82

4.5.24 Cyclic Boundary Information .......................................................................... 83

4.5.25 Conjugate Heat Transfer Information ............................................................. 83

4.6 GFSSP Execution ....................................................................................................... 84

4.6.1 Steady State Run Manager .............................................................................. 84

4.6.2 Unsteady Run Manager .................................................................................. 85

4.7 GFSSP Output File .................................................................................................... 86

4.7.1 Titles and Data Files ....................................................................................... 86

4.7.2 Logical Variables ............................................................................................. 86

4.7.3 Node and Branch Information ........................................................................ 87

4.7.4 Fluid Information ........................................................................................... 87

4.7.5 Boundary Conditions ...................................................................................... 88

4.7.6 Fluid Network Information ............................................................................. 88

4.7.7 Initial Field Information ................................................................................. 88

4.7.8 Conjugate Heat Transfer Network Information .............................................. 89

4.7.9 Solution Results .............................................................................................. 90

4.7.10 Convergence Information ................................................................................ 93

4.8 Post-Processing Simulation Data ............................................................................... 93

4.8.1 Steady State Simulation Results ...................................................................... 94
viii

TABLE OF CONTENTS (Continued)

4.8.2 Unsteady Simulation Results ........................................................................... 94

4.8.2.1 VTASC Plot ....................................................................................... 94

4.8.2.2 Winplot .............................................................................................. 97

4.8.3 Display in Flow Circuit ................................................................................... 97

5. EXAMPLES ...................................................................................................................... 100

5.1 Example 1—Simulation of a Flow System Consisting of a Pump, Valve, 

and Pipe Line .............................................................................................................. 100

5.1.1 Problem Considered ........................................................................................ 100

5.1.2 GFSSP Model ................................................................................................. 101

5.1.2.1 Branch 12 (Pump) .............................................................................. 101

5.1.2.2 Branch 23 (Gate Valve) ....................................................................... 104

5.1.2.3 Branch 34 ( Pipeline) .......................................................................... 104

5.1.3 Results ............................................................................................................ 104

5.2 Example 2—Simulation of a Water Distribution Network ......................................... 106

5.2.1 Problem Considered ........................................................................................ 106

5.2.2 GFSSP Model ................................................................................................. 106

5.2.3 Results ............................................................................................................ 108

5.3 Example 3—Simulation of Compressible Flow in a Converging-Diverging Nozzle .... 109

5.3.1 Problem Considered ........................................................................................ 109

5.3.2 GFSSP Model ................................................................................................. 109

5.3.3 Results ........................................................................................................... 111

5.4 Example 4—Simulation of the Mixing of Combustion Gases and a Cold 

 Gas Stream ................................................................................................................. 115

5.4.1 Problem Considered ........................................................................................ 115

5.4.2 GFSSP Model ................................................................................................. 115

5.4.3 Results ............................................................................................................. 116

5.5 Example 5—Simulation of a Flow System Involving a Heat Exchanger ..................... 117

5.5.1 Problem Considered ........................................................................................ 117

5.5.2 GFSSP Model ................................................................................................. 118

5.5.2.1 Branches 12, 34, 56, and 78 (Pipe Lines) ............................................ 118

5.5.2.2 Branches 23 and 67 (Heat Exchanger) ............................................... 119

5.5.3 Results ........................................................................................................... 120

5.6 Example 6—Radial Flow on a Rotating Radial Disk .................................................. 122

5.6.1 Problem Considered ........................................................................................ 122

5.6.2 GFSSP Model ................................................................................................. 123

5.6.3 Results ........................................................................................................... 123

5.7 Example 7—Flow in a Long-Bearing Squeeze Film Damper ...................................... 125

5.7.1 Problem Considered ........................................................................................ 125

5.7.2 GFSSP Model ................................................................................................. 126

5.7.3 Results ........................................................................................................... 127
ix

TABLE OF CONTENTS (Continued)

5.8 Example 8—Simulation of the Blowdown of a Pressurized Tank ............................... 130

5.8.1 Problem Considered ........................................................................................ 130

5.8.2 GFSSP Model ................................................................................................. 130

5.8.3 Results ........................................................................................................... 131

5.8.3.1 Analytical Solution ............................................................................ 131

5.9 Example 9—A Reciprocating Piston-Cylinder ............................................................ 132

5.9.1 Problem Considered ........................................................................................ 132

5.9.2 GFSSP Model ................................................................................................. 132

5.9.3 Results ........................................................................................................... 134

5.10 Example 10—Power Balancing of a Turbopump Assembly ........................................ 135

5.10.1 Problem Considered ........................................................................................ 135

5.10.2 GFSSP Model ................................................................................................. 136

5.10.3 Results ........................................................................................................... 139

5.10.4 Parametric Study ............................................................................................. 140

5.11 Example 11—Steady State and Transient Conduction Through a Circular Rod, 

With Convection ......................................................................................................... 142

5.11.1 Problem Considered ........................................................................................ 142

5.11.2 GFSSP Model ................................................................................................. 142

5.11.3 Results ........................................................................................................... 144

5.11.3.1 Analytical Solution ............................................................................ 144

5.12 Example 12—Simulation of Fluid Transient Following Sudden Valve Closure ........... 145

5.12.1 Problem Considered ........................................................................................ 145

5.12.2 GFSSP Model ................................................................................................. 146

5.12.3 Results ............................................................................................................. 147

5.13 Publications ................................................................................................................ 148

APPENDIX A—DERIVATION OF Kf FOR PIPE FLOW .................................................... 149

APPENDIX B—SUCCESSIVE SUBSTITUTION METHOD OF SOLVING COUPLED 

NONLINEAR SYSTEMS OF ALGEBRAIC EQUATIONS ...................... 150

APPENDIX C—NEWTON-RAPHSON METHOD OF SOLVING COUPLED 

NONLINEAR SYSTEMS OF ALGEBRAIC EQUATIONS ...................... 151

APPENDICES D–P (Supplemental Information on CD Inside Back Cover)
x

LIST OF FIGURES

 1. A typical flow network consists of fluid node, solid node, flow 

branches, and conductors ......................................................................................... 4

 2. Data structure of the fluid-solid network has six major elements .............................. 5

 3. Schematic of mathematical closure of GFSSP .......................................................... 6

 4. Interpropellant seal flow circuit in a rocket engine turbopump ................................. 7

 5. GFSSP’s program structure showing the interaction of three major modules ........... 10

 6. GFSSP’s graphical user interface, VTASC allows creating, running, 

and viewing results in one environment .................................................................... 11

 7. Schematic of GFSSP nodes, branches, and indexing practice ................................... 12

 8. Schematic of a branch showing gravity and rotation ................................................ 14

 9. A schematic showing the connection of a solid node with neighboring solid, 

fluid, and ambient nodes ........................................................................................... 18

10. Pipe resistance option parameters ............................................................................. 23

11. Noncircular duct cross section: (a) Rectangle, (b) ellipse, (c) concentric

cylinder, and (d) circular sector ................................................................................. 24

12. Pipe with entrance and/or exit loss resistance option parameters .............................. 26

13. Thin, sharp orifice resistance option parameters ....................................................... 26

14. Thick orifice resistance option parameters ................................................................ 27

15. Square reduction resistance option parameters ......................................................... 28

16. Square expansion resistance option parameters ........................................................ 29

17. Rotating annular duct resistance option parameters ................................................. 30

18. Rotating radial duct resistance option parameters .................................................... 31

19. Labyrinth seal resistance option parameters ............................................................. 32
xi

LIST OF FIGURES (Continued)

20. Parallel flat plate resistance option parameters ......................................................... 33

21. Control valve in a pressurization system ................................................................... 36

22. Heat exchanger core ................................................................................................. 37

23. Parallel tube .............................................................................................................. 38

24. SASS (simultaneous adjustment with successive substitution) scheme 

for solving governing equations ................................................................................ 40

25. GFSSP process flow diagram showing interaction among three modules ................. 43

26. Flowchart of nonsimultaneous solution algorithm in solver 

and property module ................................................................................................ 44

27. Flowchart of simultaneous solution algorithm in solver 

and property module ................................................................................................ 45

28. Solid temperature calculation by successive substitution method .............................. 47

29. Solid temperature calculation by Newton-Raphson method ..................................... 48

30. Calculation of residuals of energy conservation equation 

for Newton-Raphson method ................................................................................... 48

31. Main VTASC window .............................................................................................. 50

32. VTASC file menu ...................................................................................................... 51

33. VTASC edit menu ..................................................................................................... 52

34. VTASC advanced menu ............................................................................................ 52

35. Global Options dialog .............................................................................................. 54

36. General Information dialogs: (a) User Information tab, 

(b) Solution Control tab, and (c) Output Control tab ............................................... 55

37. Circuit Options dialogs: (a) Circuit tab and (b) Initial Values tab ............................. 57

38. Unsteady Options dialog .......................................................................................... 58
xii

LIST OF FIGURES (Continued)

39. Fluid Options dialog ................................................................................................. 59

40. Boundary and internal nodes on canvas ................................................................... 60

41. Node pop-up menus: (a) Boundary node and (b) internal node) .............................. 61

42. Node Properties dialog ............................................................................................. 62

43. Nodes with branch ‘handles’ ..................................................................................... 62

44. Direct line segment branch: (a) Step 1 and (b) step 2 ............................................... 63

45. Two-line segment branch: (a) Step 1 and (b) step 2 .................................................. 63

46. Example fluid circuit with complete branch connections .......................................... 64

47. Branch pop-up menu ................................................................................................ 64

48. Relocate Branch ID dialog ........................................................................................ 65

49. Change Branch Connections dialog .......................................................................... 65

50. Branch Resistance Options dialog ............................................................................ 66

51. Example fluid circuit with resistance options ............................................................ 66

52. Fittings and valve resistance option properties dialog ............................................... 67

53. Rotation/Momentum dialog ..................................................................................... 67

54. Solid Node Properties dialog .................................................................................... 68

55. Ambient Node Properties dialog .............................................................................. 68

56. Conductor pop-up menu .......................................................................................... 69

57. Conductors dialog .................................................................................................... 69

58. Transient Heat Load option dialog ........................................................................... 70

59. Heat Exchanger option dialog .................................................................................. 71
xiii

LIST OF FIGURES (Continued)

60. Tank Pressurization option dialog ............................................................................ 71

61. Turbopump option dialog ......................................................................................... 72

62. Valve Open/Close option dialog ................................................................................ 72

63. Fluid conduction dialog ............................................................................................ 73

64. GFSSP steady state run manager .............................................................................. 85

65. GFSSP unsteady run manager .................................................................................. 85

66. GFSSP steady state simulation results internal fluid node table ................................ 95

67. GFSSP results dialog for unsteady simulation .......................................................... 95

68. GFSSP VTASC plot properties dialog: (a) Data tab, (b) labeling tab, 

and (c) scale tab ........................................................................................................ 96

69. Display Results in flow circuit example ..................................................................... 98

70. Display Results/Properties dialog ............................................................................. 98

71. Display Property Units dialog .................................................................................. 99

72. Schematic of pumping system and reservoirs (Example 1) ........................................ 102

73. Manufacturer supplied pump head-flow characteristics ............................................ 102

74. GFSSP model of pumping system and reservoirs: (a) Detailed schematic 

and (b) VTASC model .............................................................................................. 103

75. Pump characteristics curve in GFSSP format ........................................................... 104

76. Fluid system operating point .................................................................................... 106

77. Water distribution network schematic (Example 2) ................................................... 107

78. GFSSP model of the water distribution network....................................................... 107

79. A flow rate comparison between GFSSP and hardy cross method predictions ......... 109
xiv

LIST OF FIGURES (Continued)

80. Converging-diverging steam nozzle schematic (Example 3) ...................................... 110

81. Converging-diverging steam nozzle model: (a) Detailed schematic 

and (b) VTASC model............................................................................................... 110

82. Predicted pressures for the isentropic steam nozzle ................................................... 112

83. Predicted temperatures for the isentropic steam nozzle ............................................. 112

84. Temperature/entropy plot comparing the isentropic steam nozzle 

with an irreversible process ....................................................................................... 115

85. Mixing problem schematic (Example 4) .................................................................... 116

86. GFSSP model of mixing problem ............................................................................. 116

87. Flow system schematic of a heat exchanger (Example 5) .......................................... 119

88. GFSSP model of the heat exchanger: (a) Detailed schematic 

and (b) VTASC model .............................................................................................. 119

89. VTASC Heat Exchanger dialog ................................................................................ 120

90. Temperature and flow rate predictions in heat exchanger .......................................... 121

91. Flow schematic of a rotating radial disk (example 6) ................................................ 123

92. GFSSP model of the rotating radial disk: (a) Detailed schematic 

and (b) VTASC model .............................................................................................. 124

93. Comparison of GFSSP model results with experimental data ................................. 125

94. Squeeze film damper schematic (example 7, view 1) .................................................. 126

95. Squeeze film damper schematic (example 7, view 2) .................................................. 126

96. Unwrapping and discretization of squeeze film damper ........................................... 127

97. GFSSP model of squeeze film damper ..................................................................... 127

98. Predicted circumferential pressure distributions in the squeeze 

film damper ............................................................................................................. 129
xv

LIST OF FIGURES (Continued)

 99. Comparison of GFSSP model results with experimental data 

for squeeze film damper ............................................................................................ 129

100. Venting air tank schematics: (a) Physical schematic, (b) detailed 

model schematic, and (c) VTASC model ................................................................... 131

101. Comparison of the predicted pressure history by GFSSP 

and the analytical solution ........................................................................................ 133

102. Piston-cylinder configuration .................................................................................... 133

103. Coordinate transformed piston-cylinder configuration ............................................. 134

104. GFSSP model of the piston-cylinder: (a) Detailed model schematic 

and (b) VTASC model .............................................................................................. 134

105. Predicted temperature history of piston-cylinder model ........................................... 135

106. Predicted pressure history of piston-cylinder model ................................................. 136

107. Simplified turbopump assembly ................................................................................ 136

108. A flow circuit containing turbopump assembly: (a) Detailed model 

schematic and (b) VTASC model .............................................................................. 137

109. Example 10 turbopump dialog ................................................................................. 139

110. Example 10 heat exchanger dialog: (a) Heat exchanger 1 and (b) heat exchanger 2 ... 139

111. GFSSP RCS model results ........................................................................................ 140

112. Parametric study results: turbopump pressure differential ........................................ 141

113. Parametric study results: turbopump hydrogen mass flow rate ................................. 141

114. Parametric study results: turbopump torque and horsepower ................................... 142

115. Schematic of circular rod connected to walls at different temperatures ..................... 143

116. GFSSP model of circular rod for example 11: (a) Detailed model schematic 

and (b) VTASC model .............................................................................................. 144
xvi

LIST OF FIGURES (Continued)

117. Comparison of GFSSP temperature prediction and closed-form solution ................ 145

118. Schematic of a propellant tank, pipeline, and valve .................................................. 146

119. GFSSP model of a propellant tank, pipeline, and valve ............................................ 146

120. Example 15 Valve Open/Close dialog ....................................................................... 147

121. Comparison of GFSSP and MOC predicted pressure oscillations ............................ 148
xvii

LIST OF TABLES

 1. Units of variables in input/output and solver module ............................................... 4

 2. Mathematical closure ................................................................................................ 6

 3. Fluids available in GASP and WASP ........................................................................ 8

 4. Fluids available in GASPAK .................................................................................... 8

 5. Resistance options in GFSSP ................................................................................... 9

 6. Poiseuille number coefficients for noncircular duct cross sections ............................. 25

 7. Constants for two-K method of hooper for fittings/valves (GFSSP 

resistance option 13) ................................................................................................. 34

 8. Mathematical closure ................................................................................................ 38

 9. Winplot comma delimited unsteady output files ....................................................... 97

10. Use of various options in example problems ............................................................ 101

11. Tabulated pump characteristics data ......................................................................... 104

12. Predicted system characteristics ................................................................................ 105

13. Water distribution network branch data ................................................................... 107

14. GFSSP predicted pressure distribution at the internal nodes .................................... 108

15. GFSSP and Hardy Cross method predicted branch flow rates .................................. 108

16. Converging-diverging nozzle branch information ..................................................... 111

17. Converging-diverging nozzle boundary conditions ................................................... 111

18. Predicted mass flow rate with varying exit pressure ................................................... 113

19. Comparison of choked mass flow rates ..................................................................... 114

20. Branch dimensions of squeeze film damper .............................................................. 128

21. Moving boundary information of squeeze film damper ............................................ 128
xviii

LIST OF ACRONYMS

CFD computational fluid dynamics

FVM finite volume method

GFSSP Generalized Fluid System Simulation Program

GUI graphical user interface

LOX liquid oxygen

MOC method of characteristics

MSFC Marshall Space Flight Center

SASS simultaneous adjustment with successive substitution

SIMPLE semi-implicit pressure linked equation

SP solver and property

TM Technical Memorandum

TPA turbopump assembly

US user subroutine

VTASC visual thermofluid dynamic analyzer for systems and components
xix

NOMENCLATURE

A area (in2)

A0 pump characteristic curve coefficient

a length (in)

B0 pump characteristic curve coefficient

b length (in)

C heat capacity (Btu/s-°R)

Cc convergence criterion

CL flow coefficient

Cp specific heat (Btu/lb °F) 

Cv flow coefficient for a valve

c clearance (in)

ci,k mass concentration of kth specie at ith node 

D diameter (in)

f Darcy friction factor

G bound mass (lb/ft2-s)

g gravitational acceleration (ft/s2); specific heat ratio

gc conversion constant (= 32.174 lb-ft/lbf-s2)

h enthalpy (Btu/lb)

hij heat transfer coefficient (Btu/ft2-s-°R)

J mechanical equivalent of heat (778 ft-lbf/Btu)

K, K1 nondimensional head loss factor 

Ke exit loss coefficient

Kf flow resistance coefficient (lbf-s2/(lb-ft)2 )

Krot nondimensional rotating flow resistance coefficient
xx

m

NOMENCLATURE (Continued)

Ki inlet loss coefficient

k thermal conductivity (Btu/ft-s-°R)

kv empirical factor

L length (in)

LΩ resistance of the Joule-Thompson device

M molecular weight

m resident mass (lb)

mass flow rate (lb/s)

mp pitch (in)

N revolutions per minute (rpm); number of iterations

NB number of branches

NE number of equations

NN number of nodes

Nu Nusselt number

n number of teeth

P pump power (hp)

Po Poiseuille number

Pr Prandtl number 

p pressure (lbf /in2)

Q, q heat source (Btu/s)

R gas constant (lbf-ft/lb-°R)

Re Reynolds number (Re = ρuD/m)

r radius (in)

S momentum source (lbf)

S entropy generation (Btu/s-°R)

s entropy (Btu/lb-°R)
xxi

NOMENCLATURE (Continued)

T fluid temperature (°F)

Ts solid temperature (°F)

u velocity (ft/s)

V volume (in3)

Vf viscosity correction factor

v specific volume (ft3/lb)

w Joule-Thomson device flow rate (lbm/hr)

x quality and mass fraction

xk molar concentrations

z compressibility factor

Greek

α multiplier for Labyrinth seal resistance; under-relaxation parameter

b ratio of orifice diameter to inside pipe diameter

g specific heat

∆h head loss (ft)

∆τ time step (s)

d distances between velocity locations (ft)

d ij distance between two solid nodes (ft)

ε absolute roughness (in), heat exchanger effectiveness, Labyrinth seal carryover factor

ε /D relative roughness

εij emissivity 

θ angle between branch flow velocity vector and gravity vector (deg); 

angle between neighboring branches for computing shear (deg)

θr time required to drain pressurized propellant tank (s)
xxii

NOMENCLATURE (Continued)

η efficiency 

μ viscosity (lb/ft-s)

n kinematic viscosity (ft2/s)

ρ density (lb/ft3)

ρ molar density (lb-mol/ft3)

s Stefan-Boltzman constant (= 4.7611×10–13 Btu/ft2-R4-s)

t time (s)

ω angular velocity (rad/s)

Subscript

amb ambient

B back

c cold

cr critical

Dis discharge 

eff effective

F front

f liquid

g vapor 

gen generation

h hot

Im impeller

i node; inner

ij branch
xxiii

NOMENCLATURE (Continued)

o outer

or orifice

rot rotational

S solid

s entropy

trans transverse

Turb turbine

u upstream branch
xxiv
1

TECHNICAL MEMORANDUM

GENERALIZED FLUID SYSTEM SIMULATION PROGRAM,

VERSION 5.0—EDUCATIONAL

1. INTRODUCTION

The need for a generalized computer program for thermofluid analysis in a flow network 

has been felt for a long time in Aerospace Industries. Designers of thermofluid systems often need 

to know pressures, temperatures, flow rates, concentrations, and heat transfer rates at different 

parts of a flow circuit for steady state or transient conditions. Such applications occur in propulsion systems for tank pressurization, internal flow analysis of rocket engine turbopumps, chilldown 

of cryogenic tanks and transfer lines, and many other applications of gas-liquid systems involving fluid transients and conjugate heat and mass transfer. Computer resource requirements to 

perform time-dependent, three-dimensional Navier-Stokes computational fluid dynamic (CFD) 

analysis of such systems are prohibitive and therefore are not practical. A possible recourse is to 

construct a fluid network consisting 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, node pressures, temperatures, and concentrations at the system boundaries are usually known. The problem is to determine all internal nodal pressures, temperatures, concentrations, and branch flow 

rates. Such schemes are known as Network Flow Analysis methods and they use mostly empirical 

information to model fluid friction and heat transfer. For example, an accurate prediction of axial 

thrust in a liquid rocket engine turbopump requires the modeling of fluid flow in a very complex 

network. Such a network involves the flow of cryogenic fluid through extremely narrow passages, 

flow between rotating and stationary surfaces, phase changes, mixing of fluids, and heat transfer. 

Propellant feed system designers are often required to analyze pressurization or blowdown processes in flow circuits consisting of many series and parallel flow branches containing various pipe 

fittings and valves using cryogenic fluids. The designers of a fluid system are also required to know 

the maximum pressure in the pipeline after sudden valve closure or opening. 

Available commercial codes are generally suitable for steady state, single-phase, incompressible flow. Because of the proprietary nature of such codes, it is not possible to extend their capability to satisfy the above-mentioned needs. In the past, specific purpose codes were developed to 

model the space shuttle main engine (SSME) turbopump. However, it was difficult to use those 

codes for a new design without making extensive changes in the original code. Such efforts often 

turn out to be time consuming and inefficient. Therefore, the Generalized Fluid System Simulation 

Program (GFSSP)1 has been developed at NASA Marshall Space Flight Center (MSFC) as a general fluid flow system solver capable of handling phase changes, compressibility, mixture thermodynamics, and transient operations. It also includes the capability to model external body forces such 

as gravity and centrifugal effects in a complex flow network. The objective of the present effort 
2

is to develop (1) a robust and efficient numerical algorithm to solve a system of equations describing a flow network containing phase changes, mixing, and rotation and (2) to implement the algorithm in a structured, easy-to-use computer program. 

This program requires that the flow network be resolved into nodes and branches. The program’s preprocessor allows the user to interactively develop a fluid network simulation consisting 

of fluid nodes and branches, solid nodes, and conductors. In each branch, the momentum equation is solved to obtain the flow rate in that branch. At each fluid node, the conservation of mass, 

energy, and species equations are solved to obtain the pressures, temperatures, and species concentrations at that node. At each solid node, the energy conservation equation is solved to calculate 

temperature of the solid. 

This Technical Memorandum (TM) documents the mathematical formulation, computer 

program, and graphical user interface (GUI). Use of the code is illustrated by 12 example problems. It also documents the verification and validation effort conducted by code developers and 

users. This section also presents an overview of the subsequent sections to provide users with a 

global perspective of the code.

1.1 Network Flow Analysis Methods

The oldest method for systematically solving a problem consisting of steady flow in a pipe 

network is the Hardy Cross method.2 Not only is this method suited for hand calculations, 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 to solve large networks.3 The Newton-Raphson method solves all the 

governing equations simultaneously and is numerically more stable and reliable than successive 

substitution methods.

The network analysis method4 has been widely used in thermal analysis codes (SINDA/G5

and SINDA/FLUINT6) using an electric analog. The partial differential equation of heat conduction is discretized into finite difference form expressing temperature of a node in terms of temperatures of neighboring nodes and ambient nodes. The set of finite difference equations is solved to 

calculate temperature of the solid nodes and heat fluxes between the nodes. There have been some 

limited applications of thermal network analysis methods to model fluid flows. Such attempts 

did not go far because of the inability of heat conduction equations to handle the nonlinear fluid 

inertia term. There has been limited success in modeling compressible and two-phase flows by such 

methods. 
3

At MSFC, another system analysis code, ROCETS, is routinely used for simulating flow in 

rocket engines.7 ROCETS has a very flexible architecture where users develop the system model by 

integrating component modules such as pumps, turbines, and valves. The user can also build any 

model of specific components to integrate into the system model. ROCETS solves the system of 

equation by a modified Newton-Raphson method.8

The finite volume method (FVM) has been widely used in solving Navier-Stokes equations in CFD.9 The FVM divides the flow domain into a discrete number of control volumes and 

determines the conservation equations for mass, momentum, energy, and species for each control 

volume. Simultaneous solutions of these conservation equations provide the pressure, velocity 

components, temperature, and concentrations representative of the discrete control volumes. The 

numerical method is called ‘pressure based’ if the pressures are calculated from the mass conservation equation and density from the equation of state. On the other hand, a ‘density based’ numerical method uses mass conservation equation to calculate density of the fluid and pressure from the 

equation of state. GFSSP uses a pressure-based FVM as the foundation of its numerical scheme. 

1.1.1 Network Definitions

GFSSP constructs a fluid network using fluid and solid nodes. The fluid circuit is constructed with boundary nodes, internal nodes, and branches (fig. 1) while the solid circuit is constructed with solid nodes, ambient nodes, and conductors. The solid and fluid nodes are connected 

with solid-fluid conductors. Users must specify conditions such as pressure, temperature, and 

concentration of species at the boundary nodes. These variables are calculated at the internal nodes 

by solving conservation equations of mass, energy, and species in conjunction with the thermodynamic equation of state. Each internal node is a control volume where there is inflow and outflow 

of mass, energy, and species at the boundaries of the control volume. The internal node also has 

resident mass, energy, and concentration. The momentum conservation equation is expressed in 

flow rates and is solved in branches. At the solid node, the energy conservation equation for solid 

is solved to compute temperature of the solid node. Figure 1 shows a schematic and GFSSP flow 

circuit of a counterflow heat exchanger. Hot nitrogen gas is flowing through a pipe, colder nitrogen 

is flowing counter to the hot stream in the annulus pipe, and heat transfer occurs through metal 

tubes. The problem considered is to calculate flow rates and temperature distributions in both 

streams.

1.2 Units and Sign Conventions

GFSSP uses British gravitational units (commonly known as engineering units). Table 1 

describes the units of variables used in the code. The units in the second column are the units 

that appear in the input and output data files. Users must specify the values in these units in their 

model. The units that are listed in the third column are internal to the code and used during the 

solution of the equations. GFSSP uses standard sign conventions for mass and heat transfer. Mass 

and heat input to a node is considered positive. Similarly, mass and heat output from a node is considered negative.
4

2

12 23

89

34

910

45 56 67

28 39 410 511 612

813

1323

1314

1422

1415

1521

1516

1620 1719

1617

2324 2223 2122 2021 1920 1819

914 1015 1116 1217

1011 1112

HTC-R

FS

1

Internal Node Boundary Node

Solid-Fluid 

Conductor

Solid 

Node

3 4 5 6 7

24 23 22 21 20 19 18

8 9 10 11 12

13 14 15 16 17

HTC-R

FS

HTC-R

FS

HTC-R

FS

HTC-R

FS

HTC-R

FS

K

SS

K

SS

K

SS

K

SS

K

SS

K

SS

K

SS

K

SS

K

SS

K

SS

K

SS

K

SS

K

SS

HTC-R

FS

HTC-R

FS

HTC-R

FS

HTC-R

FS

Solid-Solid

Conductor

Douter= 4 in

Nitrogen 70 °F

Nitrogen 250 °F

L = 2 ft

 = 2.63 lb/s

Dinner= 2 in

m m = 2.59 lb/s

Figure 1. A typical flow network consists of fluid node, solid node, flow 

branches, and conductors.

Table 1. Units of variables in input/output and solver module.

Variables Input/Output Solver Module

Length

Area

Pressure

Temperature

Mass injection

Heat source

inches

inches2

psia

°F

lbm/s

Btu/s or Btu/lbm

feet

feet2

psf

R

lbm/s

Btu/s or Btu/lbm

1.3 Data Structure

GFSSP has a unique data structure (fig. 2) that allows constructing all possible arrangements of a flow network with no limit on the number of elements. The elements of a flow network 

are boundary nodes, internal nodes, and branches. For conjugate heat transfer problems, there are 

three additional elements—solid node, ambient node, and conductor. The relationship between 

a fluid node and a branch as well as a solid node and conductor is defined by a set of relational 

geometric properties. For example, the relational geometric properties of a node are number and 

name of branches connected to it. With the help of these properties, it is possible to define any 
5

Boundary Node Internal Node Branch Solid Node Ambient Node Conductor

Solid to Solid 

Conduction

Solid to Solid 

Radiation

Solid 

to Fluid

Solid 

to Ambient

Network

Fluid Solid

Figure 2. Data structure of the fluid-solid network has six major elements.

structure of the network as it progresses through every junction of the network. The positive or 

negative flow direction is also defined locally. Unlike the structured coordinate system, there is no 

global definition of flow direction and origin. The development of a flow network can start from 

any point and can proceed in any direction. 

All elements of a network have properties. The properties can be classified into two 

categories—geometric and thermofluid. Geometric properties are again classified into two subcategories—relational and quantitative. Relational properties define the relationship of the element 

with the neighboring elements. Quantitative properties include geometric parameters such as area, 

length, and volume.

1.4 Mathematical Formulation

GFSSP solves the conservation equations of mass and momentum in internal nodes and 

branches to calculate fluid properties. It also solves for energy conservation equations to calculate 

temperatures of solid nodes. Table 2 shows the mathematical closure that describes the unknown 

variables and the available equations to solve the variables. Pressure, temperature, specie concentration, and resident mass in a control volume are calculated at the internal nodes, whereas the flow 

rate is calculated at the branch. The equations are coupled and nonlinear; therefore, they are solved 

by an iterative numerical scheme. GFSSP employs a unique numerical scheme known as simultaneous adjustment with successive substitution (SASS), which is a combination of Newton-Raphson 

and successive substitution methods. The coupling of equations is shown in figure 3. The mass and 

momentum conservation equations and the equation of state are solved by the Newton-Raphson 

method while the conservation of energy and species are solved by the successive substitution 

method.
6

Table 2. Mathematical closure.

Unknown Variables Available Equations to Solve

Pressure

Flow rate

Fluid temperature

Solid temperature

Specie concentrations

Fluid mass (unsteady flow)

Mass conservation equation

Momentum conservation equation

Energy conservation equation of fluid

Energy conservation equation of solid 

Conservation equations for species

Thermodynamic equation of state

p = Pressure

 = Flow Rate

h = Enthalpy

c = Concentration

ρ = Density

m

m

Error

Mass Momentum

State Energy

Specie

Iteration Cycle

,ρ

m,c

m,h

p, m,ρ

h

p

ρ

c

m

Figure 3. Schematic of mathematical closure of GFSSP.

The total number of equations to be solved is determined from the number of internal 

nodes and branches. Figure 4 shows a typical interpropellant flow circuit in a rocket engine turbopump. In this circuit there are five boundary nodes and seven internal nodes. These nodes are 

connected by 12 branches. There are three inlet boundary nodes (48, 66, and 22) where oxygen, 

helium, and hydrogen enter into the fluid circuit. Mixtures of helium-hydrogen and helium-oxygen 

exit the circuit through boundary nodes numbered 50 and 16, respectively. At each internal node, 

four equations are solved to calculate pressure, temperature, and two concentrations. It should be 

noted that in a mixture of three components, concentrations of two components are solved. The 

concentration of the third component is determined from the fact that the sum of all concentrations must be unity. Flow rates are calculated in 12 branches. Therefore, GFSSP solves for 40 

(=7×4+12) equations to calculate all required variables in the circuit. Since the example problem 

is at steady state, resident mass in the control volume was not calculated. The transient model of 

the same fluid circuit would require solution of 47 (=7×5+12) equations at each time step of the 

simulation. The mathematical formulation has been described in detail in section 2.
7

Notes:

• Number of internal nodes = 7

• Number of branches = 12

• Total number of equations = 7×4+12 =40

• Number of equations solved by Newton

 Raphson method = 7+12 =19

• Number of equations solved by successive

 substitution method = 3×7 =21

Helium 

151 psia 

70 °F

Atmosphere 

14.7 psia 

Hydrogen

172 psia

–174 °F 

Oxygen

550 psia

–60 °F 

Branch

Boundary Node

Internal Node

Assumed Branch Flow Direction

Atmosphere 

14.7 psia 

50 xx

xx

66

48 22

16

68

49

xx

67

23

63

46

47

58 25

60 23

142

138 129

88

25

87 86

137

Figure 4. Interpropellant seal flow circuit in a rocket engine turbopump.

1.5 Fluid Properties

GFSSP is linked with two thermodynamic property programs, GASP10 and WASP11 and 

GASPAK,12 that provide thermodynamic and thermophysical properties of selected fluids. Both 

programs cover a range of pressure and temperature that allow fluid properties to be evaluated 

for liquid, liquid-vapor (saturation), and vapor region. GASP and WASP provide properties of 12 

fluids (table 3). GASPAK includes a library of 35 fluids (table 4).
8

Table 3. Fluids available in GASP and WASP.

Index Fluid

1

2

3

4

5

6

7

8

9

10

11

12

Helium

Methane

Neon

Nitrogen

Carbon Monoxide

Oxygen

Argon

Carbon Dioxide

Fluoride

Hydrogen

Water

RP-1

Table 4. Fluids available in GASPAK.

Index Fluid Index Fluid

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

Helium

Methane

Neon

Nitrogen

Carbon Monoxide

Oxygen

Argon

Carbon Dioxide

Parahydrogen

Hydrogen

Water

R-1

ISO Utane

Butane

Deuterium

Ethane

Ethylene

Hydrogen Sulfide

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

Krypton

Propane

Xenon

R-1

R1

R2

R3

R12

R12

R12

R134

R152

Nitrogen

Ammonia

Ideal

Hydrogen Peroxide

Air

1.6 Flow Resistances

In network flow analysis code, flow resistances are modeled by empirical laws. These empirical laws have been incorporated to model flow resistances for pipe flow, orifices, valves, and various pipe fittings. GFSSP models these flow resistances in the momentum conservation equation 

as friction term. There are 20 different resistance options available to users to choose from. There 

is also a provision for introducing new resistance options through user subroutines. The available 

resistance options are shown in table 5. 
9

Table 5. Resistance options in GFSSP.

Option Type of Resistance

Input 

Parameters Option

Type 

of Resistance

Input 

Parameters

1 Pipe flow L (in), D (in), ε/D 12 Flow between parallel plates r

i

 (in), c (in), L (in)

2 Flow through restriction CL, A (in2) 13 Common fittings and valves 

(two K method)

D (in), K1, K2

3 Non-circular duct a (in), b (in) 14 Pump characteristics* A0, B0, A (in2)

4 Pipe with entrance and exit loss L (in), D (in), ε/D, Ki

, Ke 15 Pump power P (hp), η, A (in2)

5 Thin, sharp orifice D1 (in), D2 (in) 16 Valve with given Cv Cv, A

6 Thick orifice L (in), D1 (in), D2 (in) 17 Joule-Thompson device LΩ, Vf

, kv, A

7 Square reduction D1 (in), D2 (in) 18 Control valve See example 12 data file

8 Square expansion D1 (in), D2 (in) 19 User defined A (in2)

9 Rotating annular duct L (in), r

o (in), ri

 (in), N (rpm) 20 Heat exchanger core Af (in2), As (in2), Ac (in2), 

L (in), Kc, Ke

10 Rotating radial duct L (in), D (in), N (rpm) 21 Parallel tube L (in), D (in), ε/D, n

11 Labyrinth seal r

i

 (in), c (in), m (in), n, α 22 Compressible orifice CL, A (in2)

* Pump characteristics are expressed as Δp = A0 + B0m 2, Δp = pressure rise, lbf/ft2, and m = flow rate, lbm/s.

1.7 Program Structure

GFSSP has three major parts (fig. 5). The first part is the GUI visual thermofluid dynamic 

analyzer for systems and components (VTASC). VTASC allows users to create a flow circuit by 

a ‘point and click’ paradigm. It creates the GFSSP input file after the completion of the model 

building process. It can also create a customized GFSSP executable by compiling and linking user 

subroutines with the solver module of the code. Users can run GFSSP from VTASC and post 

process the results in the same environment. The second major part of the program is the solver 

and property module. This is the heart of the program that reads the input data file and generates 

the required conservation equations for all internal nodes and branches with the help of thermodynamic property data. It also interfaces with user subroutines to receive any specific inputs from 

users. Finally, output files are created for VTASC to read and display results. The user subroutine 

is the third major part of the program. This consists of several blank subroutines that are called by 

the solver module. These subroutines allow the users to incorporate any new physical model, resistance option, fluid, etc. in the model. The computer program is discussed in detail in section 4.
10

• Creates Flow Circuit

• Runs GFSSP

• Displays Results Graphically

• Time-Dependent Process

• Nonlinear Boundary

• External Source Term

• Customized Output

• New Resistance/Fluid Option

• Equation Generator

• Equation Solver

• Fluid Property Program

Output Data File

Input Data 

File

Graphical User 

Interface 

(VTASC)

Solver and 

Property Module User 

Subroutines

New 

Physics

Figure 5. GFSSP’s program structure showing the interaction of three major modules.

1.8 Graphical User Interface

GFSSP’s GUI (fig. 6) provides the users a platform to build and run their models. It also 

allows post-processing of results. The network flow circuit is first built using three basic elements—

boundary node, internal node, and branch. Then the properties of the individual elements are 

assigned. Users are also required to define global options of the model that includes input/output 

files, fluid specification, and any special options such as rotation, heat exchanger, etc. During execution of the program, a run manager window opens up and users can monitor the progress of the 

numerical solution. On the completion of the run, it allows users to visualize the results in tabular 

form for steady state solutions and in graphical form for unsteady solutions. It also provides an 

interface to activate and import data to the plotting program, Winplot13 for post-processing. The 

GUI is discussed in detail in section 5. 
11

Figure 6. GFSSP’s graphical user interface, VTASC allows creating, running, 

and viewing results in one environment.

1.9 Example Problems

Several example problems have been included to aid users in becoming familiar with different options of the code. The example problems also provide the verification and validation of the 

code by comparing the code’s predictions with analytical solution and experimental data. These 

examples include: (1) Simulation of a flow system containing a pump, valve, and pipeline, (2) flow 

network for a water distribution system, (3) compressible flow in a converging-diverging nozzle, 

(4) mixing of combustion gases and a cold gas stream, (5) flow in a counterflow heat exchanger, 

(6) radial flow in a rotating radial disk, (7) flow in a squeeze film damper, (8) blowdown of a pressurized tank, (9) a reciprocating-piston cylinder, (10) power balancing of a turbopump assembly, 

(11) steady state conduction through a circular rod, and (12) fluid transient (waterhammer) due 

to sudden valve closure. These example problems are discussed in detail in section 5. 
12

2. MATHEMATICAL FORMULATION

GFSSP assumes a Newtonian, nonreacting and one-dimensional flow in the flow circuit. 

The flow can be steady or unsteady, laminar or turbulent, incompressible or compressible, with or 

without heat transfer, phase change, mixing, and rotation. The analysis of thermofluid dynamics in 

a complex network requires resolution of the system into fluid nodes and branches, and solid nodes 

and conductors. GFSSP calculates scalar properties such as pressure, temperature, and density at 

the nodes, and vector properties such as flow rates, heat fluxes at fluid branches, and conductors, 

respectively. Fluid nodes can be either internal nodes where properties are calculated or boundary 

nodes where properties are specified. Temperatures are calculated at the solid nodes and specified 

at the ambient nodes. This section describes all governing equations and the solution procedure. 

2.1 Governing Equations

Figure 7 displays a schematic showing adjacent nodes, their connecting branches, and the 

indexing system. 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.

Node 

j =2

Node 

j =1

Node 

j =3

Node 

j =4

Node 

i

 mij

 mij

 mji

 mji

 ij m = – mji

Single 

Fluid 

k=1

Single 

Fluid 

k=2

Fluid 

Mixture

Fluid 

Mixture

Figure 7. Schematic of GFSSP nodes, branches, and indexing practice.
13

2.1.1 Mass Conservation Equation 

Following is the mass conservation equation:

mτ +Δτ − mτ

Δτ = − mij

j=1

j=n

∑ . (1)

Equation (1) requires that for the unsteady formulation, the net mass flow from a given node 

must equate to the rate of change of mass in the control volume. In the steady state formulation, 

the left 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. Each term in equation (1) has the unit of lb/s.

2.1.2 Momentum Conservation Equation

The flow rate in a branch is calculated from the momentum conservation equation (eq. (2)), 

which represents the balance of fluid forces acting on a given branch. A typical branch configuration is shown in figure 8. 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 or other external momentum source. 

(mu)τ +Δτ − (mu)

gcΔτ

+ MAX m

ij

, 0 uij − u ( u ) − MAX −m

ij

, 0 uij − u ( u ) -----Unsteady----- ------------------Longitudinal Inertia--------------------

+MAX m!trans,0 uij ! u ( p ) ! MAX !m!trans,0 uij ! u ( p ) ------------------------Transverse Inertia------------------------

"

#

$

%

$

1

gc

= i p − j ( p )Aij + ρgV cosθ

gc

− Kf mij mij Aij + ρ rot

2 K ω2

A

gc

+ μ

up − uij

gcδij, p

As

--Pressure-- --Gravity-- ----Friction---- -Centrifugal- ----Shear----

−ρAnormunormuij / gc + μd

ud − uij

δij,d

− μu

uij − uu

δij,u

⎛

⎝

⎜

⎞

⎠

⎟

Aij

gc

+ S .

-Moving Boundary- ----------Normal Stress--------- -Source-

(2)
14

Branch

Node

Node

Axis of Rotation

g

ω r

j ri

j

i

 mij

Figure 8. Schematic of a branch showing gravity and rotation.

The momentum equation consists of 11 terms. There will be no occasion when all 11 terms 

will be present in a control volume. Users have the ability to include or exclude all terms except the 

pressure term. For example, the friction and shear term will never be active at the same branch. 

The shear term will be activated for multidimensional flow modeling when the friction term must 

be set to zero. The pressure term will be active under all circumstances. The left-hand side of the 

momentum equation represents the inertia of the fluid. The surface and body forces applied in the 

control volume are assembled in the right-hand side of the equation. Each term of equation (2) has 

the unit of lbf. Following are descriptions of the 11 terms:

(1) Unsteady—This term represents rate of change of momentum with time. For steady 

state flow, time step is set to an arbitrary large value and this term is reduced to zero.

(2) Longitudinal inertia—This term is important when there is a significant change in 

velocity in the longitudinal direction due to change in area and density. An upwind differencing 

scheme is used to compute the velocity differential. Flow in a converging-diverging nozzle is an 

example where this term must be active.

(3) Transverse inertia—This term is important for multidimensional flow. It accounts 

for any longitudinal momentum being transported by a transverse velocity component. Once again, 

an upwind differencing scheme is used to compute the velocity differential.

(4) Pressure—This term represents the pressure gradient in the branch. The pressures 

are located at the upstream and downstream face of a branch.
15

(5) Gravity—This term represents the effect of gravity. The gravity vector makes an angle 

(q) with the assumed flow direction vector. At q=180°, fluid is flowing against gravity; at q=90°, 

fluid is flowing horizontally and gravity has no effect on the flow. 

(6) Friction—This term represents the frictional effect. Friction was modeled as a product 

of Kf and the square of the flow rate and area. Kf 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 Kf for different 

types of flow passages is described in section 2.1.7.

(7) Centrifugal—This 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 8. Krot is the 

factor representing the fluid rotation. Krot is unity when the fluid and the surrounding solid surface 

rotate with the same speed. This term also requires knowledge of the distances from the axis of 

rotation between the upstream and downstream faces of the branch.

(8) Shear—This term represents shear force exerted on the control volume by a neighboring branch. This term is active only for multidimensional flow. The friction term is deactivated 

when this term is present. This term requires knowledge of distances between branches to compute 

the shear stress. 

(9) Moving boundary—This term represents force exerted on the control volume by 

a moving boundary. This term is not active for multidimensional calculations. 

(10) Normal stress—This term represents normal viscous force. This term is important for 

highly viscous flows. 

(11) Source—This term represents a generic source term. Any additional force acting on 

the control volume can be modeled through the source term. In a system level model, a pump 

can be modeled by this term. A detailed description of modeling a pump by this source term, S, 

appears in sections 2.1.7.14 and 2.1.7.15. 

A simplified form of the momentum equation has also been provided to compute choked 

flow rate for compressible flow in an orifice.

m!ij = CLij

A pi!igc

2"

" #1 p( ratio )

2" 1# p( ratio )

" #1

"

$

%

&

&

'

(

)

)

, (3a)

where

pratio = pj

pi

(3b)

if
16

pj

pi

< pc! ,

where

pc! = 2

! + 1

"

#

$

%

&

'

!

! (1 (3c)

and

pj

pi

= pc! . (3d)

2.1.3 Energy Conservation Equations

GFSSP solves for energy conservation equations for both fluid and solid at internal fluid 

nodes and solid nodes. Energy conservation equation for fluid is solved for all real fluids with or 

without heat transfer. For conjugate heat transfer, the energy conservation equation for solid node 

is solved in conjunction with the energy equation of fluid node. The heat transfer between solid and 

fluid node is calculated at the interface and used in both equations as source and sink terms.

2.1.3.1 Energy Conservation Equation of Fluid. The energy conservation equation for node 

i, shown in figure 7, can be expressed following first or second law of thermodynamics. The first 

law formulation uses enthalpy as the dependent variable while second law formulation uses entropy. 

The energy conservation equation based on enthalpy is shown in equation (4a):

m h − p

ρJ

⎛

⎝

⎜

⎞

⎠

⎟

τ +Δτ

− m h − p

ρJ

⎛

⎝

⎜

⎞

⎠

⎟

τ

Δτ =

j=1

j=n

∑ MAX −mij ⎡ ,0 ⎣ ⎤

⎦ hj − MAX mij ⎡ ,0 ⎣ ⎤⎦ h { i}

+

MAX −m

ij,0 ⎡

⎣⎢ ⎤⎦⎥

m

ij

i p − j ( p ) + Kijm

ij

⎡ 2

⎣

⎢ ⎤

⎦

⎥ υij ( A) + Qi . (4a)

Equation (4a) 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 (4a), 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 
17

external heat sources (Qi) coming into the node and work done on the fluid by pressure and viscous 

forces. The MAX operator used in equation (4a) is known as an upwind differencing scheme and 

has been extensively employed in the numerical solution of Navier-Stokes equations in convective heat transfer and fluid flow applications.9 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 on 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 side of the equation. For a steady state situation, the left side of equation (4a) is zero, where as in unsteady cases, the left side of the equation 

must be evaluated. 

The energy conservation equation based on entropy is shown in equation (4b):

(ms)

τ + Δτ − (ms)τ

Δτ = MAX −mij ⎡ ,0 ⎣ ⎤

⎦ sj − MAX mij ⎡ ,0 ⎣ ⎤⎦ s { i} j = 1

j = n

∑

+

MAX −mij ⎡ ,0 ⎣ ⎤

⎦

mij

⎧

⎨

⎪

⎩

⎪

⎫

⎬

⎪

⎭

⎪

S

ij, gen +

Qi

T j = 1 i

j = n

∑ .

(4b)

The entropy generation rate due to fluid friction in a branch is expressed as 

S

ij,gen =

mijΔ pij,viscous

ρuTuJ =

K f m ( ij )

3

ρuTuJ . (4c)

Equation (4b) shows that for unsteady flow, the rate of increase of entropy in the control 

volume is equal to the rate of entropy transport into the control volume plus the rate of entropy 

generation in all upstream branches due to fluid friction plus the rate of entropy added to the control volume due to heat transfer. The first term on the right-hand side of the equation represents 

the convective transport of entropy from neighboring nodes. The second term represents the rate 

of entropy generation in branches connected to the ith node. The third term represents entropy 

change due to heat transfer. Each term in equation (4b) has the unit of Btu/R-s.

2.1.3.2 Energy Conservation Equation of Solid. Typically, a solid node can be connected 

with other solid nodes, fluid nodes, and ambient nodes. Figure 9 shows a typical arrangement 

where a solid node is connected with other solid nodes, fluid nodes, and ambient nodes. The energy 

conservation equation for solid node i can be expressed as:

∂

∂τ

mCpTs

i ( ) = q

ss

j

s =1

nss

∑ + q

sf

j f =1

nsf

∑ + q

sa

j

a =1

nsa

∑ + S

i . (5)
18

Ambient Node

j

a = 2

i

j

a = 1

j

s j = 3 s = 2

j

s j = 4 s = 1

j

f j = 4 f j = 3 f j = 2 f = 1

Solid Node

Fluid Node

Figure 9. A schematic showing the connection of a solid node with neighboring solid, 

fluid, and ambient nodes.

The left-hand side of the equation represents rate of change of temperature of the solid node, i. 

The right-hand side of the equation represents the heat transfer from the neighboring node and 

heat source or sink. The heat transfer from neighboring solid, fluid and ambient nodes can be 

expressed as

q

ss = kijs

Aijs / δijs Ts

j

s −Ts

i ( ) , (5a)

q

sf = hij f

Aij fTf

j f −Ts

⎛ i

⎝

⎞

⎠ , (5b)

 

q

sa = hija

Aija Ta

j

a −Ts

i ( ) . (5c)

The heat transfer rate can be expressed as a product of conductance and temperature differential. The conductance for equations (5a)–(5c) are as follows:

Cijs =

kijs

Aijs

δijs

, Cij f = hij fAij f, and Cij f = hijaAija , (5d)

respectively, where effective heat transfer coefficients for solid to fluid and solid to ambient nodes 

are expressed as:
19

hij f = hc,ij f+ hr,ij f

hija = hc,ija + hr,ija (5e)

and

hr,ij f =

σ Tf

j ⎛ f

⎝

⎞

⎠

2

+ Ts

i ( ) ⎡ 2

⎣

⎢ ⎤

⎦

⎥ Tf

j f +Ts

⎡ i

⎣

⎢ ⎤

⎦

⎥

1 / εij, f +1/ εij,s −1

hr,ija =

σ Ta

j

a ( )

2

+ Ts

i ( ) ⎡ 2

⎣

⎢ ⎤

⎦

⎥ Ta

j

a +Ts⎡ i

⎣ ⎤⎦

1/ εij,a +1/ εij,s −1 . (5f)

GFSSP provides users with three different options for specifying heat transfer coefficient:

(1) User can provide a constant heat transfer coefficient.

(2) User can specify Dittus-Boelter equation14 for single phase flow where the Nusselt 

number is expressed as:

hcD

kf

= 0.023(Re)

0.8 (Pr)0.33 , (5g)

where Re=ruD/mf and Pr=Cp mf /kf.

(3) User can specify modified Miropolosky’s correlation15 for two-phase flow:

Nu = 0.023 Remix ( )

0.8

Pr

v ( )

0.4 (Y )

Remix = ρuD

μv

⎛

⎝

⎜

⎞

⎠

⎟ x +

ρv

ρl

⎛

⎝

⎜

⎞

⎠

⎟ (1− x) ⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

Pr

v = Cpμv

kv

⎛

⎝

⎜

⎞

⎠

⎟

Y = 1− 0.1 ρv

ρl

− 1 ⎛

⎝

⎜

⎞

⎠

⎟

0.4

(1− x)

0.4 . (5h)
20

Equation (5h) can be rearranged to determine Ts

i

:

Ts

i =

Cijs

Ts

j

s + Cij f

Tf

j f + Cija

Ta

j

a

j

a =1

nsa

∑ +

(mCp )m

Δτ

Ts,m

i + S

j f =1

nsf

∑

j

s =1

nss

∑

mCp

Δτ

+ Cijs + Cij f+ Cija j

a =1

nsa

∑

j f =1

nsf

∑

j

s =1

nss

∑

. (6)

2.1.4 Fluid Specie Conservation Equation 

For a fluid mixture, density is a function of mass fraction of fluid species. 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 kth specie can be written as 

(mici,k )τ + Δτ − (mici, k )τ

Δτ = MAX −m ij

,0 ⎡

⎣⎢ ⎤⎦⎥ c j, k − MAX m ij

,0 ⎡

⎣⎢ ⎤⎦⎥ c { i,k} j = 1

j = n

∑ + 

Si,k . (7)

For a transient flow, equation (7) states that the rate of increase of the concentration of kth 

specie in the control volume equals the rate of transport of the kth specie into the control volume 

minus the rate of transport of the kth specie out of the control volume plus the generation rate of 

the kth specie in the control volume.

Like equation (4a), for steady state conditions, equation (7) requires that the net mass flow 

of the kth 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 side of equation (7) is zero. For the unsteady formulation, the resident 

mass in the control volume is changing and therefore the left side must be computed. Each term 

in equation (5) has the unit of lb/s.

2.1.5 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. GASP10

provides the thermodynamic and transport properties for 10 fluids. These fluids include hydrogen, 

oxygen, helium, nitrogen, methane, carbon dioxide, carbon monoxide, argon, neon, and fluorine. 

WASP11 provides the thermodynamic and transport properties for water and steam. For RP-1 fuel, 

a lookup 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. GASPAK12

provides thermodynamic properties for helium, methane, neon, nitrogen, carbon monoxide, 
21

oxygen, argon, carbon dioxide, hydrogen, parahydrogen, water, RP-1, 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, ammonia, hydrogen peroxide, and air. 

2.1.5.1 Equation of State for a Real Fluid. Transient flow calculations require the knowledge of resident mass in a control volume. The resident mass is calculated from the equation of 

state for real fluid that can be expressed as

m = pV

zRT . (8)

It may be noted that equation (8) is valid for liquid, gas, and gas-liquid mixture. For an ideal 

gas compressibility factor, z is unity. The compressibility factor for real gas is computed from the 

equation of state of real fluids using the above-mentioned thermodynamic property programs. For 

two-phase mixture, z is computed from the following relation:

z = p

ρmixRT , (8a)

where

ρmix = ρl ρg

xρl − (1− x) ρg

(8b)

and

x = s − sl

sg − s

l

. (8c)

2.1.6 Mixture Property Calculations

This section describes a procedure developed for GFSSP to estimate the density and temperature of mixtures of real fluids. Let us assume that n fluids are mixing in the ith node. At node i, 

pressure (pi) and molar concentrations (xk) are known. The problem is to calculate the density (ρi), 

temperature (Ti), specific heat (Cp), specific heat ratio (g), and viscosity (μ) of the mixture at the ith 

node.

During iterative calculations, GFSSP calculates the mixture properties using the following 

steps:

(1) Calculate ρk, μk, Cp,k and γk from pi and Ti using the thermodynamic property routines 

of the program.

(2) Calculate the compressibility factor of each component of the mixture (zk) from the 

equation of state for a real gas (eq. (9)):

zk = i p

k ρ RkTk

, (9)

where Rk is the gas constant for kth fluid.
22

(3) Calculate compressibility factor of the mixture (zi) as shown in equation (10) by taking 

the molar average of the component compressibility factors obtained in step (2):

.

zi = xkzk

k=1

k=n

∑ . (10)

(4) Calculate gas constant of the mixture (Ri) by taking the molar average of the component gas constants:

Ri = xkRk

k=1

k=n

∑ . (11)

(5) Calculate the density of the mixture (ρi) from the equation of state of a real fluid:

ρi = i p

zi

RiTi

. (12)

(6) Calculate viscosity, specific heat, and specific heat ratio of the mixture by taking the 

molar average of the component properties μk, Cp,k, and γ k as shown in equations (13)–(15):

i μ = xkμk

k=1

k=n

∑ , (13)

Cp,i = Cp,k xkMk

x k=1 kMk

k=n

∑ , (14)

and

γ i xkγ k

k=1

k=n

∑ . (15)

(7) Calculate Ti from the energy conservation equation expressed in terms of a product 

of specific heat and temperature instead of enthalpy:

τ +Δτ

T( i) =

Cp,k, j xk, jTj

MAX −mij ⎡ ,0 ⎣ ⎤

⎦

k=1

k=nf

∑

j=1

j=n

∑ + Cp,i mi

T ( i)τ/ Δτ +Qi

Cp,k, j xk, j MAX mij ⎡ ,0 ⎣ ⎤

⎦ + (Cp,i m)τ

/ Δτ

k=1

k=nf

∑

j=1

j=n

∑

, (16)

where Cp,k is the molar specific heat and xk is the mole-fraction of the kth specie. Note that 

an unsteady formulation of the energy equation (eq. (4)) was used to compute Ti. 
23

2.1.7 Friction Calculations

It was mentioned earlier in this TM that the friction term in the momentum equation is 

expressed as a product of Kf , the square of the flow rate and the flow area. Empirical information 

is necessary to estimate Kf . Several options for flow passage resistance are listed in table 5. In this 

subsection, the expression of Kf for all resistance options are described.

2.1.7.1 Pipe Flow (Branch Option 1). Figure 10 shows the pipe resistance option parameters 

that are required by GFSSP. This option considers that the branch is a pipe with length (L), 

diameter (D), and surface roughness (ε). For this option, Kf can be expressed as:

Kf = 8 fL

ρuπ2D5gc

, (17)

where ru is the density of the fluid at the upstream node of a given branch. The derivation of Kf

for pipe flow is covered in appendix A.

L

D

Detail A Detail A

D = Pipe Diameter

L = Pipe Length

 = Absolute Roughness

Figure 10. Pipe resistance option parameters.

The Darcy friction factor (f) is determined from the Colebrook equation,16 which 

is expressed as:

1

f = −2log ε

3.7D

+

2.51

Re f

⎡

⎣

⎢ ⎤

⎦

⎥ , (18) 

where e /D and Re are the surface roughness factor and Reynolds number, respectively. 

2.1.7.2 Flow Through a Restriction (Branch Option 2). This option regards the branch 

as a flow restriction with a given flow coefficient (CL), and area (A). For this option, Kf can be 

expressed as:

Kf = 1

2 cg ρu L

2 C 2 A . (19)
24

In classical fluid mechanics, head loss is expressed in terms of a nondimensional ‘K factor’:

Δh = K

2u

2g . (20)

K and CL are related as:

CL = 1

K . (21)

2.1.7.3 Noncircular Duct (Branch Option 3). This option considers a duct with a noncircular cross section. Four different types of cross sections can be modeled as shown in 

figure 11. 

a

a

b

b

a

b

a b/2

(a) (b)

(c) (d)

Figure 11. Noncircular duct cross section: (a) Rectangle, (b) ellipse,

(c) concentric cylinder, and (d) circular sector.

White describes a procedure to estimate the friction factor in a noncircular duct.17 This 

procedure consists of the following steps:

(1) Estimate the hydraulic diameter of the cross section: Dh = (4)(area)/perimeter.

(2) Estimate the Poiseuille number (Po) for a particular cross section. The Poiseuille number 

can be expressed as a polynomial function of aspect ratio as shown in equation (22). Table 6 

provides the coefficients for different geometries.

Po = A0 + A1

b

a

⎛

⎝

⎜ ⎞

⎠

⎟ + A2

b

a

⎛

⎝

⎜ ⎞

⎠

⎟

2

+ A3

b

a

⎛

⎝

⎜ ⎞

⎠

⎟

3

+ A4

b

a

⎛

⎝

⎜ ⎞

⎠

⎟

4

. (22)
25

Table 6. Poiseuille number coefficients for noncircular duct cross sections.

Coefficients Rectangle Ellipse

Concentric* 

Cylinder 

Circular

Section

A0

A1

A2

A3

A4

23.9201

–29.436

30.3872

–10.7128

–

19.7669

–4.53458

–11.5239

22.3709

–10.0874

22.0513

6.44473

–7.35451

2.78999

–

11.9852

3.01553

–1.09712

–

–

* For b/a < 0.2508 P0 = A0

b

a

⎛

⎝ ⎞⎠

A1

, where, A0 = 24.8272, A1 = 0.0479888.

(3) Calculate the friction factor for a noncircular pipe:

• Laminar flow (Re<2,300)

f = 4Po

Re . (23)

• Turbulent flow

– Compute the effective diameter:

Deff = 16Dh

Po . (24)

– Compute the effective Reynolds number:

Reeff = m

μ

Deff

A . (25)

– Compute the friction factor using the Colebrook equation (eq. (18)).

(4) Compute Kf from the following expression:

Kf = 8 fL

uρ 2 π h

5 D c g . (26)

2.1.7.4 Pipe With Entrance and Exit Loss (Branch Option 4). Figure 12 shows the pipe with 

entrance and/or exit loss resistance option parameters that are required by GFSSP. This option is 

an extension of option 1. In addition to the frictional loss in the pipe, entrance and exit losses are 

also calculated. For this option, Kf can be expressed as:
26

L

D

Detail A

Entrance Exit

Detail A

D = Pipe Diameter

L = Pipe Length

 = Absolute Roughness

Ki = Entrance Loss Coefficient

Ke = Exit Loss Coefficient

Figure 12. Pipe with entrance and/or exit loss resistance option parameters.

Kf = 8Ki

ρu

2 π D4

c g

+

8 fL

ρu

2 π D5

c g

+

8Ke

ρu

2 π D4

c g , (27)

where Ki and Ke are the entrance and exit loss coefficients, respectively.

2.1.7.5 Thin, Sharp Orifice (Branch Option 5). Figure 13 shows the thin, sharp orifice resistance option parameters that are required by GFSSP. This option considers the branch as a thin, 

sharp orifice with a pipe diameter of D1 and an orifice diameter of D2. For this option, Kf can be 

expressed as:18

Kf = K1

2 c g ρu

2 A , (28)

where A = πD1

2

4 .

D1 = Pipe Diameter

D2 = Orifice Throat Diameter D2 D1

Thin, Sharp Orifice

Figure 13. Thin, sharp orifice resistance option parameters.
27

For upstream Re ≤ 2,500,

K1 = 2.72 +

D2

D1

⎛

⎝

⎜

⎞

⎠

⎟

2

120

Re − 1 ⎛

⎝

⎜ ⎞

⎠

⎟

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ 1− D2D1

⎛

⎝

⎜

⎞

⎠

⎟

2 ⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

D1

D2

⎛

⎝

⎜

⎞

⎠

⎟

4

−1

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ . (29)

For upstream Re > 2,500,

K1 = 2.72 − D2

D1

⎛

⎝

⎜

⎞

⎠

⎟

2

4,000

Re

⎛

⎝

⎜ ⎞

⎠

⎟

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ 1− D2D1

⎛

⎝

⎜

⎞

⎠

⎟

2 ⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

D1

D2

⎛

⎝

⎜

⎞

⎠

⎟

4

−1

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ . (30)

This option is recommended for subsonic and incompressible flow.

2.1.7.6 Thick Orifice (Branch Option 6). Figure 14 shows the thick orifice resistance option 

parameters that are required by GFSSP. This option models the branch as a thick orifice with 

a pipe diameter of D1, an orifice diameter of D2, and orifice length of Lor. This option should be 

used if Lor /D2 ≤ 5. If Lor /D2 > 5, the user should use a square expansion, option 8, or a square 

reduction, option 7. For option 6, Kf can be expressed as in equation (28). However, the K1 in 

equation (28) is calculated18 from the following expressions:

For upstream Re ≤ 2,500:

K1 = 2.72 +

D2

D1

⎛

⎝

⎜

⎞

⎠

⎟

2

120

Re − 1 ⎛

⎝

⎜ ⎞

⎠

⎟

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ 1− D2

D1

⎛

⎝

⎜

⎞

⎠

⎟

2 ⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

D1

D2

⎛

⎝

⎜

⎞

⎠

⎟

4

−1

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ 0.584 +

0.0936

Lor ( / D2 )

1.5

+ 0.225

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ . (31)

D1 = Pipe Diameter

D2 = Orifice Throat Diameter

Lor = Orifice Length

D2 D1

Lor

Figure 14. Thick orifice resistance option parameters.
28

For upstream Re > 2,500:

K1 = 2.72 − D2

D1

⎛

⎝

⎜

⎞

⎠

⎟

2

4,000

Re

⎛

⎝

⎜ ⎞

⎠

⎟

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ 1− D2

D1

⎛

⎝

⎜

⎞

⎠

⎟

2 ⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

D1

D2

⎛

⎝

⎜

⎞

⎠

⎟

4

−1

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ 0.584 +

0.0936

Lor ( / D2 )

1.5

+ 0.225

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ . (32)

This option is recommended for subsonic and incompressible flow.

2.1.7.7 Square Reduction (Branch Option 7). Figure 15 shows the square reduction resistance option parameters that are required by GFSSP. This option considers the branch as a square 

reduction. The diameters of the upstream and downstream pipes are D1 and D2, respectively. For 

this option, Kf can be expressed as in equation (28). However, the K1 in equation (28) is calculated 

from the following expressions:18

For upstream Re ≤ 2,500,

K1 = 1.2 +

160

Re

⎡

⎣

⎢ ⎤

⎦

⎥

D1

D2

⎛

⎝

⎜

⎞

⎠

⎟

4

− 1

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ . (33)

For upstream Re > 2,500,

K1 = [0.6 + 0.48 f ] D1

D2

⎛

⎝

⎜

⎞

⎠

⎟

2

D1

D2

⎛

⎝

⎜

⎞

⎠

⎟

2

−1

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

2

. (34)

The Reynolds number and friction factor that are utilized within these expressions are based 

on the upstream conditions. The user must specify the correct flow direction through this branch. 

If the model determines that the flow direction is in the reverse direction, the user will have to 

replace the reduction with an expansion and rerun the model.

D1 = Upstream Pipe Diameter

D2 = Downstream Pipe Diameter D D 2 1 Flow

Figure 15. Square reduction resistance option parameters.
29

2.1.7.8 Square Expansion (Branch Option 8). Figure 16 shows the square expansion resistance option parameters that are required by GFSSP. This option considers the branch as a square 

expansion. The diameters of the upstream and downstream pipes are D1 and D2, respectively. For 

this option, Kf can be expressed as in equation (28). However, the K1 in equation (28) is calculated 

from the following expressions:18

For upstream Re ≤ 4,000,

K1 = 2 1− D1

D2

⎛

⎝

⎜

⎞

⎠

⎟

4 ⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥ . (35)

For upstream Re > 4,000:

K1 = [1+ 0.8 f ] 1− D1

D2

⎛

⎝

⎜

⎞

⎠

⎟

2 ⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

2

. (36)

The Reynolds number and friction factor that are utilized within these expressions are based 

on the upstream conditions. The user must specify the correct flow direction through this branch. 

If the model determines that the flow direction is in the reverse direction, the user will have to 

replace the expansion with a reduction and rerun the model.

D1 = Upstream Pipe Diameter

D2 = Downstream Pipe Diameter D1 D2 Flow

Figure 16. Square expansion resistance option parameters.

This option is recommended for subsonic and incompressible flow.

2.1.7.9 Rotating Annular Duct (Branch Option 9). Figure 17 shows the rotating annular 

duct resistance option parameters that are required by GFSSP. This option considers the branch as 

a rotating annular duct. The length, outer radius, and inner radius of the annular passage are L, ro, 

and ri, respectively. The inner surface is rotating at N rpm (N=30 w/p). For this option, Kf can be 

expressed as: 
30

L = Duct Length (Perpendicular to Page)

b = Duct Wall Thickness (b = r

o– ri

)

 = Duct Rotational Velocity

r

i = Duct Inner Radius

r

o = Duct Outer Radius

r

i

r

o

Figure 17. Rotating annular duct resistance option parameters.

K f = fL

ρu

2 π 2 A cg r

o − r ( i) . (37)

The friction factor (f ) in equation (37) was calculated from the following expressions:19

f

0T = 0.077(Ru)

−0.24 , (38)

where

Ru = ρuu 2 r

0r ( i)

μ

(39)

and u is the mean axial velocity, therefore:

f

f0T

= 1+ 0.7656

ωr

i

2u

⎛

⎝

⎜ ⎞

⎠

⎟

2 ⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

0.38

. (40)

2.1.7.10 Rotating Radial Duct (Branch Option 10). Figure 18 shows the rotating radial duct 

resistance option parameters that are required by GFSSP. This option considers the branch as 

a rotating radial duct. This option accounts only for the frictional losses encountered with this type 

of flow. Since centrifugal effects are also important in a rotating radial duct, the user must select 

this option and activate the rotational term in the momentum conservation equation (2). 
31

L = Duct Length

 = Duct Rotational Velocity

D = Duct Diameter

L

r

j

r

i

D

Axis of 

Rotation

Figure 18. Rotating radial duct resistance option parameters.

The length and diameter of the duct are L and D, respectively. The rotational speed is 

ω rad/s. For this option, Kf can be expressed as:

Kf = 8 fL

ρuπ2D5gc

. (41)

The friction factor (f ) in equation (41) was calculated from the following equations:20

0T f = 0.0791(Re)

−0.25 (42)

and

f

0T f = 0.942 + 0.058 ωD

u

⎛

⎝

⎜ ⎞

⎠

⎟ ωD2

ν

⎛

⎝

⎜

⎞

⎠

⎟

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

0.282

. (43)

2.1.7.11 Labyrinth Seal (Branch Option 11). Figure 19 shows the labyrinth seal resistance 

option parameters that are required by GFSSP. This option considers the branch as a labyrinth 

seal. The number of teeth, clearance, and pitch are n, c, and m respectively. For this option, 

Kf can be expressed as:
32

c = Clearance

m = Gap Length (Pitch)

r

i

 = Radius (Tooth Tip)

n = Number of Teeth

r

i

c

m

Figure 19. Labyrinth seal resistance option parameters.

Kf = 1/ 2 ( ε + 0.5) n + 1.5

2 c g ρu

2 α 2 A , (44)

where the carryover factor, e, is expressed as:

ε = 1

1− (n − 1) c / mp

n(c / m + 0.02)

. (45)

For a straight labyrinth seal, α should be set to unity. For a stepped labyrinth seal, a should 

be less than unity. A value of 0.9 has been recommended for many rocket engine turbopump 

applications.

2.1.7.12 Flow Between Parallel Plates (Branch Option 12). Figure 20 shows the parallel flat 

plate resistance option parameters that are required by GFSSP. This option considers the branch as 

having laminar flow between parallel flat plates. A face seal can be modeled using this option. The 

flow is assumed to occur between two parallel plates separated by a distance equal to the clearance 

between the shaft and the housing. The effect of curvature is neglected. The length, inner diameter, 

and clearance of the seal are L, D, and c, respectively. For this option, Kf can be expressed as:21
33

c = Seal Thickness (Clearance)

B = Passage Width (B= D)

L = Seal Length

L

L

c

c

D

Flow

Flow

B

Figure 20. Parallel flat plate resistance option parameters.

K f = 12μ Lρ

πgcDc

3

m . (46)

2.1.7.13 Common Fittings and Valves (Branch Option 13). This option considers the branch 

as common fittings or valves. The resistance in common fittings and valves can be computed 

by the two-K method.22 For this option, Kf can be expressed as:

K f = K1 / Re+ K∞ (1+1/ D)

2gc ρu

2 A , (47)

where

K1 = K for the fitting at Re =1.

K ∞ = K for the fitting at Re =∞ .

D = internal diameter of attached pipe (in).

The constants K1 and K ∞ for common fittings and valves are listed in table 7.
34

Table 7. Constants for two-K method of hooper for fittings/valves (GFSSP resistance option 13). 

Fitting Type K1 K∞

90° elbows Standard (R/D = 1), screwed

Standard (R/D = 1), flanged or welded

Long radius (R/D = 1.5), all types

800

800

800

0.40

0.25

0.20

Mitered (R/D = 1.5) 1 weld (90° angle)

2 welds (45° angle)

3 welds (30° angle)

4 welds (22.5° angle)

5 welds (18° angle)

1,000

800

800

800

800

1.15

0.35

0.30

0.27

0.25

45° elbows Standard (R/D = 1), all types

Long radius (R/D = 1.5), all types

Mitered, 1 weld, 45° angle

Mitered, 2 weld, 22.5° angle

500

500

500

500

0.20

0.15

0.25

0.15

180° elbows Standard (R/D = 1), screwed

Standard (R/D = 1), flanged or welded

Long radius (R/D = 1.5), all types

1,000

1,000

1,000

0.60

0.35

0.30

Tee, used as elbow Standard, screwed

Long radius, screwed

Standard, flanged or welded

Stub-in-type branch

500

800

800

1,000

0.70

0.40

0.80

1.00

Tee, flow through Screwed

Flanged or welded

Stub-in-type branch

200

150

100

0.10

0.50

0.0

Valves Gate, ball, plug

(β = dorifice/dpipe)

Full line size, β = 1.0

Reduced trim, β = 0.9

Reduced trim, β = 0.8

300

500

1,000

0.10

0.15

0.25

Globe, standard

Globe, angle or Y-type

Diaphragm, dam type

Butterfly

1,500

1,000

1,000

800

4.0

2.0

2.0

0.25

Check Lift

Swing

Tilting Disk

2,000

1,500

1,000

10.0

1.5

0.5

2.1.7.14 Pump Characteristics (Branch Option 14). This option considers the branch 

as a pump with given characteristics. The pump characteristics must be expressed as: 

Δp = A0 + B0m 2 , (48)

where

D p = pressure rise (lbf/ft2).

m = flow rate (lbm/s).

The momentum source (S) in equation (2) is then expressed as:

S=D p A . (49)
35

2.1.7.15 Pump Horsepower (Branch Option 15). This option considers the branch 

as a pump with a given horsepower (P) and efficiency (h). The momentum source (S) in equation (2) is then expressed as:

S = 550 u ρ Pη A

m . (50)

2.1.7.16 Valve With a Given Loss Coefficient (Branch Option 16). This option considers the 

branch as a valve with a given flow coefficient (Cv). The flow coefficient is the volume (in gallons) 

of water at 60 °F that will flow per minute through a valve with a pressure drop of 1 psi across the 

valve. The recommended formula for Cv determination with water is:

Cv = Q 1

Δp , (51)

where Q is the volumetric flow rate in gallons per minute of water at 60° F and ∆	p is the pressure 

drop in psia. For this option, Kf can be expressed as: 

Kf = 4.68 ×105

ρuCv

2 . (52)

2.1.7.17 Joule-Thompson Device (Branch Option 17). This option considers the branch as 

a Viscojet,23 which is a specific type of flow resistance with relatively large flow passages with very 

high pressure drops. The flow rate through the Viscojet is given by:

w = 10,000 kv

Vf

LΩ

Δp S (1− x) , (53)

where w is the flow rate in lbm/hr, LΩ is the resistance of the fluid device, kv is an empirical factor, 

and Vf is the viscosity correction factor.

For this option, Kf can be expressed as:

Kf = 18.6624

S

LΩ

Vf kv (1− x)

⎛

⎝

⎜

⎞

⎠

⎟

2

. (54)

2.1.7.18 Control Valve (Branch Option 18). This is an exclusively transient option that considers the branch as a control valve that monitors the pressure at some arbitrary point downstream 

of the valve and opens and closes to maintain that pressure within a prescribed tolerance. This 

option was originally developed for use with the pressurization option to model on/off, or ‘bangbang,’ pressurization systems as shown in figure 21. The valve is regarded as a flow restriction with 

a given flow coefficient (CL) and area (A), and uses equation (18) to calculate Kf for the valve. 
36

P

Ullage

Propellant

Figure 21. Control valve in a pressurization system.

The remaining inputs for the control valve option define its open/close characteristics. The 

suboption determines the type of open/close response of the valve (instantaneous, linear, or nonlinear) and the valve initial position describes whether the valve is initially open or closed. The control 

node defines the location in the model where the control valve option is to monitor and maintain 

the pressure while the pressure tolerance file name provides the code with the name of the file 

containing the pressure tolerance data for the control valve. For a linear open/close response, the 

time to open/close and the number of time steps needed to complete that response are provided as 

additional inputs. Finally, for a nonlinear open/close response, the file names for the open and close 

characteristics of the valve are required as additional inputs.

2.1.7.19 User-Defined Resistance (Branch Option 19). This option allows users to create a new resistance option that is not available in GFSSP library. Once this option is chosen, the 

user is required to supply the coding for calculating Kf for this option in the user subroutine to be 

described in the following section. In the preprocessor, the user is required to supply the crosssectional area of the branch.

2.1.7.20 Heat Exchanger Core (Branch Option 20).This option considers branch as a heat 

exchanger core. In a typical heat exchanger core (fig. 22), the fluid goes past a tube bank to allow 

heat transfer between fluids in the main duct and fluids within the tubes. Free flow area is reduced
37

1 2

Figure 22. Heat exchanger core.

and there is larger surface area of contact between fluid and solid walls. Sections 1 and 2 in figure 

22 represent inlet and outlet of the heat exchanger core, respectively. The pressure drop through the 

heat exchanger core can be expressed as:24

Δp = G2

2gc

v1 Kc + 1− σ2 ( ) + 2

v2

v1

−1 ⎛

⎝

⎜

⎞

⎠

⎟ + f

As

Ac

vm

v1

− 1− σ2 ( − Ke )

v2

v1

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

. (55)

Equation (55) can be rewritten as:

Δ p = 1

2ρ1 gc Ac

2 Kc +1− σ2 ( ) + 2 ρ1

ρ2

− 1 ⎛

⎝

⎜

⎞

⎠

⎟ + f

As

Ac

ρ1

ρavg

− 1− σ2 ( − Ke ) ρ1

ρ2

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

m 2 . (56)

Therefore, Kf can be expressed as:

Kf =

Kc + 1− σ2 ( ) + 2 ρ1

ρ2

−1 ⎛

⎝

⎜

⎞

⎠

⎟ + f

As

Ac

ρ1

ρavg

− 1− σ2 ( − Ke ) ρ1

ρ2

2ρ1gcAc

2 . (57)

2.1.7.21 Parallel Tube (Branch Option 21). This option considers branch as a parallel tube 

where fluid flows through n number of tubes (fig. 23). The flow is assumed to be distributed uniformly in all tubes. This resistance option calculates the pressure drop in the parallel tube. For this 

option, Kf can be expressed as:

Kf = 8 fL

ρuπ2D5gcn

2 . (58)

2.1.7.22 Compressible Orifice (Branch Option 22). This option considers branch as an orifice for compressible flow. In this option, unlike other options, flow rate is calculated from a simplified momentum equation (eq. (3c)). There is no need to calculate Kf

 for this option. The input to 

this option is identical to option 2 (flow-through restriction).
38

Figure 23. Parallel tube.

2.2 Solution Procedure

GFSSP numerically solves the governing equations described in the previous section to 

compute pressure, temperature, flow rate, and other fluid properties in a given flow circuit. The 

mathematical closure is described in table 8 where each variable and the designated governing equation to solve that variable are listed. It may be noted that the pressure is calculated from the mass 

conservation equation although pressure does not explicitly appear in equation (1). This is, however, possible in the iterative scheme where pressures are corrected to reduce the residual error in 

mass conservation equation. This practice was first implemented in a semi-implicit pressure linked 

equation (SIMPLE) algorithm proposed by Patankar and Spalding25 and commonly referred to as 

‘pressure based’ algorithm in computational fluid dynamics literature. The momentum conservation equation (eq. (2)) which contains both pressure and flow rate is solved to calculate the flow 

rate. The strong coupling of pressure and flow rate requires that mass and momentum conservation 

equations are solved simultaneously.

Table 8. Mathematical closure.

Variable 

Number

Variable 

Name

Designated Equation 

to Solve the Variable

1

2

3

4

5

6

Pressure

Flow rate

Fluid enthalpy or entropy

Solid temperature

Specie concentration

Fluid mass

Mass conservation

Momentum conservation

Energy conservation of fluid

Energy conservation of solid

Species conservation

Thermodynamic state

The energy conservation equation can either be expressed in terms of enthalpy or entropy. 

The temperature, density, and other thermophysical properties, such as viscosity, specific heats 

are computed either from pressure and enthalpy or from pressure and entropy using thermodynamic property programs, GASP10/WASP11 or GASPAK.12 In flow circuits where solid to fluid 

heat transfer is present, the energy conservation equation for solid is solved to calculate the solid 

temperature. The rate of heat transfer between solid to fluid appears as source or sink term in the 

energy conservation equations of fluid and solid. 
39

For a mixture, the conservation equations of species are solved to compute the mass fraction of species. Unlike single fluid, the energy equation is solved in terms of temperature (eq. (15)) 

instead of enthalpy or entropy. For each specie density and other thermophysical properties such 

as viscosity, specific heats are computed from pressure and temperature using the above-mentioned 

thermodynamic property programs. 

For a transient problem, fluid mass is required in mass and momentum conservation equations (eqs. (1) and (2)). GFSSP uses the thermodynamic equation of state (eq. (8)) to calculate 

resident mass in an internal node where density, compressibility factor, and temperature are 

computed. 

There are two types of numerical methods available to solve a set of nonlinear coupled algebraic equations: (1) the successive substitution method and (2) the Newton-Raphson method. In 

the successive substitution method, each equation is expressed explicitly to calculate one variable. 

The previously calculated variable is then substituted into the other equations to calculate another 

variable. In one iterative cycle, each equation is visited. The iterative cycle is continued until the 

difference in the values of the variables in successive iterations becomes negligible. The advantages 

of the successive substitution method are its simplicity to program and its low code overhead. The 

main limitation, however, is finding an optimum order for visiting each equation in the model. 

This visiting order, which is called the information flow diagram, is crucial for convergence. Underrelaxation (partial substitution) of variables is often required to obtain numerical stability. Details 

of the successive substitution method appear in appendix B.

In the Newton-Raphson method, the simultaneous solution of a set of nonlinear equations 

is achieved through an iterative guess and correction procedure. Instead of solving for the variables 

directly, correction equations are constructed for all of the variables. The intent of the correction 

equations is to eliminate the error in each equation. The correction equations are constructed in 

two steps: (1) The residual errors in all of the equations are estimated and (2) the partial derivatives of all of the equations, with respect to each variable, are calculated. The correction equations 

are then solved by the Gaussian elimination method. These corrections are then applied to each 

variable, which completes one iteration cycle. These iterative cycles of calculations are repeated 

until the residual error in all of the equations is reduced to a specified limit. The Newton-Raphson 

method does not require an information flow diagram. Therefore, it has improved convergence 

characteristics. The main limitation to the Newton-Raphson method is its requirement for a large 

amount of computer memory. Details of the Newton-Raphson method appear in appendix C. 

In GFSSP, a combination of the successive substitution method and the Newton-Raphson 

method is used to solve the set of equations. This method is called SASS (simultaneous adjustment 

with successive substitution) (fig. 24). In this scheme, the mass and momentum conservation equations are solved by the Newton-Raphson method. The energy and specie conservation equations 

are solved by the successive substitution method. The underlying principle for making such a division was that the equations that are more strongly coupled are solved by the Newton-Raphson 

method. The equations that are not strongly coupled with the other set of equations are solved by 

the successive substitution method. Thus, the computer memory requirement can be significantly 

reduced while maintaining superior numerical convergence characteristics.
40

Iteration Loop Governing Equations Variables

Mass Conservation Pressure

Simultaneous 

Solution

Successive 

Substitution

Property 

Calculation

Convergence 

Check

Flow Rate

Resident Mass

Enthalpy

Solid Temperature

Temperature, Density, 

Compressibility, Factor 

Viscosity, etc.

Momentum Conservation

Equation of State

Energy Conservation of Fluid

Energy Conservation of Solid

Thermodynamic Property 

Program

Figure 24. SASS (simultaneous adjustment with successive substitution) scheme 

for solving governing equations.

SASS has two options available. In one option, there are two iterative loops—inner and 

outer. In the inner iterative loop, mass and momentum conservation equations are solved by the 

Newton-Raphson scheme. For unsteady formulation, the equation of state is also solved by the 

Newton-Raphson scheme in addition to mass and momentum conservation equations. In the 

outer loop, the energy and specie conservation equations are solved by the successive substitution 

method. The outer loop also calculates the density and other thermodynamic and thermophysical properties and the flow resistance coefficient (Kf

) which is a function of density. This option is 

called nonsimultaneous option. The total number of iterations in this option can be expressed as:

Ntotal = ni

i=1

n0

∑ , (59)

where n0 is the number of outer iterations and ni is the number of inner iterations. The inner iterative cycle is terminated when the normalized maximum correction (D max), is less than the convergence criterion (Cc). Dmax is determined from 

Δmax = MAX

Φ′

i

i=1 Φi

NE

∑ , (60)

where NE is the total number of equations solved by the Newton-Raphson scheme 
41

NE = Number of nodes + number of branches (steady flow)

NE = Number of nodes × 2 + number of branches (unsteady flow).

The outer iteration is terminated when ∆o

max is less than the convergence criterion (Cc). 

Δo

max is determined from 

Δ max ° = MAX ΔKf , Δρ,Δh, or Δs , (61)

where

ΔKf = MAX

K′

f

i=1 Kf

NB

∑

Δρ = MAX ρ′

ρ i=1

NN

∑

Δh = MAX h′

h i=1

NN

∑

 Δs = MAX s′

s i=1

NN

∑ . (62)

NB and NN are the number of branches and nodes, respectively, in a flow circuit.

In the second option, there is only one iterative loop. During the iterative cycle mass, 

momentum and the equation of state are first solved by the Newton-Raphson scheme. Then 

the energy and specie conservation equations are solved by the successive substitution method. 

The iterative cycle is terminated when the normalized maximum correction, D max, is less than 

the convergence criterion (Cc). 

This option is called simultaneous option and is more efficient than nonsimultaneous 

option. The nonsimultaneous option, however, is more numerically stable. With the help of 

a logical variable. 
42

3. COMPUTER PROGRAM

The purpose of this section is to describe the structure and major functions of the program. 

The main objective of the computer program is to implement the numerical algorithm described in 

the previous section in a way that is easy to follow, modular to allow for future extension, robust, 

and free of errors. There are seven major functions of the computer program:

(1) Development of a flow circuit with fluid and solid nodes with branches and conductors.

(2) Development of an indexing system or data structure to define a network of fluid 

and solid nodes with branches and conductors.

(3) Generation of conservation equations of fluid mass, momentum, energy and species, 

and solid temperatures in respective nodes and branches.

(4) Calculation of thermodynamic and thermophysical properties of the fluid and solid 

in nodes.

(5) Numerical solution of conservation equations.

(6) Input/output.

(7) User-defined modules.

GFSSP consists of three major modules: 

(1) GUI.

(2) Solver and property (SP) module.

(3) User subroutine (US) module (fig. 25). 

Functions (1) and (6) are done in the GUI, functions (2)–(6) are done in the SP module, 

and function (7) is done in the US module. A distinct boundary is maintained among GUI, SP, 

and US modules. The GUI and US modules supply the information to the SP module through an 

input data file and USs. The SP module returns output data file and plot files for graphical and text 

display of results. The maintenance of a strict boundary among three modules is a key feature of 

GFSSP that makes the code easy to use, maintain, and upgrade. Users are not required to know 

the details of the computational method to become a proficient user of the code. The modularity 

also helps the developer to add new capabilities with minimum impact to the existing code. This 

section describes the SP and US modules. The GUI module is described in section 4.
43

• Creates Flow Circuit

• Runs GFSSP

• Displays Results Graphically

• Time-Dependent Process

• Nonlinear Boundary

• External Source Term

• Customized Output

• New Resistance/Fluid Option

• Equation Generator

• Equation Solver

• Fluid Property Program

Output Data File

Input Data 

File

Graphical User 

Interface 

(VTASC)

Solver and 

Property Module User 

Subroutines

New 

Physics

Figure 25. GFSSP process flow diagram showing interaction among three modules.

3.1 Process Flow Diagram

Figure 25 shows GFSSP’s process flow diagram to describe the interaction among GUI, 

SP, and US modules. Users create a flow circuit in the GUI VTASC by a ‘point, drag, and click’ 

method. VTASC creates an input data file that is read by the SP module. The user runs the SP 

module from VTASC, which also reads the output data file generated by the SP module to display 

the results in the GUI. The VTASC also allows users to plot time-dependent results in Winplot.13

Specialized input to the model can be supplied through user subroutines that also interact with the 

SP module. Such specialized input includes time-dependent processes; nonlinear boundary conditions; external mass, momentum, and energy sources; customized output; and new resistance and 

fluid options. 

3.2 Solver and Property Module

The main routine and the associated set of subroutines perform seven major functions 

that include the following: 

(1) Reading the input data file generated by VTASC, GFSSP’s graphical user interface.

(2) Generation of trial solution based on the initial guess.

(3) Supply time-dependent boundary conditions for unsteady flow.

(4) Numerical solution of conservation equations by SASS scheme.

(5) Interaction with thermodynamic property programs to calculate properties at nodes.

(6) Calculation of flow resistances in the branches.

(7) Create text output and plot files.

The flow charts of the SP module for nonsimultaneous and simultaneous schemes 

are shown in figures 26 and 27, respectively.
44

Obtain a Trial 

Solution

Simultaneous (S) 

solution of 

pressure, flow rate, 

and resident mass

NEWTON—Controls 

the iteration loop of 

Newton-Raphson 

scheme

ENTHALPY—Solves 

enthalpy 

by SS method

MASSC—Solves 

concentration 

by SS method

TSOLID or 

TSOLIDNR—Solves 

temperature by SS 

or S method

RESIST—Calculates 

resistances for all 

branches

KFACT1–KFSCT21 Calculates 

branch resistances

DENSITY—Calculates 

density and other 

properties

GASP and WASP or 

GASPAK—A routine for 

thermodynamic properties

UPDATE—After applying 

corrections update each 

variable (11)

Specie 

Conservation

WINPLOT—Generates 

plot files for WINPLOT

PLTVTASC—Generates 

plot files for VTASC

BOUND—Supply 

time-dependent 

boundary conditions

READIN—Reads 

input data from data 

file generated 

by VTASC

COEF—Constructs 

the coefficient matrix 

of the correction 

equation

SOLVE—Solve correction 

equation by Gaussian 

elimination method

EQNS—Calculates 

residuals of three 

governing 

equations

GASP and Wasp or 

GASPAK—A routine for 

thermodynamic properties

INIT—Generate 

trial solution

Successive 

substitution (SS) 

solution of enthalpy 

and concentration

Solution of 

temperature 

of solid

Calculate flow 

resistances

Convergence 

check

Time check

Stop

Start

Time Loop

Iteration Loop

Energy 

Conservation

Solid 

Temperature

Equation 

of State

Momentum 

Conservation

Mass 

Conservation

Figure 26. Flowchart of nonsimultaneous solution algorithm in solver and property module.
45

Obtain a Trial 

Solution

Implement SASS 

scheme by calling 

Subroutine 

NEWTON

NEWTON—Controls 

the iteration loop of

simultaneous 

solution scheme

 Newton-Raphson 

solution of pressure, 

flow rate, and 

resident mass

Successive 

substitution (SS) 

solution of enthalpy 

and concentration

Solution of 

temperature of solid

Calculate flow 

resistances

ENTHALPY—Solves enthalpy 

by SS method

MASSC—Solves concentration 

by SS method

TSOLID or TSOLIDNR—Solves 

temperature by SS or S method

RESIST—Calculates 

resistances for all branches

KFACT1–KFSCT21 Calculates 

branch resistances

DENSITY—Calculates 

density and other 

properties

GASP and WASP or 

GASPAK—A routine for 

thermodynamic properties

UPDATE—After applying corrections 

update each variable (11)

Specie 

Conservation

WINPLOT—Generates 

plot files for WINPLOT

PLTVTASC—Generates 

plot files for VTASC

BOUND—Supply 

time-dependent 

boundary conditions

READIN—Reads 

input data from data 

file generated 

by VTASC

COEF—Constructs 

the coefficient matrix 

of the correction 

equation

SOLVE—Solve 

correction equation by 

Gaussian elimination 

method

EQNS—Calculates 

residuals of three 

governing 

equations

GASP and WASP or 

GASPAK—A routine for 

thermodynamic properties

INIT—Generate 

trial solution

Time check

Stop

Start

Time Loop

Energy 

Conservation

Solid 

Temperature

Equation 

of State

Momentum 

Conservation

Mass 

Conservation

Figure 27. Flowchart of simultaneous solution algorithm in solver and property module.
46

3.2.1 Nonsimultaneous Solution Scheme

The flowchart of nonsimultaneous solution algorithm is shown in figure 26. In this scheme, 

there are two iteration cycles—outer and inner. In the inner iteration cycle, the mass and momentum conservation equation and resident equation of state are solved by the Newton-Raphson 

scheme until convergence. During this iteration cycle, property and resistance coefficients are not 

updated. In the outer iteration cycle, the energy and concentration equations are solved and density 

and resistance coefficients are updated and a new set of Newton-Raphson iteration is started. The 

outer iteration loop is repeated until the fractional change in density, resistance coefficient, and 

enthalpy is negligible.

The subroutine READIN reads the input data file. The subroutine INIT generates a trial 

solution by interacting with the thermodynamic property codes GASP, WASP, and GASPAK, 

or the property tables. Subroutine BOUND reads any applicable time-dependent boundary conditions from the model history files. Subroutine NEWTON conducts the Newton-Raphson solution 

of the mass conservation, flow rate, and energy conservation equations with the help of the subroutines EQNS, COEF, SOLVE, and UPDATE. The subroutine EQNS generates the equations. 

The coefficients of the correction equations are calculated in COEF. The correction equations are 

solved by the Gaussian elimination method in SOLVE. After applying the corrections, the variables are updated in subroutine UPDATE. This cycle of calculations is repeated until the corrections are negligible. The energy conservation equation is then solved in subroutine ENTHALPY or 

ENTROPY by the successive substitution scheme. For problems involving fluid mixture subroutine, 

MASSC is called to solve specie conservation equations. For conjugate heat transfer problems, the 

energy conservation equations for solid nodes are solved in subroutine TSOLID or TSOLIDNR. 

The resistance for each branch is calculated in RESIST following the calculation of fluid densities 

at each node in the subroutine DENSITY. The flow resistance coefficients (Kf ) for each branch 

are computed in subroutines KFACT1 through KFACT18, depending upon the resistance option 

selected for a particular branch. The convergence of the numerical scheme is checked to determine 

if the cycle of calculation needs to be repeated. The solver module also calls 15 user subroutines 

from various subroutines. 

3.2.2 Simultaneous Solution Scheme

The flowchart of the simultaneous solution algorithm is shown in figure 27. The functionality of subroutine READIN, INIT, and BOUND is identical to the nonsimultaneous scheme. In this 

scheme, there is only one iteration loop. The enthalpy (or entropy), concentrations, density, and 

resistance coefficient are updated in each Newton-Raphson iteration. Therefore, in each NewtonRaphson iteration, subroutine ENTHALPY or ENTROPY, MASSC, TSOLID or TSOLIDNR, 

RESIST, and DENSITY are called to compute and update all variables. The iteration loop is controlled in subroutine NEWTON. The interaction of SP and user subroutine module is identical to 

the nonsimultaneous scheme.

3.2.3 Conjugate Heat Transfer

GFSSP can model solid to fluid heat transfer which is commonly known as conjugate heat 

transfer. There are two solution options for solving the conservation equation for solid nodes—
47

successive substitution and Newton-Raphson. The successive substitution scheme is implemented 

in subroutine TSOLID and the flowchart of the subroutine is shown in figure 28. Subroutine 

TSOLID calls STCOND, CONVHC, RADHCF, RADHCSA, RADHCSSR, and SLDCP to estimate different terms of the conservation equation. STCOND determines the thermal conductivity 

from the property table. CONVHC determines the heat transfer coefficient from GFSSP’s built-in 

correlation. CONVHC also calls subroutine USRHCF to allow the user to provide a problemspecific heat transfer coefficient. RADHCF, RADHCSA, and RADHCSSR are three subroutines 

called from TSOLID to compute radiation heat transfer from solid to fluid, solid to ambient, and 

solid to solid, respectively. SLDCP determines the specific heat from property table for computing 

the transient term. Subroutine QDOTSSCR calculates solid to solid conduction and radiation heat 

transfer.

The Newton-Raphson scheme is implemented in subroutine TSOLIDNR. The flowchart of 

this subroutine is shown in figure 29. TSEQNS, TSCOEF, and GAUSSY are three subroutines that 

perform the major functions of Newton-Raphson scheme. The residuals are calculated in TSEQNS. The flowchart of subroutine TSEQNS is shown in figure 30. TSCOEF calculates the coefficient matrix. The correction equations are solved in GAUSSY. After convergence of the numerical 

scheme, subroutine QDOTSSCR calculates solid to solid conduction and radiation heat transfer.

 Subroutine TSOLID

Solid to solid heat transfer

Solid to fluid heat transfer

 Convection 

 Radiation

Solid to ambient 

heat transfer

Solid to solid radiation 

heat transfer

Transient term

Calculate solid temperature 

by successive substitution 

method

Check for 

convergence

Calculate solid to solid 

conduction and radiation heat

STCOND

Conductivity of solid

CONVHC

Solid to fluid heat transfer coefficient

RADHCF

Solid to fluid radiation heat transfer coefficient

RADHCSA

Solid to ambient radiation heat transfer coefficient

RADHCSSR

Solid to solid radiation heat transfer

SLDCP

Specific heat of solid

QDOTSSCR

Solid to solid conduction and radiation heat transfer

Figure 28. Solid temperature calculation by successive substitution method.
48

 Subroutine TSOLIDNR

Calculate residuals 

of conservation equations

Constructs coefficient matrix 

of correction equation

Solves the correction equation 

by Gaussian elimination method

Calculate new temperature 

by applying correction

Check for convergence

Calculate solid to solid 

conduction and radiation 

heat transfer

TSEQNS

TSCOEF

GAUSSY

QDOTSSCR

Solid to solid conduction 

and radiation heat transfer

Figure 29. Solid temperature calculation by Newton-Raphson method.

 Subroutine TSEQNS

Transient term

Solid to solid heat 

transfer

Solid to fluid heat 

transfer

 Convection

 Radiation

Solid to ambient heat 

transfer

Solid to solid radiation 

heat transfer

SLDCP

Specific heat of solid

STCOND

Conductivity of solid

CONVHC

Solid to fluid heat transfer coefficient

RADHCF

Solid to fluid radiation heat transfer coefficient

RADHCSA

Solid to ambient radiation heat transfer coefficient

RADHCSSR

Solid to solid radiation heat transfer

Figure 30. Calculation of residuals of energy conservation equation 

for Newton-Raphson method.
49

3.2.4 Thermodynamic Property Package

The thermodynamic property packages included in GFSSP consist of three separate programs—GASP,10 WASP,11 and GASPAK.12 It also includes tabulated data for RP-1. The GASP 

and WASP programs consist of a number of subroutines. GASP provides the thermodynamic 

properties for 10 fluids: helium, methane, neon, nitrogen, carbon monoxide, oxygen, argon, carbon 

dioxide, fluorine, and hydrogen. WASP provides the thermodynamic properties of water. RP-1 

properties are provided in the form of tables. Subroutine RP-1 searches for the required property 

values from these tables. GASPAK 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. 

The thermodynamic property subroutines are called from two GFSSP subroutines—INIT 

and DENSITY. In subroutine INIT, enthalpies and densities are computed from given pressures 

and temperatures at the boundary and internal nodes. In subroutine DENSITY, density, temperatures, specific heats, and specific heat ratios are calculated from given pressures and enthalpies 

at each node.
50

4. GRAPHICAL USER INTERFACE

This section introduces the VTASC, a unique GUI designed to simplify the model building 

process for GFSSP. VTASC allows the user to design GFSSP models using an interactive ‘point 

and click’ paradigm. The program seeks to eliminate some of the more tedious, error prone, and 

time-consuming operations associated with the model building process such as the selection of 

unique numbers for nodes, branches, and the explicit specification of the upstream and downstream nodes for every branch. The models may be easily modified both in terms of additional 

nodes, branches, and the model-specific data. Figure 31 shows the main VTASC window, which 

consists of menu and toolbar options and a blank canvas.

Figure 31. Main VTASC window.
51

4.1 Menus

4.1.1 File Menu

The file pull-down menu, shown in figure 32, contains the functions to begin a new model, 

open an existing model, save the model, save the model with an alternate location and name, print 

the model to a printer, print an image of the model to a bitmap (.bmp extension) file, write an input 

file for GFSSP based on the current model, and exit the application. The most commonly used of 

these functions are available, as shortcuts, from the file input/output toolbar. In addition, the File 

menu contains a listing of the nine most recently saved models.

Figure 32. VTASC file menu.

The new model function will reinitialize the application to a clean state without having to 

exit and then restart the application. If model data are present, then the user will be prompted to 

continue the operation. The open model function will present a file dialog to allow the user to select 

a previous model; note that all model files have a ‘.vts’ extension. The VTASC model files are not 

synonymous with the GFSSP input files. GFSSP compatible input files may be generated, based on 

the current model, as described below. The save function allows the user to save the current model 

to a desired location. In the case where the model has not been previously read or saved, a ̨file 

dialog will appear and allow the user to save the current model to a given location. The ‘save as’ 

function works identically except that a file dialog will appear in all instances. The print function 
52

produces a postscript file, which allows the user to print the current circuit to a printer or to a file. 

The print image function will save an image of the circuit in a bitmap file. The write GFSSP file 

function will become active once the user has input the required data; this is covered in the following section. The save and print functions are not available until at least one node is present.

4.1.2 Edit Menu

The edit pull-down menu, shown in figure 33, contains the functions to delete a selected 

item(s), activate the Global Options dialog, open an existing GFSSP output file using Notepad, 

and select all elements on the canvas. The delete function, which is not available until at least one 

node is present, also appears as a shortcut from the file input/output toolbar. The Global Options 

dialog is discussed in detail in section 4.2.

Figure 33. VTASC edit menu.

4.1.3 Advanced Menu

The advanced pull-down menu, shown in figure 34, contains the functions to activate dialogs for GFSSP’s advanced options such as Transient Heat, Heat Exchanger, Tank Pressurization, 

Turbopump, Valve Open/Close, Fluid Conduction, and Conjugate Heat Transfer. With the exception of Conjugate Heat Transfer, these functions are not available unless the advanced option has 

been activated through the Global Options dialog. The advanced option dialogs are discussed in 

detail in section 4.4.

Figure 34. VTASC advanced menu.
53

4.1.4 Run and Module Menus

The run pull-down menu contains the functions to call and run GFSSP and Winplot. Both 

of these functions also appear as shortcuts from the file input/output toolbar. Note that Winplot 

is not part of the GFSSP installation package and must be obtained separately. Also note that 

Winplot is not made available by VTASC unless a model is defined as an unsteady model. The 

module menu contains the function to activate the user executable Build dialog. The user executable Build dialog is discussed in detail in section 4.6.

4.1.5 Display, Canvas, and Help Menus

The display pull-down menu contains the functions to activate the Display Results/Properties dialog and clear any results/properties displayed on the canvas. The Display Results/Properties 

dialog is discussed in detail in section 4.8.2. The canvas pull-down menu contains the functions to 

toggle between a single- or a double-page canvas. The double-page canvas configuration is useful 

for larger models that will not easily fit on a single-page canvas. Note that in this configuration, 

a print margin is defined between the two canvas pages. The area within the print margin will not 

appear when a model hard copy is printed, so the user should avoid placing model elements in this 

area. The help pull-down menu contains the function to activate a pop-up window with additional 

information about that particular version of VTASC.

4.2 Global Options

From the Edit menu, select the menu option labeled Options... to display the Global 

Options dialog shown in figure 35. A left mouse click on items listed to the extreme left allows 

access to the desired information within the right pane. As shown, selecting the Instructions option 

gives general instructions on the use of this dialog.

4.2.1 General Information

Selecting the General Information item displays the following dialog (fig. 36) to access 

User Information, Solution Control, and Output Control. The User Information tab shown in 

figure 36(a) allows specification of the title, name of the analyst, working directory, GFSSP compatible input file, the output file to be generated by GFSSP, and the name and location of the 

GFSSP executable that will be used to run the model. The working directory, which is the directory 

where VTASC will write the GFSSP input and output files associated with the model, is assigned 

by specifying a file path for the GFSSP input file. The installation version of GFSSP is the default 

executable defined for any new model. Note that the file menu option ‘Write GFSSP File’ will 

become active only when the input and output GFSSP files have been specified. 

The Solution Control tab shown in figure 36(b) allows specification of certain characteristics 

of the solution procedure of a particular model. The user can choose either a simultaneous solution procedure or the original hybrid solution scheme. The user can also choose between the first 

law and second law of thermodynamics-based solution procedures for the energy equation. For the 

solid energy equation, the user can select either a Newton-Raphson or Successive Substitution 
54

Figure 35. Global Options dialog.

(a) 
55

(b) 

(c) 

Figure 36. General Information dialogs: (a) User Information tab, 

(b) Solution Control tab, and (c) Output Control tab.
56

solution scheme. In addition, the user can specify the convergence criteria, relaxation parameters, 

and choose the method by which the initial guess is made. Note that while there is an input box 

where the maximum number of iterations can be changed, this value is not currently being used 

by GFSSP. The user can also specify if they wish to use restart files by checking the appropriate 

box. Checking the Save Information box indicates that the user wishes to save the final solution 

to use as an initial guess in another model. Checking the Read Information box indicates that 

the user wishes to read in a previously saved solution as an initial guess. Note that if both boxes 

are checked, GFSSP will overwrite the initial restart files during the simulation. Two restart files 

are used by GFSSP for both saving and reading information. One file is used for node information and one for branch information. The user may name the restart files using the designated text 

boxes.

The Output Control tab shown in figure 36(c) allows specification of the type of data to 

output during the GFSSP simulation as well as offering the user the option to perform a sanity 

check on the GFSSP solution and flag any unreasonable values to the user. The solution check can 

be activated by choosing Check Values. The output options consist of inclusion (default) or suppression of (1) network information, (2) extended thermodynamic and thermophysical information 

at the nodes, (3) the initial flow field, and (4) Winplot data. For Winplot data, the user can define 

how often GFSSP writes output to the Winplot data files. (The default is 1, which writes output at 

every time step) as well as the type of file that is generated (comma separated value files or binary 

files.) 

Note that the Accept button must be pressed to accept modifications to the data in any of 

the tabs. To reset the data to default on all three tabs press Default and then the Accept buttons. 

The Cancel button closes the General Information item without accepting any changes. 

4.2.2 Circuit Options

Selecting the Circuit Options item displays the pane in figure 37, allowing access to Circuit 

Options and Initial Values. The Circuit tab shown in figure 37(a) allows specification of the options 

that will be activated for this circuit. These options include Inertia. (Note that while this option 

allows you to choose to supply relative angles between adjacent branches, no mechanism currently 

exists to define those angles using VTASC.) The user must manually edit the GFSSP input file to 

supply those angles. This is also discussed in section 4.5.13: Rotation, Gravity, Buoyancy (requires 

Gravity to be activated), Fluid Conduction (activates Fluid Conduction advanced option), Fluid 

Mass Injection, Heat Exchangers (activates Heat Exchanger advanced option), Turbopumps (activates Turbopump advanced option), Moving Boundary, Momentum Source, Axial Thrust, Normal 

Stress, Heat Source (with optional units), Cyclic Boundary, Solid-Fluid Heat Transfer Coefficient 

with the ability to choose method of calculation (only active for models using conjugate heat transfer), and Dalton’s Law of Partial Pressure (only active for multiple fluid models). Note that activating these options may require additional inputs in other areas of VTASC.

The Initial Values tab shown in figure 37(b) allows the user to set initial values of Pressure, 

Temperature, etc. for both boundary and interior nodes as well as solid nodes. The initial guess for 

all nodes may be changed, after the model has been built, by modifying the desired data and pressing the Apply to All button.
57

(a) 

(b) 

Figure 37. Circuit options dialogs: (a) Circuit tab and (b) Initial Values tab.
58

4.2.3 Unsteady Options

Selecting Unsteady Options displays the dialog shown in figure 38. The Unsteady Options 

dialog allows users to choose from various levels of unsteady modeling. The options in this window 

include: Steady (default), Quasi-Steady, Time Step, Start Time (a relative time—does not have to 

be zero), End Time, Print Interval (controls print interval to all output files except Winplot files), 

Unsteady, Variable Rotation (user specifies variable rotation file name), Variable Heat Load (activates Transient Heat advanced option), Variable Geometry (user specifies the variable geometry 

file name), Tank Pressurization (activates Tank Pressurization advanced option), and Valve Open/

Close (activates Valve Open/Close advanced option). Note that activating the unsteady options 

may require additional inputs in other areas of VTASC.

Figure 38. Unsteady Options dialog.

4.2.4 Fluid Options

Finally, selecting Fluid Options displays the dialog shown in figure 39. The Fluid Options 

dialog allows users to choose the thermodynamic property approach used in the model. The user 

can choose from the embedded thermodynamic property packages (1) GASP and WASP or

(2) GASPAK. Additionally, the user can choose a constant density fluid (the energy equation is not 
59

Figure 39. Fluid Options dialog.

calculated with this option and this option cannot be used with fully unsteady modeling); an ideal 

gas by specifying the fluid gas constant, specific heat, viscosity, thermal conductivity, reference 

pressure, reference temperature, reference enthalpy, and reference entropy; hydrogen peroxide, 

with the capability to define the mole fraction of water present in the fluid; and user-defined fluids 

through user-defined property files. Note that choosing the constant density fluid option causes the 

program to set the flow to steady and all unsteady options to be turned off. 

Desired fluids from the library of fluids may be selected and added to the Selected Fluids 

list by pressing the button. Fluids may be deleted from the Selected Fluids list by selecting the 

unwanted fluids and pressing the button. Note that a number enclosed in parentheses appears 

by each selected fluid. This is the GFSSP index number for that fluid. If a user fluid is selected, 

VTASC will prompt the user to double-click the fluid name to supply the fluid property file names 

and fluid molecular weight. Note that all user-defined property files must reside in the model’s 

working directory. The Switch Package button allows the user to switch between the two available 

thermodynamic property packages. The Switch Package button will only work when all selected 

fluids are common to both fluid libraries. Note that manually switching between the two thermodynamic property packages will delete all fluids from the selected fluids list.
60

4.3 Fluid Circuit Design

4.3.1 Boundary and Internal Node Properties

The boundary node addition tool located on the left border of VTASC is used to place 

boundary nodes on the drawing area, henceforth called the canvas. Upon selection of this tool 

the user may add boundary nodes by moving the mouse to the canvas and pressing the left 

mouse button. Placing a node will activate the delete function in the toolbar, edit menu, or by typing CTRL+D. Similarly, interior nodes may be added by selecting the interior node addition tool 

. Note that the nodes are automatically given unique numeric identifiers. Figure 40 shows 

a canvas with a number of boundary and interior nodes.

2

1 7

8

9

10

3 5

4

6

Figure 40. Boundary and internal nodes on canvas.

The selection tool is used to select the desired node and either modify its location or 

enable its deletion. Positioning the mouse and pressing the left button performs the selection; 

upon selection, the selected node will be shown with a red border. Repositioning a node is simply 

performed by pressing and holding the left mouse button over a node, moving the mouse to the 

desired location, and releasing the left mouse button. Multiple nodes may be selected for deletion 

by using the Ctrl keyboard button in conjunction with the mouse. The nodes may then be deleted 

using the Delete toolbar button as long as they are not attached to any branches. A press of the left 

mouse button, within the canvas, away from any nodes, will deselect any previously selected nodes. 

A right mouse button press upon a node will select the node and present a pop-up menu 

(fig. 41), allowing the user to delete the node (will not work on multiple node selections and is not 

available if the node is connected to any branch), set the properties for the indicated node, or align 

the node either horizontally or vertically with its neighboring elements (horizontal alignment aligns 
61

(a) 

(b) 

Figure 41. Node pop-up menus: (a) Boundary node and (b) internal node.

all elements to the right of the selected node, while vertical alignment aligns all elements below the 

selected node). Once a model has been run, the internal node pop-up menu also allows the user to 

activate the results dialogs, which will be discussed in section 4.8. 

Choosing the Properties ... option will present the dialog shown in figure 42. The appropriate inputs will be activated dependent upon the choices present within the Global Options dialog 

and the type of the selected node, whether boundary or interior. The user may modify the desired 

data within this dialog. To modify the concentration of a given fluid, select the desired fluid and 

type in the desired concentration. Note that directly upon selecting a fluid, the user may type without having to reposition the mouse. If the user wants to change the numeric identifier for a node, 

simply type in the desired numeric identifier (maximum of five numbers). The user may enter any 

desired descriptive text into the Node Description input box. Pressing the OK button will accept 

and adjust the revised data and the cancel button will reject the revised data. Also, in the case of 

an unsteady flow, each boundary node will be automatically assigned a unique Node History File 

name that is subject to user modification.

4.3.2 Branch Properties

The branch addition tool is used to specify the branches between the nodes. Selection of 

this tool immediately causes each interior and boundary node to be drawn with a series of ‘handles’ 

as denoted by the green squares in figure 43. 

The ‘handles’ serve to clearly identify eight possible locations on a node where the initial 

(upstream) or terminal (downstream) points of a branch may be located. Note that an unlimited 

number of branches may initiate or terminate at each ‘handle.’ There are two different types of 

branches that may be created: the first being a directed line segment between any two nodes while 

the second type is two possibly discontinuous line segments. In either case, a left mouse button click 

on a ‘handle’ will initiate a branch. Once an initial ‘handle’ (specifies the upstream node) has been 

selected, further movement of the mouse will draw a directed line segment from that ‘handle’ to the 
62

Figure 42. Node Properties dialog.

2

1 7

8

9

10

3 5

4

6

Figure 43. Nodes with branch ‘handles.’

current location of the mouse, as shown in figure 44(a). For the first type of branch, selecting 

another ‘handle’ completes the branch, as this second ‘handle’ effectively specifies the downstream 

node (fig. 44(b)).
63

(a) 

(b) 

Figure 44. Direct line segment branch: (a) Step 1 and (b) step 2.

The second type of branch is initiated identically; however, after selecting an initial ‘handle,’ 

an additional anchor point may be set at any location on the canvas by a press of the left mouse 

button at the desired location. The branch is then completed as usual by selecting another ‘handle.’ 

This series of steps is shown in figure 45(a) and (b).

(a) 

(b) 

Figure 45. Two-line segment branch: (a) Step 1 and (b) step 2.

Completing a branch will activate the delete function in the toolbar, edit menu, or by typing CTRL+D. As shown in figures 44(b) and 45(b), a unique numeric branch identifier formed by 

concatenating adjacent node numbers is automatically generated and the directed arrowhead visually defines the branch upstream and downstream node relationship. A click on the branch addition tool may be used to clear a branch currently under design and reinitialize the process. The 

selection tool may be used to select the desired branch and either modify its location or enable 

its deletion. The Ctrl button may be used in conjunction with the mouse to select multiple branches 

and nodes for deletion using the Delete toolbar button. Note that a node can only be deleted using 

this method if all branches attached to it are deleted as well.

Figure 46 shows an example fluid circuit complete with nodes and branches.
64

Specification of a new branch shows an image, which visually indicates that the resistance for this branch has not been specified. A right mouse click on the image will present the popup menu shown in figure 47, which allows the user to specify the Component; i.e., the resistance for 

the branch; align the branch either horizontally or vertically with its neighboring elements (horizontal alignment aligns all elements to the right of the selected branch, while vertical alignment 

aligns all elements below the selected branch); activate the Relocate ID dialog; activate the Change 

Branch Connection dialog; or delete the branch. The Properties option will be activated once a 

resistance has been selected for the branch and the Rotation/Momentum Data option will be activated in the case where either Rotation or Momentum is selected in the Global Options and 

2

1 7

8

9

10

3

5

4

6

17

78

79

39

910

810

410

510

? 106

?

? ?

? ?

?

?

?

?

Figure 46. Example fluid circuit with complete branch connections.

Figure 47. Branch pop-up menu.
65

for the specific branch under consideration. Once a model has been run, the branch pop-up menu 

also allows the user to activate the results dialogs, which will be discussed in section 4.8. 

Figure 48 shows the Relocate Branch ID dialog. This dialog gives the user a choice of eight 

locations where the branch identifier can be placed in relation to the branch element on the canvas. 

Figure 49 shows the change Branch Connections dialog. This dialog allows the user to 

change the nodes that the branch connects and/or the ‘handles’ where the branch attaches to each 

node. Note that changing the nodes connected to a branch does not automatically change the 

branch identifier.

Figure 48. Relocate Branch ID dialog.

Figure 49. Change Branch Connections dialog.

Choosing the ‘Component ...’ option will present the Resistance Options dialog shown in 

figure 50. This dialog shows pictorial representations for each of the 22 branch Resistance Options 

currently allowed in GFSSP. Note that the Control Valve option will not be available for steady 

flows. The available resistance options are discussed in section 2.1.7. To assign a resistance option, 

left click on the desired component and press the Accept button. Figure 51 shows an example 

where each of the branches has been assigned a resistance option. The user may change the resistance option for a branch at any time without deleting the branch.
66

Figure 50. Branch Resistance Options dialog.

2

1 7

8

9

10

3

5

4

6

17

78

79

39

910

810

410

510

106

2K

Figure 51. Example fluid circuit with resistance options.

A right click of the mouse on a branch where the resistance has been specified will present 

the pop-up menu (fig. 47) with the Properties... option activated. Choosing the Properties... option 

will present a dialog that is specifically tailored to receive input for that resistance option. In all 

instances, the properties specification dialogs behave in an identical fashion; however, the Fitting & 

Valves - dialog is somewhat different and will be shown as an example. Choosing the Properties... option, for a fitting and valve, will present the dialog shown in figure 52. The user can input 

the desired data or use the tree structure to the right to select a desired fitting or valve. Selection of 

a fitting or valve, from the tree, will load its specific data into the fields to the left; these data may 

then be edited as desired. The Accept button must be pressed to apply the data.
67

Figure 52. Fittings and valve resistance option properties dialog.

In general, the following applies to every resistance option. If the user wants to change the 

numeric identifier for a branch, simply type in the desired numeric identifier (maximum of five 

numbers). For the unsteady case, the initial flow rate may be specified, and depending upon the 

selected Global options, the Rotation, Momentum Source, and Inertia checkboxes may be active. 

Notice, in this case, that the Rotation checkbox is active and has been selected. If Rotation has 

been checked, a right click of the mouse button will present the pop-up menu with the Rotation/

Momentum Source option active and selecting this option will present the dialog shown in 

figure 53 to allow input of the relevant information.

Figure 53. Rotation/Momentum dialog.

4.3.3 Conjugate Heat Transfer

The Solid Node , Ambient Node , and Conductor addition tools are used to 

specify the conjugate heat transfer portion of a GFSSP model. These tools are inactive until the 
68

user activates Conjugate Heat Transfer by selecting it from the Advanced menu. In practice, adding 

solid and ambient nodes to the VTASC canvas is analogous to adding fluid boundary and internal 

nodes, while adding conductors is an identical process to adding fluid resistance branches. 

A right mouse click on an ambient node reveals a pop-up menu identical to the fluid boundary node pop-up menu shown in figure 41(a), while the solid node pop-up menu is identical to the 

fluid internal node pop-up menu shown in figure 41(b). The Properties... dialog for a solid node is 

shown in figure 54. For a solid node, all inputs are required at each node. A list of available materials is shown at the right of the dialog. When the user left mouse clicks the desired material from the 

list, the GFSSP index number for that material is automatically written to the Material input box. 

The ambient node properties dialog, shown in figure 55, requires only temperature as a modeling 

input.

The conductor pop-up menu, shown in figure 56, is very similar to the fluid branch popup menu (fig. 47). The only differences are that a Conductors option replaces the Components 

option—there is no need for a Rotation/Momentum option, and the Change Branch Connection

Figure 54. Solid Node Properties dialog.

Figure 55. Ambient Node Properties dialog.
69

Figure 56. Conductor pop-up menu.

Figure 57. Conductors dialog.

option is here named the Change Conductor Connection option. (The Functionality is identical to 

the Change Branch Connection option.) 

When the user first adds a conductor to the canvas, VTASC indicates that the conductor 

type is undefined using this symbol . The user defines the type of conductor by selecting the 

Conductors option from the conductor pop-up menu. This opens the Conductors dialog shown in 

figure 57. The user must select the appropriate type of conductor and click the Accept button. Note 

that VTASC will not allow the user to select a conductor type that is inconsistent with the types of 

nodes attached to that particular conductor; i.e., a solid-ambient convection conductor type cannot be applied to a conductor connecting two solid nodes. As with fluid resistance branches, once 

the conductor type has been defined, the Properties option becomes active in the conductor pop-up 

menu. Selecting the Properties option activates a Properties dialog where the user supplies characteristics that are specific to that conductor type.
70

4.4 Advanced Options

GFSSP contains many advanced features. These advanced features include: Transient Heat 

Load, Heat Exchanger, Tank Pressurization, Turbopump, Valve Open/Close, and Fluid Conduction. If any or all of the advanced features are selected, the user can input the appropriate information by selecting the option corresponding to the feature from the Advanced pull-down menu (see 

fig. 34).

The dialogs for each of the advanced options operate in an identical fashion. The user may 

add any number of components to that option by pressing the Add button. To modify the data for 

a particular component, the user must select the component of interest, modify the data, and press 

the Accept button. To delete a component, press the Delete button after a component has been 

selected.

4.4.1 Transient Heat

The Transient Heat Load option dialog is shown in figure 58. This option is activated from 

the Unsteady Options pane on the Global Options menu (see sec. 4.2.3). The user provides the fluid 

node where the heat load is applied and the name and location of a history file containing the heat 

load as a function of time.

Figure 58. Transient Heat Load option dialog.

4.4.2 Heat Exchanger

The Heat Exchanger option dialog is shown in figure 59. This option is activated from the 

Circuit Options pane on the Global Options menu (see sec. 4.2.2). The user has the option of modeling a counterflow or parallel flow heat exchanger. The user supplies the branch numbers that will 

be identified as the ‘hot’ and ‘cold’ branches as well as a value for heat exchanger effectiveness. If 

the user enters a heat exchanger effectiveness value between zero and 1, GFSSP will perform calculations based on that effectiveness. If the user enters a value >1, this prompts GFSSP to internally 

calculate the effectiveness. In the latter case, the user must also supply a value of UA (the product 

of the overall conductance for heat transfer and the surface area on which that conductance is 

based).
71

Figure 59. Heat Exchanger option dialog.

4.4.3 Tank Pressurization

The Tank Pressurization option dialog is shown in figure 60. This option is activated from 

the Unsteady Options pane on the Global Options menu (see sec. 4.2.3). The user has the option 

of modeling a vertically oriented cylindrical tank or a spherical tank. The user identifies the fluid 

nodes and fluid resistance branches in the model that represent the tank’s ullage and propellant. 

The user is also asked to provide the initial surface areas where the ullage is interacting with (a) the 

propellant and (b) the tank wall. The user must also supply relevant tank characteristics (density, 

specific heat, and thermal conductivity of the tank wall material; tank wall thickness; and initial 

tank wall temperature). Finally, VTASC provides default values for the constants used in GFSSP’s 

tank pressurization heat transfer calculations, but the user may modify these constants through this 

dialog if desired.

Figure 60. Tank Pressurization option dialog.
72

4.4.4 Turbopump

The Turbopump option dialog is shown in figure 61. This option is activated from the 

Circuit Options pane on the Global Options menu (see sec. 4.2.2). The user supplies the fluid 

resistance branches that will represent the pump and the turbine. The user also supplies some 

characteristics of the turbine (speed, efficiency, diameter, and design point velocity ratio) as well as 

the name and location of a history file containing the pump characteristics (the quantities (Head/

Speed2) and (Torque/(Density*Speed2)) as a function of (Flow Rate/Speed)).

Figure 61. Turbopump option dialog.

4.4.5 Valve Open/Close

The Valve Open/Close option dialog is shown in figure 62. This option is activated from the 

Unsteady Options pane on the Global Options menu (see sec. 4.2.3). The user provides the fluid 

resistance branch that represents the valve and the name and location of a history file containing 

the cross-sectional flow area of the valve as a function of time.

Figure 62. Valve Open/Close option dialog.
73

4.4.6 Fluid Conduction

The Fluid Conduction option dialog is shown in figure 63. This option is activated from 

the Circuit Options pane on the Global Options menu (see sec. 4.2.2). The user may populate the 

list of internal nodes in two ways. First, left mouse clicking the Load Nodes button will automatically populate the list with each internal node. Second, an individual internal node may be added 

to the list by typing the node identifier into the New Node input box and left mouse clicking the 

Add button. If the user wants to remove an internal node from the list, select that node from the 

list and left mouse click the Delete button. Selecting an internal node from the list reveals the list 

of upstream and downstream neighbors for that node. The user supplies the area and distance for 

each neighbor node by selecting that node from the neighbor nodes list.

Figure 63. Fluid conduction dialog.

4.5 GFSSP Input File

The primary interface between VTASC and GFSSP is the GFSSP input data file generated 

by VTASC. While it should not be necessary for the user to directly access the GFSSP input file for 

most modeling activities, knowledge of the content and format of the GFSSP input data file may 

be helpful in some circumstances. This section is intended to provide the user with a helpful reference for better understanding the GFSSP input data file.

4.5.1 Title Information

The GFSSP input data file always begins with the title and documentary information for the 

model, which are shown below. The first two entries are internally defined by VTASC. The first, the 

GFSSP version, tells GFSSP what features and formatting to expect while reading the input data 

file. The second entry defines where the user has installed GFSSP. The user defines the remaining 

four entries in VTASC. They include the analyst’s name, the working directory path and input data 

file name, output data file name, and a descriptive title for the model.
74

GFSSP VERSION

503

GFSSP INSTALLATION PATH

C:\Program Files\GFSSP\

ANALYST

ALOK MAJUMDAR

INPUT DATA FILE NAME

C:\Program Files\GFSSP\Examples\Ex1\Ex1.dat

OUTPUT FILE NAME

Ex1.out

TITLE

Simulation of a Flow System Consisting of a Pump, Valve and Pipe Line 

4.5.2 Logical Variables

The GFSSP input data file includes all 44 of the logical variable definitions based on the 

user’s choices in VTASC. The user is referred to appendix D (on CD inside back cover) for a specific definition of each logical variable.

USETUP

F

DENCON GRAVITY ENERGY MIXTURE THRUST STEADY TRANSV SAVER

F T T F F T F F

HEX HCOEF REACTING INERTIA CONDX ADDPROP PRINTI ROTATION

F F F F F F T F

BUOYANCY HRATE INVAL MSORCE MOVBND TPA VARGEO TVM

F T F F F F F F

SHEAR PRNTIN PRNTADD OPVALVE TRANSQ CONJUG RADIAT WINPLOT

F T T F F F F T

PRESS INSUC VARROT CYCLIC CHKVALS WINFILE DALTON

F F F F F T F

NORMAL SIMUL SECONDL NRSOLVT

F T T F

4.5.3 Node, Branch, and Fluid Information

This section of the GFSSP input data file defines the basic scope of the model, including 

(1) the total number of nodes, (2) the number of internal nodes, (3) the number of branches, 

and (4) the number of fluids.

NNODES NINT NBR NF

4 2 3 1

4.5.4 Solution Control Variables

The next section of the GFSSP input data file defines the numerical parameters chosen by 

the user, including the three under-relaxation parameters, the convergence criteria, and the maximum number of iterations.

RELAXK RELAXD RELAXH CC NITER

1 0.5 1 0.0001 500
75

4.5.5 Time Control Variables

This section of the GFSSP input data file is applicable only for unsteady models. It defines 

the time step, initial time, final time, output file print step, and the Winplot file print step.

DTAU TIMEF TIMEL NPSTEP NPWSTEP

1 0 200 25 1

4.5.6 Fluid Designation

This section of the GFSSP input data file lists the appropriate fluid definition information 

based on the user’s selections in VTASC. 

For a general fluid (GASP/WASP or GASPAK), the fluid designation lists the GFSSP index 

number for each selected fluid. 

NFLUID(I), I=1, NF

1 6 12

For a constant property fluid, the fluid designation lists the reference density and viscosity.

RHOREF EMUREF

62.4 0.00066

For an ideal gas fluid, the fluid designation lists the index number for an ideal gas and the 

reference properties associated with the ideal gas fluid.

NFLUID(I), I = 1, NF

 33

RREF CPREF GAMREF EMUREF AKREF PREF TREF HREF SREF

 53.34 0.24 1.3999 1.26e-05 4.133e-06 14.7 -459 0 0

For hydrogen peroxide, the fluid designation lists the index number and the mole fraction of 

water for the fluid.

NFLUID(I), I = 1, NF

 34

MFRAC

 0.5

Finally, for a user-defined fluid, the fluid designation lists the molecular weight of the fluid 

and the property table file names supplied by the user.

FLUID 1 PROPERTY FILES

28.0

AKFL1.DAT

RHOFL1.DAT

EMUFL1.DAT

GAMFL1.DAT

HFL1.DAT

SFL1.DAT

CPFL1.DAT
76

4.5.7 Node Numbering and Designation

The next section of the GFSSP input data file lists each node, designates whether that node 

is a boundary node (INDEX=2) or an internal node (INDEX=1), and includes a user-supplied text 

description of the node. Nodes are listed in the order that they are created in VTASC, which may 

not be numerical order.

NODE INDEX DESCRIPTION

 1 2 “Node 1”

 2 1 “Node 2”

3 1 “Node 3”

 4 2 “Node 4”

If the user has chosen to activate buoyancy, the reference node will be defined in this section 

right below the node listing.

REFERENCE NODE FOR DENSITY

2

4.5.8 Node Variables

The next section of the GFSSP input data file lists the initial properties at each node based 

on the user’s selections in VTASC.

For a steady state model, the model boundary conditions are listed along with the internal 

node initial guesses. Concentrations are listed at each node in the same order the fluids are listed 

(see sec. 4.5.6). If the user has chosen the Constant Property Fluid option, the temperature will not 

appear in this listing.

NODE PRES (PSI) TEMP(DEGF) MASS SOURC HEAT SOURC THRST AREA CONCENTRATION

 1 500 1500 0 0 0 0.1 0.9

 2 500 80 0 0 0 1 0

 3 338.2 1500 0 0 0 0.1 0.9

 4 14.7 80 0 0 0 0.5 0.5

For an unsteady model, the internal node initial solution values are listed first in the same 

order they were created. These properties include the node volume property, which does not appear 

in a steady state model. After all of the internal nodes have been listed, each boundary node history file is listed (again, in the order they were created). Each boundary node requires a separate 

history file. 

NODE PRES (PSI) TEMP(DEGF) MASS SOURC HEAT SOURC THRST AREA NODE-VOLUME CONCENTRATION

 1 100 80 0 0 0 17280 

 ex8hs2.dat
77

4.5.9 Transient Heat/Variable Geometry Information

This section of the GFSSP input data file is applicable only for unsteady models where the 

user has activated either the Transient Heat or Variable Geometry options. 

If the user has elected to use the Variable Geometry option, the Variable Geometry file 

name will appear right below the last boundary node history file name. There is no description line 

associated with the variable geometry listing in the input data file.

The transient heat section of the GFSSP input data file first lists the number of nodes identified by the user as having a transient heat load. Each identified node is then listed along with the 

corresponding heat load file name in the order the user added their information in the Transient 

Heat dialog (see sec. 4.4.1).

Transient Heat Load Information

Number of Nodes with Transient Heat Loads

1

Transient Heat Node Number

2

Corresponding Heat Load History File Name

qdot.dat

4.5.10 Node-Branch Connections

The next section of the GFSSP input data file identifies which branches are attached to each 

internal node. Each internal node (variable INODE) is listed in the same order they were created. 

The variable NUMBR defines how many branches are attached to that node. The array NAMEBR 

identifies which branches are attached to that node in the order they were created.

INODE NUMBR NAMEBR

2 2 12 23

3 2 23 34

4.5.11 Branch Flow Designation and Resistance Options

The next section of the GFSSP input data file describes the characteristics of each fluid 

branch. This section consists of two subsections. The first subsection identifies the upstream node, 

downstream node, and branch resistance option chosen by the user. This subsection also includes 

any text description of that branch supplied by the user. 

BRANCH UPNODE DNNODE OPTION DESCRIPTION

12 1 2 14 “Pump 12”

23 2 3 13 “Valve 23”

34 3 4 1 “Pipe 34”

The second subsection lists the properties of each fluid branch as defined by the user. The 

branch properties are specific to each branch resistance option. They are discussed in detail in section 2.1.7.
78

BRANCH OPTION-14 PUMP CONST1 PUMPCONST2 PUMPCONST3 AREA

12 30888 0 -0.0008067 201.06

BRANCH OPTION-13 DIA K1 K2 AREA

23 6 1000 0.1 28.274

BRANCH OPTION-1 LENGTH DIA EPSD ANGLE AREA

34 18000 6 0.005 95.74 28.274

4.5.12 Unsteady Information

The 12th section of the GFSSP input data file is applicable only for unsteady models. This 

section defines the initial mass flow rates in each fluid branch in the order the branches were 

created.

INITIAL FLOWRATES IN BRANCHES FOR UNSTEADY FLOW

1001 0.803

1002 0.803

1003 0.423

1004 0.423

4.5.13 Inertia Information

This section of the GFSSP input data file is applicable only for models where the user has 

activated the Inertia option. This section is divided into four subsections. The first subsection 

defines a fluid branch’s relationship with any upstream branches. Branches are listed in the order 

they were created. The variable NOUBR defines the number of upstream branches connected to 

a particular branch. The array NMUBR identifies which branches are attached upstream of that 

branch in the order they were created. 

The second subsection defines a fluid branch’s relationship with any downstream branches. 

Each branch is listed in the order they were created. The variable NODBR defines the number 

of downstream branches connected to a particular branch. The array NMDBR identifies which 

branches are attached downstream of that branch in the order they were created. 

BRANCH NOUBR NMUBR

12 0

23 1 12

34 1 23

BRANCH NODBR NMDBR

12 123

23 1 34

The third subsection allows the user to define relative angles between branches if desired. 

As discussed in section 5.2.2, the user must edit this section of the input data file manually. 

VTASC supplies a template for each branch in the order they were created. The template lists each 

upstream branch and corresponding angle first, then each downstream branch and corresponding angle. VTASC defines each angle as a placeholder value of zero degrees, which the user must 

replace with the appropriate angles. 
79

BRANCH

12

UPSTRM BR. ANGLE

DNSTRM BR. ANGLE

23 0.00000

BRANCH

23

UPSTRM BR. ANGLE

12 0.00000

DNSTRM BR. ANGLE

34 0.00000

The fourth subsection identifies the fluid branches where the user has activated inertia. The 

first number (16 in the example below) defines the number of branches where the user has activated 

inertia. The subsequent lines list each branch where the user has activated inertia in the order they 

were created.

NUMBER OF BRANCHES WITH INERTIA

16

12

23

34

45

56

67

78

89

910

1011

1112

1213

1314

1415

1516

1617

4.5.14 Fluid Conduction Information

This section of the GFSSP input data file is applicable only for models where the user has 

activated the Fluid Conduction option. This section is divided into two subsections. The first subsection identifies how many fluid conduction nodes the user has selected. For each of these fluid 

nodes, the number of upstream and downstream nodes connected to that node (identified here 

as neighbors) is defined and each neighbor node is listed.

NUMBER OF FLUID CONDUCTION NODES

2

NODE NO. OF NEIGHBORS NEIGHBOR NODES

2 2 1 3 

3 2 2 4 
80

The second subsection defines the fluid conduction properties (area and internode distance) 

for the interaction between a particular node and each of its neighbors. For each node, the property of interest between that node and each neighbor is listed in the order the neighbor nodes are 

listed in the first subsection.

NODE CONDUCTION AREAS

2 0.04909 0.04909 

3 0.04909 0.04909 

NODE INTERNODE DISTANCES

2 960 960 

3 960 960 

4.5.15 Rotation Information

This section of the GFSSP input data file is applicable only for models where the user has 

activated the Rotation option. This section first defines the number of branches where the user has 

activated rotation. Then, each branch is identified along with the rotational information for that 

branch (upstream and downstream node radii from center of rotation, rotational speed, and the 

rotational ‘slip’ factor). 

NUMBER OF ROTATING BRANCHES

9

BRANCH UPST RAD DNST RAD RPM K ROT

23 1.25 2.25 5000 0.8671

34 2.25 3.625 5000 0.8158

45 3.625 4.6875 5000 0.763

56 4.6875 5.375 5000 0.7252

67 5.375 5.5 5000 0.7076

89 5.5 5.375 5000 0.7129

910 5.375 4.6875 5000 0.7349

1011 4.6875 3.625 5000 0.7824

1112 3.625 2.65 5000 0.8376

4.5.16 Valve Open/Close Information

This section of the GFSSP input data file is applicable only for models where the user has 

activated the Valve Open/Close option. This section first defines the number of fluid branches 

where the user will be modeling a valve transient. Next, for each valve, the branch that will represent the valve and the valve history file name are listed.

NUMBER OF CLOSING/OPENING VALVES IN THE CIRCUIT

1

BRANCH

67

FILE NAME

ex15vlv.dat
81

4.5.17 Momentum Source Information

This section of the GFSSP input data file is applicable only for models where the user has 

activated the Momentum Source option. This section first defines the number of fluid branches 

where the user wishes to add a momentum source. Next, each branch where the user has defined 

a momentum source is listed along with the momentum source.

NUMBER OF BRANCHES WITH MOMENTUM SOURCE

1

BRANCH MOMENTUM SOURCE

12 100

4.5.18 Heat Exchanger Information

This section of the GFSSP input data file is applicable only for models where the user has 

activated the Heat Exchanger option. First, the number of heat exchangers identified by the user 

is defined. Then, the characteristics of each heat exchanger are listed as defined by the user, including the ‘hot’ and ‘cold’ branches, the type of heat exchanger (Counterflow=1, Parallel Flow=2), 

the ‘hot’ and ‘cold’ surface areas, UA, and the heat exchanger effectiveness. Note that the ‘hot’ and 

‘cold’ surface areas are not currently recommended for use in GFSSP and cannot be modified using 

VTASC.

NUMBER OF HEAT EXCHANGERS

1

IBRHOT IBRCLD ITYPHX ARHOT ARCOLD UA HEXEFF

23 67 1 0 0 1.1038 1.5

4.5.19 Moving Boundary Information

This section of the GFSSP input data file is applicable only for models where the user has 

activated the Moving Boundary option. This section defines the number of nodes identified as 

having moving boundary, and lists each identified node.

NUMBER OF NODES WITH MOVING BOUNDARY

2

NODE

1

2

4.5.20 Turbopump Information

This section of the GFSSP input data file is applicable only for models where the user has 

activated the Turbopump option. First, the number of turbopumps in the model is listed. Then, the 

characteristics for each turbopump (fluid branch representing the pump, fluid branch representing 

the turbine, speed, turbine efficiency, turbine diameter, design point velocity ratio, and the pump 

characteristics curve file name) are listed.
82

NUMBER OF TURBOPUMP ASSEMBLY IN THE CIRCUIT

1

IBRPMP IBRTRB SPEED(RPM) EFFTURB DIATRB PSITRD

23 1213 80000 0.5 3.435 0.4

PUMP CHARACTERISTICS CURVE DATA FILE

ex11pmp23.dat

4.5.21 Tank Pressurization Information

This section of the GFSSP input file is applicable only for models where the user has activated the Tank Pressurization option. First, the number of pressurized propellant tanks in the 

model is defined. Next, the characteristics of each tank are listed including the tank type (spherical=0, cylindrical=1); fluid node representing the ullage; ullage-propellant interface pseudo boundary node; fluid node representing the propellant; fluid branch representing the propellant surface; 

the initial tank wall surface area exposed to the ullage; the tank wall thickness; the tank wall material density, specific heat, and thermal conductivity; the ullage-propellant interface surface area; the 

heat transfer coefficient adjustment factor; the initial tank wall temperature; the heat transfer correlation ullage-propellant constants; and the heat transfer correlation ullage-tank wall constants.

NUMBER OF PRESSURIZATION PROPELLANT TANKS IN CIRCUIT

1

TNKTYPE NODUL NODULB NODPRP IBRPRP TNKAR TNKTH TNKRHO TNKCP

1 2 3 4 34 6431.9 0.375 170 0.2 

TNKCON ARHC FCTHC TNKTM CIP FNIP CIW FNIW

0.0362 4015 1 -264 0.27 0.25 0.54 0.25 

4.5.22 Variable Rotation Information

This section of the GFSSP input file is applicable only for models where the user has activated the Variable Rotation option. The variable rotation history file name is listed in this section.

ROTATION DATA FILE

varrot.dat

4.5.23 Restart Information

This section of the GFSSP input file is applicable only for models where the user has elected 

to read from and/or write to restart files. The section lists the node and branch restart file names.

RESTART NODE INFORMATION FILE

FNDEX15.DAT

RESTART BRANCH INFORMATION FILE

FBREX15.DAT
83

4.5.24 Cyclic Boundary Information

This section of the GFSSP input file is applicable only for models where the user has activated the cyclic boundary option. The section lists the boundary node where the cyclic boundary 

option has been activated and the node that is upstream of the cyclic boundary node.

CYCLIC BNDARY NODE UPSTREAM NODE

1 22

4.5.25 Conjugate Heat Transfer Information

This section of the GFSSP input file is applicable only for models where the user has 

activated Conjugate Heat Transfer. The section is divided into seven subsections.

The first subsection identifies how many solid and ambient nodes are present in the model 

as well as how many conductors of each type (solid-solid conduction, solid-fluid, solid-ambient, 

and solid-solid radiation) are present in the model.

NSOLID NAMB NSSC NSFC NSAC NSSR

8 2 7 8 2 1

The second subsection defines the characteristics of each solid node in the order they were 

created. The properties of each solid node (material, mass, and initial temperature) are listed first. 

Then, the number of conductors of each type attached to the solid node is listed along with a text 

description of the solid node provided by the user. Finally, for each type of conductor, every conductor attached to that solid node is listed.

NODESL MATRL SMASS TS NUMSS NUMSF NUMSA NUMSSR DESCRIPTION

2 41 1.00000 70.00000 1 1 1 1 “S Node 2”

NAMESS

23

NAMESF

122

NAMESA

12

NAMESSR

24

The third subsection lists each ambient node in the order they were created and lists the 

temperature at that node along with a text description of the ambient node provided by the user.

NODEAM TAMB DESCRIPTION

1 32.00000 “A Node 1”

The fourth subsection lists each solid-solid conduction conductor along with its characteristics (‘upstream’ solid node, ‘downstream’ solid node, surface area, distance, and a user-supplied text 

description) in the order they were created.
84

ICONSS ICNSI ICNSJ ARCSIJ DISTSIJ DESCRIPTION

23 2 3 3.14159 3.00000 “Conductor 23”

The fifth subsection lists each solid-fluid conductor along with its characteristics including: 

solid node; fluid node; heat transfer coefficient model (User Supplied=0, Dittus-Boelter=1, and 

Miropolskii=2); surface area, user-supplied heat transfer coefficient (if Model=0); emissivity 

of the solid; emissivity of the fluid; and user-supplied text description.

ICONSF ICS ICF MODEL ARSF HCSF EMSFS

122 2 12 0 1.88500e+01 3.17000e-04 0.00000e+00

EMSFF DESCRIPTION

0.00000e+00 “Convection 122”

The sixth subsection lists each solid-ambient conductor along with it characteristics (solid 

node, ambient node, surface area, heat transfer coefficient, emissivity of the solid, emissivity 

of the ambient, and a user-supplied text description).

ICONSA ICSAS ICSAA ARSA HCSA EMSAS

12 2 1 3.14159e+00 2.00000e-02 0.00000e+00

EMSAA DESCRIPTION

0.00000e+00 “Convection 12”

The seventh subsection lists each solid-solid radiation conductor along with its characteristics (‘upstream’ solid node, ‘downstream’ solid node, ‘upstream’ surface area, ‘downstream’ 

surface area, view factor, ‘upstream’ emissivity, ‘downstream’ emissivity, and a user-supplied text 

description).

ICONSSR ICNSRI ICNSRJ ARRSI ARRSJ VFSIJ EMSSI

24 2 3 3.14159 3.14159 1.00000 0.70000

EMSSJ DESCRIPTION

0.70000 “Conductor 24”

4.6 GFSSP Execution

As noted in section 4.1.4, GFSSP can be executed directly from the VTASC environment 

using either the Run menu or the shortcut on the file Input/Output toolbar. When the user activates 

the Run GFSSP command, VTASC automatically writes the GFSSP text input file before executing GFSSP. If the text file already exists, VTASC will ask if the user wishes to overwrite the file. 

VTASC then executes GFSSP and opens the GFSSP Run Manager window. The appearance and 

function of the Run Manager depends on whether the model is a steady state or transient model.

4.6.1 Steady State Run Manager

Figure 64 shows the GFSSP Run Manager appearance for a steady state simulation. If the 

user wishes to stop the GFSSP simulation for any reason, the button in the upper left corner of 

the Run Manager will stop GFSSP execution. During execution, the Run Manager will display 

GFSSP-generated messages in the GFSSP display pane. After execution is complete, the GFSSP 
85

messages may be printed by clicking the Print button in the lower right corner of the Run Manager. 

The GFSSP-generated text output file may also be viewed after execution is complete by clicking 

the Edit Output button in the lower left corner of the Run Manager. Clicking the Close button will 

exit the Run Manager and return the user to the VTASC window.

Figure 64. GFSSP steady state run manager.

4.6.2 Unsteady Run Manager

Figure 65 shows the GFSSP Run Manager appearance for an unsteady simulation. The only 

difference between the unsteady and steady Run Managers is that the GFSSP display pane appears 

in the top half of the Run Manager and a real-time updated plot of GFSSP’s convergence behavior 

as a function of time appears in the bottom half of the Run Manager.

Figure 65. GFSSP unsteady run manager.
86

4.7 GFSSP Output File

The basic GFSSP output for any simulation is the text output file. As mentioned in the previous section, once a simulation is complete the Run Manager gives the user the option of viewing 

the output file in a text editor. The content of the output file is dependent on the options selected 

by the user during VTASC model development. This section is intended to give the user an understanding of the format and general layout of the text output file.

4.7.1 Title and Data Files

Each GFSSP output file begins with the header shown below, which identifies the version of 

GFSSP that was used for the model simulation.

G F S S P (Version 5.0)

Generalized Fluid System Simulation Program

September, 2006

Developed by NASA/Marshall Space Flight Center

Copyright © by Marshall Space Flight Center

A generalized computer program to calculate flow

rates, pressures, temperatures, and concentrations

in a flow network.

 

Directly below the GFSSP header, the model title, analyst name, model working directory 

and text input file name, and text output file name as defined by the user are supplied.

TITLE:Simulation of a Flow System Consisting of a Pump, Valve and Pipe Line

ANALYST:ALOK MAJUMDAR

FILEIN:C:\Program Files\GFSSP\Examples\Ex1\Ex1.dat

FILEOUT:Ex1.out

4.7.2 Logical Variables

This section of the GFSSP output file lists the logical variable definitions as used in the 

simulation.

LOGICAL VARIABLES

DENCON = F

GRAVITY = T

ENERGY = T

MIXTURE = F

THRUST = F

STEADY = T

TRANSV = F

SAVER = F

HEX = F

HCOEF = F

REACTING = F
87

INERTIA = F

CONDX = F

TWOD = F

PRINTI = T

ROTATION = F

BUOYANCY = F

HRATE = T

INVAL = F

MSORCE = F

MOVBND = F

TPA = F

VARGEO = F

TVM = F

SHEAR = F

PRNTIN = T

PRNTADD = T

ADDPROP = F

PRESS = F

INSUC = F

VARROT = F

NORMAL = F

SECONDL = T

CONJUG = F

NRSOLVT = F

4.7.3 Node and Branch Information

This section of the output file documents the size and scope of the model. It lists the total 

number of nodes as well as the number of internal nodes, the number of branches, the number of 

fluids, the number of variables (or equations) in the model (this is the sum of NBR and NINT), 

and finally the enthalpy reference node, which is hard-coded in GFSSP as the second node.

NNODES = 4

NINT = 2

NBR = 3

NF = 1

NVAR = 5

NHREF = 2

4.7.4 Fluid Information

This section of the output file documents the fluids that were used in the simulation. 

For a constant property fluid, the fluid information lists the reference density and viscosity. 

RHOREF = 62.4000 LBM/FT**3

EMUREF = 0.6600E-03 LBM/FT-SEC
88

For all other fluid options, the fluid information lists each fluid in the order they were 

entered by the user.

FLUIDS: O2 H2O

4.7.5 Boundary Conditions

This section of the output file documents the boundary conditions of the model, which were 

supplied at each boundary node by the user. For a model with multiple fluids, the pressure, temperature, density, thrust surface area of the node, and concentration of each fluid at that node are 

listed. A single-fluid model provides the same listing with the exception of the fluid concentrations. 

A constant property fluid model lists only the pressure and surface area.

BOUNDARY NODES

NODE P T RHO AREA CONCENTRATIONS

(PSI) (F) (LBM/FT^3) (IN^2) O2 H2O

1 0.5000E+03 0.1500E+04 0.3931E+00 0.0000E+00 0.1000E+00 0.9000E+00

2 0.5000E+03 0.8000E+02 0.2819E+01 0.0000E+00 0.1000E+01 0.0000E+00

4 0.1470E+02 0.8000E+02 0.4725E+02 0.0000E+00 0.5000E+00 0.5000E+00

4.7.6 Fluid Network Information

This section of the output file is only active if the user has selected to print network information from the Global Options dialog (see sec. 4.2.1). For each internal node, it lists the thrust 

surface area, mass source, and heat source designated by the user. For each branch, it reprints the 

branch flow designation and resistance option information from the input text file (see sec. 4.5.11).

INPUT SPECIFICATIONS FOR INTERNAL NODES 

NODE AREA MASS HEAT

(IN^2) (LBM/S) (BTU/S)

2 0.0000E+00 0.0000E+00 0.0000E+00

3 0.0000E+00 0.0000E+00 0.0000E+00

BRANCH UPNODE DNNODE OPTION 

12 1 2 14

23 2 3 13

34 3 4 1

BRANCH OPTION -14: PUMP CONST1 PUMP CONST2 PUMP CONST3 AREA

12 0.309E+05 0.000E+00 -0.807E-03 0.201E+03

BRANCH OPTION -13: DIA, K1, K2, AREA

23 0.600E+01 0.100E+04 0.100E+00 0.283E+02

BRANCH OPTION -1: LENGTH, DIA, EPSD, ANGLE, AREA

34 0.180E+05 0.600E+01 0.500E-02 0.957E+02 0.283E+02

4.7.7 Initial Field Information

This section of the output file is only active if the user has selected to print the initial field 

from the Global Options dialog (see sec. 4.2.1). 
89

For each internal node of a single-fluid model, it lists the initial guesses for pressure and 

temperature as well as the resulting compressibility, density, and quality from the thermodynamic 

property calculations. For a multiple fluid model, the list is the same except that the quality is 

replaced with the initial guesses for the mass concentration of each fluid. For a constant property 

fluid model, only the initial guess for pressure is listed.

INITIAL GUESS FOR INTERNAL NODES

NODE P(PSI) TF(F) Z(COMP) RHO QUALITY

 (LBM/FT^3)

2 0.1470E+02 0.6000E+02 0.7616E-03 0.6237E+02 0.0000E+00

3 0.1470E+02 0.6000E+02 0.7616E-03 0.6237E+02 0.0000E+00

For each branch, the trial solution for the pressure drop across the branch and the mass 

flow rate in the branch is listed.

TRIAL SOLUTION

BRANCH DELP(PSI) FLOWRATE(LBM/SEC)

12 0.0000 0.0100

23 0.0000 0.0100

34 0.0000 0.0100

4.7.8 Conjugate Heat Transfer Network Information

This section of the output file is only active for models where the user has activated conjugate heat transfer and selected to print network information from the Global Options dialog (see 

sec. 4.2.1). For each solid node and conductor, the conjugate heat transfer information from the 

input text file is printed (see sec. 4.5.25).

CONJUGATE HEAT TRANSFER

NSOLIDX = 8

NAMB = 2

NSSC = 7

NSFC = 8

NSAC = 2

NSSR = 0

NODESL MATRL SMASS TS NUMSS NUMSF NUMSA

2 41 1.0000 70.0000 1 1 1

NAMESS

23

NAMESF

122

NAMESA

12

NODEAM TAMB

 1 32.0000

10 212.0000

ICONSS ICNSI ICNSJ ARCSIJ DISTSIJ

23 2 3 3.1416 3.0000

34 3 4 3.1416 3.0000
90

ICONSF ICS ICF ARSF EMSFS

122 2 12 18.8500 0.0000 0.0000

123 3 12 18.8500 0.0000 0.0000

ICONSA ICSAS ICSAA ARSA HCSA EMSAS EMSAA

12 2 1 0.3142E+01 0.2000E-01 0.0000E+00 0.0000E+00

910 9 10 0.3142E+01 0.2000E-01 0.0000E+00 0.0000E+00

4.7.9 Solution Results

This section of the output file documents the solution results of the GFSSP model. If the 

model is unsteady, a solution will be output at each time step the user has chosen to print (defined 

by the Print Frequency as discussed in sec. 4.2.3). The first line in the solution will list the current 

time step and the time at this step.

ISTEP = 25 TAU = 0.25000E+02

Next, the unsteady model will print out the boundary conditions at each boundary node 

for that time step. The format is identical to that discussed in section 4.8.5.

BOUNDARY NODES

NODE P(PSI) TF(F) Z(COMP) RHO QUALITY

 (LBM/FT^3)

2 0.1470E+02 0.8000E+02 0.1000E+01 0.7355E-01 0.0000E+00

After this line, the solution will be output in the same format for a steady model or the time 

step of interest in an unsteady model. For each internal node in a single fluid model, the calculated pressure, temperature, compressibility, density, resident mass, and fluid quality are listed. The 

listing is identical for a multiple fluid model except that the quality is replaced with the calculated 

mass concentration of each fluid at that node. For a constant property fluid model, only the calculated pressure and resident mass are listed. Note that for a steady model, the resident mass will 

always be zero. 

SOLUTION

INTERNAL NODES

NODE P(PSI) TF(F) Z RHO EM(LBM) QUALITY

 (LBM/FT^3)

2 0.2290E+03 0.6003E+02 0.1186E-01 0.6241E+02 0.0000E+00 0.0000E+00

3 0.2288E+03 0.6003E+02 0.1185E-01 0.6241E+02 0.0000E+00 0.0000E+00

If the user elects to print extended information in the Global Options dialog (see sec. 4.2.1), 

the output file will next list the calculated enthalpy, entropy, viscosity, thermal conductivity, specific 

heat, and specific heat ratio for each internal node. Note that this information will not be printed 

for constant property fluid models.

NODE H ENTROPY EMU COND CP GAMA

BTU/LB BTU/LB-R LBM/FT-SEC BTU/FT-S-R BTU/LB-R

2 0.2869E+02 0.5542E-01 0.7542E-03 0.9523E-04 0.1000E+01 0.1003E+01

3 0.2869E+02 0.5542E-01 0.7542E-03 0.9523E-04 0.1000E+01 0.1003E+01
91

For each branch, the calculated resistance factor, pressure drop, mass flow rate, velocity, 

Reynolds number, Mach number, entropy generation, and lost work are listed. 

BRANCHES

BRANCH KFACTOR DELP FLOW RATE VELOCITY REYN. NO. MACH NO.

(LBF-S^2/ (PSI) (LBM/SEC) (FT/SEC)

(LBM-FT)^2)

12 0.000E+00 -0.214E+03 0.191E+03 0.219E+01 0.241E+06 0.183E-02

23 0.764E-03 0.193E+00 0.191E+03 0.156E+02 0.644E+06 0.130E-01

34 0.591E+00 0.214E+03 0.191E+03 0.156E+02 0.644E+06 0.130E-01

ENTROPY GEN. LOST WORK

BTU/(R-SEC) LBF-FT/SEC

0.000E+00 0.000E+00

0.210E-03 0.848E+02

0.162E+00 0.657E+05

If the second law is used to solve the energy equation, the total entropy generation and work 

lost will be listed directly below the branch solution information.

***** TOTAL ENTROPY GENERATION = 0.163E+00 BTU/(R-SEC) *****

***** TOTAL WORK LOST = 0.120E+03 HP *****

If the user has activated conjugate heat transfer for a model, the conjugate heat transfer results will 

be listed next in the output file. 

For each solid node, the specific heat that was used is listed along with the calculated solid 

temperature. Note that the specific heat will be zero for steady models.

SOLID NODES

NODESL CPSLD TS

BTU/LB F F

2 0.000E+00 0.423E+02

3 0.000E+00 0.569E+02

For each solid-to-solid conductor, the thermal conductivity that was used is listed along 

with the calculated heat transfer rate.

SOLID TO SOLID CONDUCTOR

ICONSS CONDKIJ QDOTSS

BTU/S FT F BTU/S

23 0.261E-02 -0.333E-02

34 0.261E-02 -0.279E-02

For each solid-to-fluid conductor, the calculated heat transfer rate is listed along with the 

convection and radiation heat transfer coefficients that were used.

SOLID TO FLUID CONDUCTOR

ICONSF QDOTSF HCSF HCSFR

BTU/S BTU/S FT**2 F BTU/S FT**2 F

122 -0.115E-02 0.317E-03 0.000E+00

123 -0.544E-03 0.317E-03 0.000E+00
92

For each solid-to-ambient conductor, the calculated heat transfer rate is listed along with 

the convection and radiation heat transfer coefficients that were used.

SOLID TO AMBIENT CONDUCTOR

ICONSA QDOTSA HCSA HCSAR

BTU/S BTU/S FT**2 F BTU/S FT**2 F

12 0.448E-02 0.200E-01 0.000E+00

910 -0.136E-01 0.200E-01 0.000E+00

For each solid-to-solid radiation conductor, the calculated heat transfer rate is listed along 

with the effective conductivity that was used.

SOLID TO SOLID RADIATION CONDUCTOR

ICONSSR QDOTSSR EFCSSR

BTU/S BTU/S F

79 -0.113E-06 0.421E-08

If the user has requested that axial thrust be calculated in the Global Options dialog (see 

sec. 5.2.2), the calculated axial thrust will be listed next in the output file.

AXIAL THRUST = -527.30169 LBF

If the user has activated the turbopump advanced option (see sec. 4.4.4), the turbopump 

output will be listed next in the output file. First, the number of turbopumps in the model is listed. 

(Note that this value is not labeled in the output file.) Then, the pump branch, turbine branch, 

speed, the turbine efficiency at the design point, the turbine velocity ratio at the design point, the 

required torque, and the horsepower are listed for each turbopump in the model. 

1

IBRPMP IBRTRB SPEED(RPM) ETATRB PSITR TORQUE(LB-IN) HPOWER 

23 1213 0.800E+05 0.578E+00 0.269E+00 0.511E+02 0.649E+02

If the user has activated the pressurization advanced option (see sec. 4.4.3), the pressurization output will be listed next in the output file. First, the number of pressurization tanks 

in the model is listed. Then, the ullage node, propellant node, ullage to propellant heat transfer 

rate, ullage to tank wall heat transfer rate, tank wall conduction heat transfer rate, tank wall temperature, propellant volume and ullage volume are listed for each pressurization tank. Note that 

the labels for this output do not include the units. The units are Btu/s for the heat transfer rates, 

degrees Rankine for the tank wall temperature, and cubic feet for the volumes.

NUMBER OF PRESSURIZATION SYSTEMS = 1

NODUL NODPRP QULPRP QULWAL QCOND TNKTM VOLPROP VOLULG 

2 4 0.6644 2.1888 0.0000 195.6238 473.0886 26.9114
93

4.7.10 Convergence Information

The final section of the output file contains information on the convergence of the solution. 

It is important to remember that GASP/WASP allows extrapolation outside the stated limits 

of its fluid property relationships. While this allows for flexibility during the iterative process, it can 

occasionally lead to a final solution based on extrapolated properties. For models where the user 

has selected GASP/WASP, GFSSP checks the pressure and temperature at each node and prints 

a warning in the output file if they are outside of GASP/WASP’s stated limits so that the user can 

verify that the results are reasonable.

WARNING! CHKGASP: T out of fluid property range at node 1

WARNING! CHKGASP: T out of fluid property range at node 3

GFSSP also prints a statement indicating whether or not the solution converged. For an 

unsteady model, this statement is printed at each time step. If the solution converges, the statement 

lists the convergence criteria and the number of iterations needed to reach convergence.

SOLUTION SATISFIED CONVERGENCE CRITERION 0.100E-03 IN 5 ITERATIONS

If the solution does not converge, the statement lists the convergence criteria, the number of iterations that were performed, and the maximum difference after the last iteration. 

SOLUTION DID NOT SATISFY CONVERGENCE CRITERION 0.100E-02 IN 541 ITERATIONS

DIFMAX IN SUCCESSIVE ITERATION = 0.175E-02

If the model includes the cyclic boundary option (see sec. 4.2.2), the number of adjustment 

iterations and the final temperature difference are listed next in the output file.

ITERADJC = 3 DIFTEM = 1.634E-16

This section of the output file is only active if the user has selected to print the initial field 

from the Global Options dialog (see sec. 4.2.1). The time and time step will be listed after the convergence information. For a steady model, the time will read 100,000,000 s and the time step will 

be 1.

TAU = 100000000.000000 ISTEP = 1

The final section of the output file lists the CPU time to complete the model simulation.

TIME OF ANALYSIS WAS 1.001440000000000E-002 SECS

4.8 Post-Processing Simulation Data

While the GFSSP output file provides a comprehensive summary of model simulation 

results, it is not always a practical source of information to meet the user’s needs. Therefore, 

VTASC provides alternative methods of viewing GFSSP output for both steady and unsteady 

simulations.
94

4.8.1 Steady State Simulation Results

As mentioned in section 4.3, each GFSSP solution model element (fluid internal nodes 

and branches, conjugate heat transfer solid nodes and conductors) has a Results... dialog option 

located on their respective pop-up menus. After running a steady state simulation, if the user selects 

the Results option for a particular element, a table of results at that location will be displayed. 

Figure 66 shows the internal fluid node results table for a steady state simulation with multiple fluids. For this case, the table includes the calculated pressure, temperature, compressibility, 

density, resident mass, and the mass concentration of each fluid at that node. The internal fluid 

node table contents will vary just like the fluid node solution results discussed in section 4.8.9 based 

on the user’s selections (multiple fluid, single fluid, constant property fluid, print extended information). Clicking the OK button will close the Results option and return the user to VTASC.

The results tables for the fluid branches, solid nodes, and conductors are the same in appearance and function as the internal fluid node table discussed above. The parameters that are listed in 

each table are the same parameters discussed in section 4.8.9 for each respective element.

4.8.2 Unsteady Simulation Results

There are two options available in VTASC for generating plots of unsteady simulations, 

depending on the user’s needs—VTASC Plot and Winplot.

4.8.2.1 VTASC Plot. VTASC Plot is a built-in plotting capability within VTASC. As with 

the steady state results tables, it is accessed by selecting the Results... dialog from the desired model 

element’s pop-up menu. The appearance and function of VTASC Plot is the same for each model 

element. The only difference will be the parameters available to plot, which are the same parameters discussed for each element in section 4.7.9. 

Figure 67 shows the Results dialog for an unsteady simulation. Initially, the plot canvas 

space will be blank. The user can generate a hard copy or a bitmap of the desired plot by clicking 

the Print or Print to Bitmap button, respectively. The user creates the desired plot using the Properties dialog, which is activated by clicking the Properties... button. Once the user has finished, clicking the Close button will end the Results dialog and return the user to VTASC. 

Figure 68 shows the Properties dialog used to create a plot of unsteady results. The dialog 

consists of three tabs, as well as Apply, OK, and Cancel buttons. Clicking the Apply button accepts 

any changes that have been made. The OK button closes the Properties dialog and returns to the 

Results dialog. The Cancel button fulfills the same function as the OK button. 

The first tab (fig. 68(a)) allows the user to define the data they wish to plot. The user may 

select parameters for plotting using the Data Selection list at the top of the tab. The user selects 

plot parameters by highlighting the desired parameters in the Available Data list and clicking the 

 button. These parameters will then be added to the Selected Data list. If the user wishes to 

remove parameters from the Selected Data list, highlight those parameters and click the button. 
95

Figure 66. GFSSP steady state simulation results internal fluid node table.

Figure 67. GFSSP results dialog for unsteady simulation.

Note that VTASC Plot does not have Multi-Y axis plotting capability so scale should be considered 

when plotting multiple parameters on a single plot. The Data Properties list at the bottom of the 

tab can be used to design the line style of each plot parameter. The user highlights the parameter 

whose line style they wish to design and then selects the line type, width, and color from the available selections on the right. Once all changes have been made to a particular parameter, click the 

Apply button to accept the changes. 

The second tab (fig. 68(b)) allows the user to define the labeling parameters. Titles may be 

written or modified for the X and Y axes as well as the overall plot. The user may select whether 

or not they wish to include a grid or a legend on the plot by clicking the appropriate check box. 

The user may also define the number of minor tick marks they wish to see for each axis. Once all 

changes have been made, click the Apply button to accept the changes. 
96

(a) 

(b) 

(c) 

Figure 68. GFSSP VTASC plot properties dialog: (a) Data tab, (b) labeling tab, 

and (c) scale tab.
97

The third tab (fig. 68(c)) allows the user to modify the scale of each axis. By default, VTASC 

will auto scale a plot for the user. Deselecting the Auto Scale check box allows the user to define 

the minimum and maximum values for each axis. The user also has the option of converting the 

Y axis to a log scale. Once all changes have been made, click the Apply button to accept the 

changes. 

4.8.2.2 Winplot. If the user selects the Winplot plotting option from the Global Options 

dialog (see sec. 5.2.1), unsteady plot files will be generated in either comma delimited or binary 

formats. If the user selects the comma delimited option, several files are generated. The naming 

convention, description, and available parameters for each file are listed in table 9. Note that the 

conjugate heat transfer related output files are only written if a node or conductor of that type is 

present in the model. If the user selects the binary format, a single file with the name convention 

‘filename.WPL’ will be generated. This file will contain all of the available parameters for each 

model element as discussed in section 4.7.9 as well as DIFMAX, RSDMAX, and ITER.

Table 9. Winplot comma delimited unsteady output files.

Naming 

Convention Description Parameters

filenameFN.CSV Fluid node results P (psia), T (°F), Z, ρ (lbm/ft3), x of fluid or c of each fluid, μ (lbm/ft-s), 

k (Btu/ft-s-R), V (ft3), DIFMAX, RSDMAX, ITER

filenameB.CSV Fluid branch results v (ft/s), DP (psid), m (lbm/s), S

gen(Btu/R-s)

filenameSN.CSV Solid node results Cp,s (Btu/lbm-°F), Ts (°F)

filenameSF.CSV Solid-fluid conductor results hcsf (Btu/s-ft2-°F), hcsf,rad (Btu/s-ft2-°F) if needed, Q

sf (Btu/s)

filenameSS.CSV Solid-solid conductor results kss (Btu/ft-s-R), Q

ss (Btu/s)

filenameSA.CSV Solid-ambient conductor results hcsa (Btu/s-ft2-°F), hcsa,rad (Btu/s-ft2-°F) if needed, Q

sa (Btu/s)

filenameSSR.CSV Solid-solid radiation results keff,ssr (Btu/ft-s-R), Q

ssr (Btu/s)

As mentioned in section 4.1.4, Winplot must be obtained separately by the user. If Winplot is installed on the user’s computer, the user may open Winplot using the Run Menu’s Winplot 

selection or the Run Winplot button on the VTASC toolbar. If the unsteady plot files already exist, 

they will automatically be loaded into Winplot. Otherwise, the user must reselect the Run Menu 

Winplot selection or toolbar Winplot button to load the plot files into Winplot. For plotting and 

manipulating data in Winplot, the user is referred to the Winplot user’s manual.13

4.8.3 Display in Flow Circuit

An option exists within VTASC to display results for a particular model element on the 

VTASC canvas itself. The user can then observe how certain solution parameters change as adjustments are made to the model. This also allows the user to print the model with the results of interest as a hard copy or bitmap for reports or presentations. The user activates this option by first 

selecting the model element(s) using the selection tool. Then, as discussed in section 4.1.5, the user 

activates the Display Results/Properties dialog from the Display menu. Figure 69 shows an example 

of displaying the results on the VTASC canvas.
98

Figure 69. Display Results in flow circuit example.

Figure 70 shows the Display Results/Properties dialog. The user can select a maximum of 

three different parameters to display for each model element. Note that display parameters cannot 

be varied between like model elements (i.e., the user cannot display the pressure at one fluid node 

and the temperature at another fluid node). The selected parameters are displayed at all selected

elements of that type. Also note that the conjugate heat transfer parameter selections are not active 

unless the user has activated conjugate heat transfer. Once the user has selected all desired display 

parameters, the Apply button is used to accept the changes. The OK button closes the Display 

Results/Properties dialog and returns the user to VTASC. The Cancel button fulfills the same 

function as the OK button. 

The display units may be changed through the Display Property Units dialog, which is activated by clicking the Change button. Figure 71 shows the Display Property Units dialog. Alternative display units are available here for certain display parameters. The user toggles between the 

available units choices for each parameter. The OK button closes the Display Property Units dialog 

and returns the user to the Display Results/Properties dialog. The Cancel button fulfills the same 

function as the OK button.

Figure 70. Display Results/Properties dialog.
99

Figure 71. Display Property Units dialog.
100

5. EXAMPLES

This section demonstrates the major features of the code through 12 example problems, 

selected to serve two purposes. First, these problems will instruct the user on how to use the various 

options available in the code to analyze different fluid engineering problems. The other purpose of 

the examples contained within this section is to verify the code’s predictions. This verification was 

accomplished by comparing the GFSSP solutions with analytical solutions, other numerical solutions, or with test data. The included demonstration problems are: 

 (1) Simulation of a flow system consisting of a pump, valve, and pipe line.

 (2) Simulation of a water distribution network.

 (3) Simulation of compressible flow in a converging-diverging nozzle.

 (4) Simulation of the mixing of combustion gases and a cold gas stream.

 (5) Simulation of a flow system involving a heat exchanger.

 (6) Radial flow on a rotating radial disk.

 (7) Flow in a long-bearing squeeze film damper.

 (8) Simulation of the blowdown of a pressurized tank.

 (9) A reciprocating piston-cylinder.

(10) Power balancing of a turbopump assembly.

(11) Steady state and transient conduction through a circular rod, with convection.

(12) Simulation of fluid transient following sudden valve closure.

The selection of the order for these problems is primarily determined by their complexities. 

The first seven problems consider steady state flows and use relatively simple flow networks. Each 

example demonstrates use of a special option. Example 8 is a two-node model to demonstrate the 

use of the unsteady option. More complex unsteady flow examples are illustrated in examples 9, 10, 

and 12. Options for time-dependent geometry and moving boundary are shown in example 9. The 

application of the turbopump and heat exchanger options in a typical gas turbine system is illustrated in example 10. Example 11 introduces GFSSP’s conjugate heat transfer capability. Example 

11 is an adaptation of a classical heat transfer problem into GFSSP. Example 12 emphasizes 

GFSSP’s ability to predict fluid transient phenomena. Table 10 describes a matrix of the example 

problems and their use of various options to model the necessary physical processes. Several technical papers37–39 also illustrate the application of GFSSP in computing secondary flow and axial 

thrust in rocket engine turbopumps for both steady and unsteady flow. 

5.1 Example 1—Simulation of a Flow System Consisting of a Pump, Valve, and Pipe Line

5.1.1 Problem Considered

A problem commonly encountered in fluid engineering is to match a pump’s characteristics 

with the operating system’s characteristics. The designer needs to know the flow rate in the system 

and the power consumed by the pump. The following example problem demonstrates how GFSSP

can be used to simulate such problems. 
101

Table 10. Use of various options in example problems.

Example

Option

Conjugate Heat TransferConstant PropertyGravityHeat ExchangerIdeal GasLong InertiaMixtureMoving BoundaryNon Circular DuctPhase ChangePumpRotationTurbopumpUnsteadyVariable GeomWater Hammer

1 1 1

2 2

3 3

4 4

5 5

6 6 6

7 7 7 7

8 8 8

9 9 9 9

10 10 10

11 11

12 12 12

The system considered for this example is shown in figure 72. It consists of two reservoirs 

connected by 1,500 ft of 6-in-diameter pipe with a roughness factor (e/D) of 0.005. The receiving 

reservoir is located at an elevation that is 150 ft higher than the supply reservoir. The head-flow 

characteristics of the pump considered in this problem are shown in figure 73. This pump should 

be used to transport water, at 60 oF, from the supply reservoir to the receiving reservoir. GFSSP 

will be used to determine the system flow rate and required pump horsepower. 

5.1.2 GFSSP Model

The fluid system shown in figure 72 can be simulated with a GFSSP model consisting of 

four nodes and three branches as shown in figure 74(a). Nodes 1 and 4 are the boundary nodes 

representing the supply and receiving reservoirs that are both at 14.7 psia and 60 oF. Node 2 is an 

internal node representing the pump exit and the inlet to the gate valve. Node 3 is an internal node 

representing the exit from the gate valve and the inlet to the pipeline that connects the valve to the 

receiving reservoir. Branches 12, 23, and 34 represent the pump, gate valve, and pipeline, respectively. Figure 74(b) shows how this model appears in VTASC.

Once the boundary conditions have been established, the next step is to obtain the necessary 

information required to model the resistances and the momentum source located in the branches.

5.1.2.1 Branch 12 (Pump). Option 14 was selected to represent the pump because this 

option allows the user to model a pump with a given characteristics curve using either two or three 

constants. For this problem, two constants, A0 and B0, were input as well as the pump area (area 
102

Receiving Reservoir

Supply Reservoir

Gravity Vector

Pump

Gate Valve

L = 1,500 ft

D = 6 in

ε/D = 0.005

θ(wrt Gravity Vector)= 95.74°

Pipe

Flow

Water

θ

150 ft

Figure 72. Schematic of pumping system and reservoirs (Example 1).

Flow Rate (gpm)

Total Head

500

450

400

350

300

250

200

150

100

50

0

0 5,000 10,000 15,000 20,000

Figure 73. Manufacturer supplied pump head-flow characteristics.
103

Branch

Boundary Node

Internal Node

xx

12 23 34

xx

1 4

xx

2 3

1

12 23 34

C 2 3 4

L = 18,000 in

D = 6 in

ε/D = 0.03 in

θ(WRT Gravity Vector)= 95.74°

Gate Valve

D = 6 in

K1 = 1,000

K2 = 0.1

Pump

Intercept = 30,888

Slope = –8.067×10–4

A = 201.06177 in2

Supply Reservoir

P = 14.7 psia

T = 60 °F

Receiving Reservoir

P = 14.7 psia

T = 60 °F

Valve Outlet

Pipe Inlet

Pipe Outlet

Valve Inlet

Pipe

2 K

(a)

(b)

Figure 74. GFSSP model of pumping system and reservoirs: (a) Detailed schematic 

and (b) VTASC model.

is only required for the velocity calculation). These constants represent the slope and the intercept 

of the m 2 versus ∆p curve. The following procedure is used to obtain these constants: 

(1) Construct a table, based on user-selected points, from the pump characteristics curve 

shown in figure 73 to develop a relationship between m 2 and ∆p. These data are shown in table 11, 

where m =ρQ and Δ p = ρg

gc

H.

(2) Plot the ∆p and m 2 data from table 11 as shown in figure 75. Note that the relationship 

is linear (i.e., ∆p=A0+B0 m 2). Therefore, the pump characteristic curve can be prescribed with two 

constants, A0 and B0, and the optional third constant is not necessary.

(3) Determine the constants A0 and B0 from figure 75.

The intercept (∆p at m 2=0) is: 

A0=30,888. (63)

The slope (taken about Q=10,000 gpm) is:

B0 = 28,080 − 29,328

(2.784E + 6) − (1.2372E + 6)

= −(8.067E − 4) . (64)
104

Table 11. Tabulated pump characteristics data.

Q

(gpm)

m

(lb/s) 

Head 

(ft)

Δp

(psf)

m 2

(lb/s)2

–

4,000

8,000

12,000

16,000

20,000

–

556.13

1,112.3

1,668.4

2,224.5

2,781

495

485

470

450

425

385

30,888

30,264

29,328

28,080

26,520

24,024

–

3.093×105

1.2372×106

2.784×106

4.9484×106

7.734×106

0 1×106 2×106 3×106 4×106 5×106 6×106 7×106 8×106

35,000

30,000

25,000

20,000

15,000

10,000

5,000

0

2 (lb2/s2 m )

Delp (psf)

Figure 75. Pump characteristics curve in GFSSP format.

Since the curve is linear, the slope can be determined at any point without sacrificing accuracy. In 

VTASC, A0 is entered in the intercept input box and B0 is entered in the second order input box.

5.1.2.2 Branch 23 (Gate Valve). Option 13 was used to represent the gate valve. Option 13 

requires two constants, K1 and K∞ (two-K method), and the internal diameter to model various 

pipe fittings. The two required constants were obtained from table 7, assuming a reduced trim 

(b=0.8) gate valve).

5.1.2.3 Branch 34 ( Pipeline). Option 1 was used to represent the pipeline. This branch 

resistance option requires the user to supply the length, diameter, roughness factor (e/D), and 

the angle between the gravity vector and the pipeline.

5.1.3 Results

The example 1 GFSSP input and output data files (ex1.dat and ex1.out) are included 

in appendix D (on CD inside back cover). The example 1 GFSSP model predicts a flow rate 

of 191 lb/s and the pressure rise across the pump is 214 psi. Interpolation in the table 11 data shows 

that the pump pressure rise for the GFSSP predicted flow rate is 213 psi, which indicates 
105

that the model is working as expected. This example demonstrates that GFSSP can accurately 

predict the operating point of a fluid system consisting of a pump and a pipeline with a valve. 

When selecting a pump, the mass flow rate in the attached fluid system is generally 

unknown. By generating a system characteristic curve and plotting this curve against the pump 

characteristic curve, the operating point of the system can be determined. Using the example 1 

GFSSP model, the system characteristic curve can be generated in the following manner:

(1) Eliminate the pump by setting A0 and B0 to zero.

(2) Set boundary pressures P1 and P4 to desired values.

(3) Run the model.

(4) Repeat steps (2) and (3) to cover the desired range.

Table 12 shows typical system characteristics generated by performing the parametric study 

described above. The inlet pressures were arbitrarily selected to cover the expected operating range. 

The GFSSP predicted mass flow rates are shown in the third column. The flow resistance coefficients for the valve and the pipe are shown in the next two columns. The method of calculating 

these values is discussed in section 2. The next three columns of table 12 show hand-calculated 

pressure drops over the gate valve and the pipe and the pressure difference that is associated with 

the elevation change that exists in the example 1 model. The sum of these three losses is tabulated 

in the next column. The last column shows the overall pressure drop for the system for the given 

mass flow rate. A comparison of the values contained in the last two columns shows good 

agreement.

Table 12. Predicted system characteristics.

P1

(psia)

P4

(psia)

m

(lbm/s)

Kf, gate

lbf–s2/

(lb–ft)2

Kf, pipe

lbf–s2/

(lb–ft)2

Δp, gate

(psia)

Δp, pipe

(psia)

Δp, gate

(psia)

Δp, cal

(psia)

Δp, pres

(psia)

150

200

250

300

14.7

14.7

14.7

14.7

131

171

203

231

0.0019

0.0019

0.0019

0.0019

0.593

0.592

0.591

0.591

0.225

0.385

0.544

0.704

70.67

120.21

169.13

219.00

65.00

65.00

65.00

65.00

135.89

185.59

234.67

284.71

135.3

185.3

235.3

285.3

To determine the operating point of the system with a given pump, the system characteristics can be plotted together with the pump characteristics. The system operating point will be 

determined by the intersection of these two curves. Figure 76 shows the example 1 fluid system 

characteristics, generated previously, plotted along with the pump characteristics data that are 

shown in table 11. As seen from this figure, the operating point predicted by this method occurs at 

a pressure drop of approximately 214 psid and a mass flow rate of 190 lbm/s, which compares well 

with the predicted results from the GFSSP model containing the pump.
106

300

250

200

150

100

50

0

0 500 1,000 1,500 2,000 2,500 3,000

Pressure (psid)

Mass Flow Rate (lbm/s)

Pump Characteristics

Fluid System Characteristics

Figure 76. Fluid system operating point.

Finally, the power input to the pump can be calculated from the following relationship:

Horse power =

m

ρ

Δ p = (191 lbm/s)

62.4 lbm/ft

3 ( )

214 lbf/in2 ( ) 144 in2 /ft

2 ( )

(550 ft-lbf/hp) = 171 hp . (65)

5.2 Example 2—Simulation of a Water Distribution Network

5.2.1 Problem Considered

In Example 1, a single line pipe flow problem commonly encountered by pipeline designers is analyzed. In this example, an example associated with multipath systems, commonly known 

as flow networks, is considered. In general, water supply systems are considered as flow networks, 

since nearly all such systems consist of many interconnecting pipes. A 10-pipe (commercial steel) 

distribution system is shown in figure 77. Water at 50 psia enters the circuit at boundary node 1. 

Water is removed from the circuit at boundary nodes 3, 4, and 9 where pressures are also known. 

Use GFSSP to determine pressures at all of the remaining nodes and the flow rates in all of the 

pipes. The length, diameter, and roughness factors for each pipe are given in table 13. Note that 

pipes are designated as branches in the subsequent discussions.

5.2.2 GFSSP Model

The system shown in figure 77 is modeled by GFSSP in VTASC using 9 nodes and 

10 branches as shown in figure 78. The fluid was assumed incompressible and therefore a constant 
107

48 psia 45 psia

50 psia

1

46 psia

2

7

3 4

4

8

5 6

Figure 77. Water distribution network schematic (Example 2).

Table 13. Water distribution network branch data.

Branch

Length

(in)

Diameter

(in)

Roughness 

Factor

12

25

27

57

53

56

64

68

78

89

120

2,400

2,400

1,440

120

2,400

120

1,440

2,400

120

6

6

5

4

5

4

4

4

4

5

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

0.0018

48 psia 45 psia

50 psia

1

46 psia

2

7

3 4

4

8

5 6

Figure 78. GFSSP model of the water distribution network.
108

density (DENCON = .TRUE.) option was used. Nodes 1, 3, 4, and 9 are boundary nodes where 

the pressures are prescribed. Node 1 represents the inlet boundary node. Nodes 3, 4, and 9 are 

outlet boundary nodes. All of the remaining nodes (2, 5, 6, 8, and 7) are internal nodes where the 

pressures are calculated. All of the branches in this circuit simulate pipes. Therefore, each branch 

uses branch resistance option 1. The length, diameter, and roughness factors of all branches are 

given in table 13.

5.2.3 Results

The input and output data files (ex2.dat and ex2.out) are shown in appendix E (on CD 

inside back cover). The GFSSP-predicted results are shown in tables 14 and 15. Table 15 lists the 

predicted pressures at the internal nodes in psia and feet of water. Table 15 lists the predicted flow 

rates in lbm/s and ft3/s. Table 15 also provides a comparison between the GFSSP predicted results 

and values predicted by the Hardy Cross2,26 method.

Table 14. GFSSP predicted pressure distribution at the internal nodes.

P2 P5 P6 P7 P8

psia ft psia ft psia ft psia ft psia ft

49.80 114.92 48.11 111.02 45.34 104.63 48.35 111.58 46.01 106.18

Table 15. GFSSP and Hardy Cross method predicted branch flow rates.

Flow Rate

GFSSP Hardy Cross

lb/s ft3/s lb/s ft3/s

m 12 100.16 1.605 100.16* 1.605*

m 25

63.1 1.011 63.59 1.019

m 27

37.0 0.593 36.58 0.5862

m 53

44.43 0.7115 44.43* 0.7115*

m 56

29.1 0.466 29.11 0.4665

m 57

–10.4 –0.167 –9.93 –0.1592

m 64

47.07 0.7548 47.07* 0.7548*

m 68

–18.0 –0.288 –17.99 –0.2883

m 78

26.7 0.428 26.64 0.4270

m 89

8.66 0.1387 8.66* 0.1387*

* Boundary flow rates are prescribed from GFSSP predictions.

Figure 79 shows a comparison between GFSSP and Hardy Cross predicted flow rates. The comparison appears reasonable considering the fact that the Hardy Cross method assumes 

a constant friction factor in the branch while GFSSP computes the friction factor for all branches. 

Therefore, as the flow rates change, the friction factor also changes.
109

1.6

1.4

1.2

1

0.8

0.6

0.4

0.2

0

–0.2

–0.4

Flow Rate (ft3/s)

Q - GFSSP Solution (ft3/s)

Q - Hardy Cross Prediction (ft3/s)

12 25 27 53 56 57 64 68 78 89

Branch Number

Figure 79. A flow rate comparison between GFSSP and hardy cross method predictions.

5.3 Example 3—Simulation of Compressible Flow in a Converging-Diverging Nozzle

5.3.1 Problem Considered

In the previous examples, incompressible flows in fluid systems were considered. In this 

example, compressible flow in a converging-diverging nozzle demonstrating GFSSP’s capability to 

handle compressibility will be considered. One of the characteristics of the compressible flow in 

a duct is that the flow rate becomes independent of exit pressure after reaching a threshold flow 

rate. This threshold value is known as the choked flow rate and it is a function of inlet pressure and 

temperature. Flow in a confined duct becomes choked when the flow velocity equals the local velocity of sound. The purpose of this example is to investigate how accurately GFSSP can predict the 

choked flow rate in a converging-diverging nozzle.

The converging-diverging nozzle considered for this example is shown in figure 80. 

The nozzle is 6.3 in long with a 0.492-in-diameter throat. The inlet diameter of the nozzle 

is 0.6758 in and throat of the nozzle is located 0.158 in downstream of the inlet. The fluid considered was steam at 150 psia and 1,000 oF. The nozzle backpressure was varied from 134 to 45 psia. 

GFSSP is required to predict the flow rate and the pressure distribution for different exit 

pressures. The predicted flow rate will also be compared with the isentropic solution. 

5.3.2 GFSSP Model

The fluid system shown in figure 80 can be simulated with a GFSSP model consisting 

of 17 nodes and 16 branches as shown in figure 81(a). Nodes 1 and 17 are the boundary nodes 

representing the inlet and outlet of the nozzle. All of the remaining nodes are internal nodes connected in series. GFSSP can be used to construct an isentropic model by selecting branch resistance 

option 2 (flow through a restriction), using a flow coefficient, CL set equal to zero and by setting 

the logical flag INERTIA = .TRUE. This option eliminates friction from the momentum equation, 

which represents a balance between the inertia and pressure force with the inclusion of the inertia 

term. Each branch assumes a constant flow area that was determined from the nozzle geometry 

at the midpoint of the branch location. The branch information is listed in table 16. Figure 81(b) 

shows how this model appears in VTASC. It may also be noted that the energy equation has been 

solved using the second law option.
110

Flow

6.142 in

0.158 in

Not to Scale

0.492 in 1.453 in

Figure 80. Converging-diverging steam nozzle schematic (Example 3).

Inlet Plane Throat

0.0296 in

×2

0.0395 in

×4

0.4424 in

×2

0.8848 in

×5

0.1658 in

×2

0.4424 in

Node 

Number

Exit Plane

0.246 in

0.7265 in

R = 0.158 in Not to Scale

Flow

6.142 in

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

1 17

2

12 23

3

R R

4

34 45

5

R R

6

56 67

7

R R

8

78 89

9

R R

10

910 1011

11

R R

12

1112 1213

13

R R

14

1314 1415

15

R R

16

1516 1617

R R

(a)

(b)

Figure 81. Converging-diverging steam nozzle model: (a) Detailed schematic 

and (b) VTASC model.
111

Table 16. Converging-diverging nozzle branch information.

Branches Option Area (in2)

12

23

34

45

56

67

78

89

910

1,011

1,112

1,213

1,314

1,415

1,516

1,617

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

0.3587

0.2717

0.2243

0.2083

0.1901

0.1949

0.2255

0.2875

0.3948

0.5640

0.7633

0.9927

1.2520

1.4668

1.5703

1.6286

It may be mentioned here that the temperature at the outlet boundary node has no influence 

on the solution. However, the code must be supplied with a value to satisfy the input requirements 

for a boundary node. In this example the temperature was specified arbitrarily to 1,000 oF at 

node 17. The boundary conditions are specified in table 17.

Table 17. Converging-diverging nozzle boundary conditions.

P1

(psia)

T1

(°F)

P17

(psia)

T17

(°F)

150

150

150

150

150

1,000

1,000

1,000

1,000

1,000

134

100

60

50

45

1,000

1,000

1,000

1,000

1,000

5.3.3 Results

The outlet boundary node pressures were varied to include 134, 100, 60, 50, and 45 psia. 

The input and output files (ex3.dat and ex3.out) from the example case using an exit pressure of 

60 psia are included in appendix F (see CD inside back cover). The predicted pressure distributions 

for five test cases are shown in figure 82. Figure 83 shows the predicted temperature distributions 

for the same cases. 

Table 18 lists the model predicted mass flow rates with varying exit pressures. As expected, 

the mass flow rate increased as the exit pressure was decreased until the pressure ratio decreased 

below the critical pressure ratio. At this point and below, the mass flow rate remained constant due 

to choking of the flow at the nozzle throat.
112

160

140

120

100

80

60

40

20

0

Pressure (psia)

0 1 2 3 4 5 6 7

Axial Position (in)

Exit Pressure = 130 psia

Exit Pressure = 100 psia

Exit Pressure = 60 psia

Exit Pressure = 50 psia

Exit Pressure = 45 psia

Figure 82. Predicted pressures for the isentropic steam nozzle.

1,100

1,000

900

800

700

600

500

Temperature (°F)

0 1 2 3 4 5 6 7

Axial Position (in)

Exit Pressure = 130 psia

Exit Pressure = 100 psia

Exit Pressure = 60 psia

Exit Pressure = 50 psia

Exit Pressure = 45 psia

Figure 83. Predicted temperatures for the isentropic steam nozzle.
113

Table 18. Predicted mass flow rate with varying exit pressure.

Pexit

(psia)

m

(lbm/s)

134

100

60

50

45

0.279

0.329

0.336

0.337

0.337

The isentropic flow rate was calculated from equation (66). This equation assumes that the 

inlet pressure is a stagnation pressure. GFSSP’s formulation assumes that the prescribed boundary 

conditions are taken at a static condition. The nozzle inlet velocity head component must be added 

to the GFSSP static inlet boundary pressure to obtain the correct nozzle inlet stagnation pressure 

to use in equation (66) as shown in equation (67): 

m = Athroat Pinlet

gc γ

R T inlet

γ +1

γ −1

⎛

⎝

⎜

⎞

⎠

⎟ 2

γ + 1

⎛

⎝

⎜

⎞

⎠

⎟ (66)

and

Pinlet = Pstatic 1 + γ −1

2

⎛

⎝

⎜ ⎞

⎠

⎛ ⎟ M2

⎝

⎜ ⎞

⎠

⎟

γ

γ −1 . (67)

A value of the specific heat ratio at the nozzle inlet boundary node was obtained from 

GASP. Then, using the case where the exit pressure is 60 psia (the model output data file contained 

in appendix F (on CD inside back cover)), the inlet Mach number is obtained from the first branch 

connected to the inlet node. Substituting these values into equation (67) gives the following total 

pressure at the nozzle inlet:

Pinlet = (150 psia) 1 + 1.2809 −1

2

⎛

⎝

⎜ ⎞

⎠

⎟(0.342) ⎛ 2

⎝

⎜ ⎞

⎠

⎟

1.2809

1.2809−1 = 161.6 psia . (68)

Substituting the calculated total inlet pressure from equation (68) into equation (66) 

and solving for the choked mass flow rate gives a calculated isentropic choked mass flow rate 

of 0.327 lbm/s as shown:

m = 0.19012 in2 ( ) 161.6 lbf

in2

⎛

⎝

⎜ ⎞

⎠

⎟

32.174 lbm-ft

lbf-s

2 (1.281)

85.83 lbf-ft

lbm-oR

1,460 oR

2

1.281+1

⎛

⎝

⎜ ⎞

⎠

⎟

2.281

0.281 = 0.327 lbm

s . (69)
114

As the reader can see by comparing the results shown in table 18 and equation (69), there is 

good agreement (approximately a 3% difference) between the GFSSP predicted choked flow rate 

and the calculated isentropic choked flow rate. The prediction from the first law-based formulation 

has also been compared with the second law formulation. In the first law formulation, enthalpy 

was used as the dependent variable and equation 4(a) was used instead of equation 4(b). Unlike 

the second law formulation where density is computed from pressure and entropy, the first law 

formulation calculates density from pressure and enthalpy. The comparison of choked mass flow 

rate for both formulations with the isentropic solution is shown in table 19. It may be noted that 

the entropy-based formulation more accurately predicts the mass flow rate than the enthalpy-based 

formulation. 

Table 19. Comparison of choked mass flow rates.

Parameter

Second Law 

Formulation

First Law 

Formulation

GFSSP predicted isentropic mass flow rate 0.337 lbm/s .308 lbm/s

Calculated stagnation inlet pressure 161.6 psia 159.71 psia

Calculated isentropic mass flow rate 0.323 lbm/s 0.323 lbm/s

Percent difference between GFSSP predicted 

and calculated isentropic mass flow rates 

3.06% – 4.64%

To validate GFSSP’s ability to predict temperatures, the GFSSP predicted temperature at 

the nozzle throat was also compared with a hand-calculated temperature using equation (70). This 

equation assumes an isentropic process. In equation (70), the pressure P1 was assumed to be the 

total pressure at the nozzle inlet that was calculated from equation (68) and the throat pressure was 

taken from the GFSSP output:

T2 = T1

P2

P1

⎛

⎝

⎜

⎞

⎠

⎟

γ −1

γ

 = (1,460 °R) 82.13 psia

161.6 psia

⎛

⎝

⎜

⎞

⎠

⎟

1.2809–1

1.2809 = 1,258.6 °R = 798.6 °F . (70)

The GFSSP predicted nozzle throat temperature is 819.9 °F. This temperature compares 

well, within 3%, with the value calculated from equation (70).

The isentropic steam nozzle model represents a reversible process in which no entropy 

is generated. In an actual steam nozzle there will be frictional losses that result in an increase in 

entropy. To demonstrate GFSSP’s ability to predict the irreversibility (entropy generation) of 

a process, the isentropic steam nozzle model was modified to allow for frictional losses. The branch 

resistance option for each model branch was changed from Option 2 with a CL =0, to Option 1 

using the branch lengths and diameters shown in figure 81 and assuming an absolute roughness 

of 0.01 in.

Figure 84 shows a temperature/entropy comparison between the GFSSP isentropic steam 

nozzle model and the GFSSP steam nozzle model with frictional losses. As one can see in figure 84, 
115

1,100

1,000

900

800

700

600

500

Temperature (°F)

1.8 1.82 1.84 1.86 1.88 1.9 1.92 1.94 1.96 1.98 2

Entropy (Btu/lbm-R)

Isentropic

Irreversible

Figure 84. Temperature/entropy plot comparing the isentropic steam nozzle 

with an irreversible process.

the isentropic model predicts no change in entropy through the nozzle while the irreversible process 

predicts a increase in entropy of 0.019 Btu/(lbm °R). 

5.4 Example 4—Simulation of the Mixing of Combustion Gases and a Cold Gas Stream

5.4.1 Problem Considered

In the previous examples, the fluid systems that were considered employed a single fluid. 

In this example, simulation of multiple fluids in a mixing process are considered. How to use the 

MIXTURE logical option in the code to simulate the mixing of combustion gases and a cold gas 

stream by utilizing the flow system shown in figure 85 will be demonstrated. A mixture of hot 

combustion products, consisting of water vapor and oxygen, is mixed with cooler oxygen gas. The 

mixture temperature and composition are required to be calculated. 

A mixture consisting of 90% water vapor and 10% oxygen (by mass) at 500 psia and 

1,500 oF mixes with pure oxygen at 500 psia and 80 oF. GFSSP is required to predict the flow rate, 

mixture temperature, and composition of the mixture. A hand calculation will also be performed 

of the mixture temperature and the composition of the mixture to verify GFSSP’s predictions. 

5.4.2 GFSSP Model

The mixing chamber shown in figure 85 can be simulated with a GFSSP model consisting of 

four nodes and three branches as shown in figure 86. Nodes 1, 2, and 4 are the boundary nodes 
116

500 psia

1,500 °F

XH2O = 0.9

XO2

 = 0.1

CpH2O = 0.5696 (btu/(lbm °F))

CpO2

 = 1.2512 (btu/(lbm °F))

500 psia

80 °F

XO2

 = 1.0

CpO2

 = 0.2324 (btu/(lbm °F))

14.7 psia

80 °F

23

3

1

4

2

m

34 m

13 m

Figure 85. Mixing problem schematic (Example 4).

1

2

13

23

34

3 4

R

R

C O

Figure 86. GFSSP model of mixing problem.

representing the inlet and outlet of the mixing chamber and node 3 is the internal node representing the mixing chamber. Branches 13 and 23 are represented by option 2 using a flow coefficient 

of 0.6 and area of 1 in2. Branch 34 also uses a flow coefficient of 0.6 and area of 1 in2, but it is 

modeled using option 22. The reason option 22 is used for branch 34 is to account for the possibility of choked flow in the branch.

5.4.3 Results

The input and output files of this example are included in appendix G as ex4.dat and ex4.

out (see CD inside back cover). The predicted flow rates in branches 13, 23, and 34 are 1.16, 3.10, 

and 4.26 lbm/s, respectively. The predicted temperature at the outlet of the mixing chamber, node 3, 

is 689 oF and the composition is 24.47% water vapor and 75.53% oxygen. The mixture will not vary 

between nodes 3 and 4.
117

Now verify the predicted results by performing hand calculations of the mixing process. The 

temperature of the mixture can be calculated from the energy conservation equation written for the 

mixing chamber. The energy conservation equation for node 3 can be written as:

x

H2O m13Cp,H2O T1 + xO2

m

13Cp,O2

T

1 + m23 Cp,O2

T2 = m

34Cp,mixT3 . (71)

Equation (71) can be rearranged to find T3:

T3 =

xH2O m13Cp,H2OT1 + xO2

m13 Cp,O2

T1 + m 23 Cp,O2T2

m 34 cp,mix

 

=(0.9)(1.16 lb/s)(0.5696 Btu/lb-R)(1,960R) + (0.1)(1.16 lb/s)(1.2512 Btu/lb-R)(1,960R)

(4.26 lb/s)(0.3757 Btu/lb-R)

 + (3.1 lb/s)(0.2324 Btu/lb-R)(540R)

(4.26 lb/s)(0.3757 Btu/lb-R) = 1149.1 °R or 689.5 °F . (72)

Equation (72) calculates the temperature to be 689.5 °F, which compares well with the GFSSP 

prediction of 688.6 °F.

The mass concentration of the species can be calculated from the species conservation 

equation for node 3. The concentration of water vapor and oxygen can be expressed as:

X3,H2O =

X1,H2O m13

m 34

= (0.9)(1.16 lb/s)

4.26 lb/s = 0.2451 (73)

and

X3,O2 =

X1,O2

m13 + m 23

m 34

= (0.1) (1.16 lb/s) + 3.1 lb/s

4.26 lb/s = 0.7549 . (74)

The concentration of water vapor and oxygen from the above equations are 0.2451 and 

0.7549, respectively, which also compare well with the GFSSP predictions.

5.5 Example 5—Simulation of a Flow System Involving a Heat Exchanger

5.5.1 Problem Considered

In dealing with fluid system analysis, engineers often encounter systems that contain a heat 

exchanger. It is important that the thermal behavior of a heat exchanger is correctly accounted for 

in any system simulation. Otherwise, temperature discrepancies in the fluid property calculations 

will result in inaccurate system characteristics being predicted. The following example demonstrates 
118

GFSSP’s ability to accurately predict fluid temperatures in a heat exchanger system using effectiveness calculations.

GFSSP has the ability to calculate temperatures downstream of a heat exchanger for three 

different cases. For the first case, a known heat exchanger effectiveness is used by GFSSP to calculate the flow temperatures downstream from the heat exchanger. The second case involves requiring 

GFSSP to calculate the effectiveness of the heat exchanger using the counterflow heat exchanger 

equations. That calculated effectiveness is then used to calculate the heat exchanger downstream 

temperatures. For the third case, the heat exchanger effectiveness is calculated using the parallel 

flow heat exchanger equations and the heat exchanger downstream temperatures are calculated 

from that effectiveness. 

A simple counterflow heat exchanger system configuration, as shown in figure 87, was chosen for this example. As shown in figure 87, counterflow occurs when the hot branch of the heat 

exchanger has flow that is propagating in a direction opposite to the cold branch. This counterflow 

heat exchanger configuration consists of hot water, at 50 psi and 100 °F, flowing through 10 in of 

0.25-in inner diameter pipe, through a 10-in-long heat exchanger and out through another 10-inlong section of 0.25-in inner diameter pipe at 25 psi. Also, cold water, at 50 psi and 60 °F, flows 

through a 10-in section of 0.5-in inner diameter pipe, through the heat exchanger, and out through 

another 10-in section of 0.5-in inner diameter pipe at 25 psi. All of the pipes are assumed to have 

an absolute roughness of zero. At the conditions described above, Cp=0.9978 Btu/lbm-R is calculated for the hot water and Cp=1.0014 Btu/lbm-R is calculated for the cold water.

5.5.2 GFSSP Model

A GFSSP model consisting of eight nodes and six branches can represent the counterflow 

heat exchanger system shown in figure 87. This model is shown in figure 88(a). Nodes 1, 4, 5, and 

8 are boundary nodes. Nodes 1 and 5, the inlet boundary nodes for the hot and cold flow, respectively, both have a pressure of 50 psi. In addition, the boundary temperature at Node 1 is 100 °F 

while the boundary temperature at Node 5 is 60 °F. Nodes 4 and 8, which are the downstream 

boundary nodes for the hot and cold flow, both have a boundary pressure of 25 psi. Downstream 

boundary temperatures are not used in GFSSP calculations so ‘dummy’ temperature values of 

80 and 70 °F, respectively, are used for the hot and cold flow downstream boundary nodes. 

Nodes 2 and 6 are internal nodes that represent the entrances to the heat exchanger for the hot and 

cold flow, respectively. In the same manner, nodes 3 and 7 are internal nodes that represent the hot 

and cold flow heat exchanger exits. Branches 12, 34, 56, and 78 represent the pipes leading into and 

out of the heat exchanger for both the hot and cold flows. Finally, branches 23 and 67 represent the 

hot and cold sides of the heat exchanger. Figure 88(b) shows how this model appears in VTASC.

5.5.2.1 Branches 12, 34, 56, and 78 (Pipe Lines). Option 1 was used to represent each of the 

pipe sections in the heat exchanger model. The user is required to provide the length, inner diameter, and relative roughness factor (ε/D) for this branch resistance option.
119

Pipe

L = 10 in

D = 0.25 in

Water

P = 50 psi

T = 100 °F

Water

P = 50 psi

T = 60 °F

P = 25 psi

P = 25 psi

Heat Exchanger

Pipe

L = 10 in

D = 0.5 in

Pipe

L = 10 in

D = 0.5 in

Pipe

L = 10 in

D = 0.25 in

Figure 87. Flow system schematic of a heat exchanger (Example 5).

1 2 4

1 12 2 4

12

3

23 34

8 7 5

78

6

67 56

23 3 34

8 78 7 67 6 56 5

Water

Pi = 50 psi

Ti = 100 °F

Water

Pi = 50 psi

Ti = 60 °F

Heat Exchanger

Hot Branch

L = 10 in

D = 0.25 in

L = 10 in

D = 0.25 in

L = 10 in

D = 0.25 in

Heat Exchanger

Cold Branch

L = 10 in

D = 0.50 in

L = 10 in

D = 0.50 in

L = 10 in

D = 0.50 in

Pe = 25 psi

Pe = 25 psi

(a)

(b)

Figure 88. GFSSP model of the heat exchanger: (a) Detailed schematic 

and (b) VTASC model.

5.5.2.2 Branches 23 and 67 (Heat Exchanger). Option 1 was also used to simulate the two 

heat exchanger branches. In addition to providing the length, inner diameter, and ε/D for the two 

branches, the user must designate additional information in the Heat Exchanger dialog window 

shown in figure 89. First, the user must add a heat exchanger using the Add button. Next, the user 

must input which branch represents the hot flow side of the exchanger and which branch represents 

the cold flow side of the exchanger. Also, the user may either designate an effectiveness between 

zero and 1 or enter an effectiveness >1 and designate the type of heat exchanger in the system and 
120

Figure 89. VTASC Heat Exchanger dialog.

a value for UA. UA is the product of the overall conductance for heat transfer and the surface area 

on which that conductance is based. 

5.5.3 Results

The example 5 GFSSP input and output data files (ex5.dat and ex5.out) are included in 

appendix H (on CD inside back cover). The output file includes all of the input data, the trial solution for the internal nodes, and the final model solution.

The first area of interest is the counterflow heat exchanger effectiveness calculations. During the GFSSP preprocessor input, a value for UA and an effectiveness are input to define the heat 

exchanger’s characteristics. For these verifications, a value for UA was assumed with the following 

process. It was assumed that the heat exchanger had an effectiveness of ε =0.7. Then, a GFSSP 

model was run with that effectiveness to obtain the mass flow rates. The mass flow rate for the cold 

branch was calculated to be 5.41 lbm/s and the flow rate for the hot branch was calculated to be 

0.885 lbm/s. Equations (75) and (76) were used to calculate the hot and cold fluid capacity rates:24

Ch = (m Cp )hot branch = (0.885 lbm/s)(0.9978 Btu/lbm –oR) = 0.883 Btu/s-oR (75)

and

Cc = (m Cp )cold branch = (5.41 lbm/s)(1.0014 Btu/lbm –oR) = 5.418 Btu/s–oR . (76)

Based on the previously calculated values, Cmax = Cc and Cmin = Ch. For the counterflow 

heat exchanger case, it was assumed that Cmax, Cmin and ε would remain the same as the values above. Next, a counterflow exchanger performance table was used, along with the previously 

calculated values, to estimate the number of heat transfer units, Ntu=1.25, for the counterflow 

exchanger.24 Then, equation (77) was used to calculate UA in equation (78):24
121

Ntu = UA

Cmin

(77)

and

UA = NtuCmin = (1.25)(0.883 Btu/s −oR) = 1.10375 Btu/s −oR) . (78)

A value >1 was used for the effectiveness, which instructed GFSSP to calculate the effectiveness instead of employing a user input value. For this case, the counterflow heat exchanger option 

was chosen. Figure 90 shows the predicted temperatures and mass flow rates from this model.

1 2 4

12 23 34

8 5

78

6

67 56

M = 0.885 1×10 lb/s

M = 0.541 1×101 lb/s

T = 0.7180 1×102 °F

T = 0.6451 1×102 °F

3

7

Figure 90. Temperature and flow rate predictions in heat exchanger.

The mass flow rates predicted by this model (0.885 lbm/s at branch 12 and 5.41 lbm/s at 

branch 56) were the same as predicted for the counterflow exchanger that was used for the assumed 

case. Therefore, Cc and Ch remained the same. The heat exchanger effectiveness is not included 

in the data written to the GFSSP output file so it was necessary to calculate the effectiveness that 

GFSSP used. GFSSP’s effectiveness was calculated using equation (79):24

ε = Ch (Th,in −Th,out )

Cmin(Th,in −Tc,in )

= Cc (Tc,out −Tc,in )

Cmin(Th,in −Tc,in )

ε = (0.883 Btu/s–oR)(560 oR − 531.8 oR)

(0.883 Btu/s–oR)(560 oR − 520 oR) = 0.705 . (79)

For comparison, a hand-calculated counterflow effectiveness was determined using 

equation (80):24

ε = 1− e

−Ntu 1−Cmin /Cmax ( )

1− Cmin ( / Cmax )e

−Ntu 1−Cmin /Cmax ( )
122

ε = 1− e

−1.25(1− 0.883 Btu/s–oR

5.418 Btu/s–oR

⎛

⎝

⎜ ⎞

⎠

⎟

1− 0.883 Btu/s–oR

5.418 Btu/s–oR

⎛

⎝

⎜ ⎞

⎠

⎟ e

−1.25 1− 0.883 Btu/s–oR

5.418 Btu/s–oR

⎛

⎝

⎜ ⎞

⎠

⎟

= 0.688 . (80)

Good agreement can be seen in a comparison between the hand-calculated value and the GFSSP 

value. 

The second area of interest is the accuracy of GFSSP’s temperature predictions at the nodes 

downstream of the hot and cold heat exchanger branches. Equation (79) was manipulated to come 

up with equations (81) and (82) which were used to hand calculate Th,out and Tc,out, respectively:24

Th,out = Th,in −

Cmin Th,in −T ( c,in ) ε

Ch

Th,out = 560 oR − (0.883 Btu/s–oR)(560 oR–520 oR) (0.688)

(0.883 Btu/s–oR) = 532.48 oR = 72.48 oF  . (81)

and

Tc,out = Tc,in +

Cmin Th,in −T ( c,in ) ε

Cc

  .

Tc,out = 520 oR − (0.883 Btu/s–oR)(560 oR–520 oR) (0.688)

(5.418 Btu/s–oR) = 524.49 oR = 64.49 oF  . (82)

Comparing the hand-calculated values with GFSSP’s temperature results of 71.8 °F 

at node 3 and 64.51 °F at node 7 in figure 90 shows very good agreement, verifying the heat 

exchanger temperature calculation process used by the GFSSP code.

5.6 Example 6—Radial Flow on a Rotating Radial Disk

5.6.1 Problem Considered

This example illustrates the rotational effect (centrifugal force contribution) capability 

of GFSSP by modeling the flow of water through a closed impeller.27 The impeller has a diameter 

of 11 in, it uses water as the operating fluid, and is running at 5,000 rpm. The ‘slip’ of the fluid 

is described by the rotational K-factor (Krotation). Krotation is defined as the ratio of the mean 

circumferential fluid speed divided by the impeller speed: (Krotation=uθ/rw). (Higher Krotation
123

factors translate to a higher pressure rise for radially outward flow.) A Krotation for each side 

of the impeller has been proposed and it is the purpose of this example to validate these Krotation

correlations. The proposed correlations are: Kfront face=0.8455– 0.1403(r– ri /ro –ri) for the front 

face; Kback face=0.8857–0.1762(r –ri /ro –ri) for the back face. It is desired to compare the results 

with experimental data. The impeller is schematically shown in figure 91. For this example, the 

effects of friction will be neglected for the rotating branches.

Figure 91. Flow schematic of a rotating radial disk (example 6).

5.6.2 GFSSP Model

The GFSSP model circuit is shown in figure 92(a). All branches are modeled using option 2. 

In the model, branches 23, 34, 45, 56, 67, 89, 910, 1011, and 1112 are rotating at 5,000 rpm. The 

inlet and outlet radii are defined in the preprocessor for each of the rotating branches. The area 

of each of the radial branches is calculated as the average cross-sectional area for each branch. 

Figure 92(b) shows how this model appears in VTASC.

5.6.3 Results

The example 6 GFSSP input and output data files (ex6.dat and ex6.out) are included in 

appendix I (on CD inside back cover).

The pressure distribution predicted by GFSSP for the front and back faces of the impeller 

is shown in figure 93. As is seen in figure 93, the model results show excellent agreement with the 

experimental data.
124

78

67

7

r = 5.5 in

r = 5.375 in

r = 4.6875 in

r = 3.625 in

Rotating 

Branches

Rotating 

Branches

r = 2.65 in

r = 2.25 in

r

CL

r = 1.25 in

56

45

34

89

910

1011

1112

1213

23

12

8

6 9

5 10

4

11

3

12 13

1 2

1

13

12

23

34

45

56

67

78

89

910

1011

1112

1213

R 2

R

R

3

R

4

R

5

R

6

7

R

12

R

11

R

10

R

9

8

R

R

(a)

(b)

Figure 92. GFSSP model of the rotating radial disk: (a) Detailed schematic 

and (b) VTASC model.
125

(R/Ro)

2

∆p(ρ(ω/Ro)

2)

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

0.25

0.2

0.15

0.1

0.05

0

–0.05

GFSSP model using new correlation method-front face

GFSSP model using new correlation method-back face

Pratt and Whitney water experimental data-front face

Pratt and Whitney water experimental data-back face

Figure 93. Comparison of GFSSP model results with experimental data.

5.7 Example 7—Flow in a Long-Bearing Squeeze Film Damper

5.7.1 Problem Considered

Squeeze film dampers are used in turbomachinery to dampen out unstable behavior. The 

damper is installed at the bearing supports of a rotor-stator system on the outer race of a rollingelement bearing. The squeeze film damper consists of inner and outer elements separated by fluid 

(usually an oil). The inner element is mounted to the outer race of the rolling-element bearing, and 

the outer element is mounted to the bearing support. The arrangement is similar to a journal bearing except that the inner damper element does not rotate; it only translates. In order to calculate the 

effect of the squeeze film damper on the system, the forces generated by the squeeze film damper in 

the radial and tangential directions must be estimated. The forces are estimated by integrating the 

pressure distribution of the fluid in the damper. The difficulty for the designer/analyst is the estimation of the pressure distribution. The following example problem demonstrates how GFSSP can be 

used to predict this pressure distribution.

The squeeze film damper considered for this example is shown schematically in figures 94 

and 95. Since the damper has sealed ends, the axial flow is neglected. The diameter (d) of the bearing is 5 in, the width (w) of the bearing is 0.94 in, the clearance (c) is 0.0625 in, and the ratio of the 

dynamic eccentricity (ε, radius of orbit to the clearance) is 0.82. The fluid density (ρ) is 57.806 lbm/

ft3. The fluid viscosity (μ) is 5.932×10–3 lbm/(ft•s). The rotational speed is set at 1,770 rpm 

(ω=185.354 rad/s). GFSSP will be used to determine the pressure distribution around the damper 

and the results will be compared with experimental data.28
126

Bearing

Rotar Shaft

Squeeze Film Region

Inner Segment

Outer Segment

Figure 94. Squeeze film damper schematic (example 7, view 1).

Fluid

Inner Segment of Squeeze Film Damper

Outer Segment of Squeeze Film Damper

Figure 95. Squeeze film damper schematic (example 7, view 2).

5.7.2 GFSSP Model

A GFSSP model consisting of 20 nodes and 19 branches can approximate the fluid contained within the squeeze film damper system shown in figures 94 and 95. The fluid to be used is 

not contained in the standard library of fluids and is assumed to be incompressible; therefore, the 

constant density feature of GFSSP must be used. In order to model the squeeze film damper, the 

damper will be ‘unwrapped.’ Figure 96 shows the unwrapping of the damper and the discretization 

of the flow region. The GFSSP model is shown in figure 97. As is shown in figure 97, nodes 1 and 

20 are the boundary nodes. The branches will use branch resistance option 3 - Non Circular Duct, 
127

Housing

Fluid

Bearing

e ω

x

y

y

x

y

x

Figure 96. Unwrapping and discretization of squeeze film damper.

1 20

12

2

23

3

34

4

45

5

56

6

67

7

78

8

89

9

910

10

1011

11

1112

12

1213

13

1314

14

1415

15

1516

16

1617

17

1718

18

1819

19

1920

Figure 97. GFSSP model of squeeze film damper.

Sub-option 1 - Rectangular Duct. The heights of the branches are given in table 20. The motion of 

the inner element will be simulated using the moving boundary option in GFSSP. (Only the motion 

of the damper normal to the inner element is modeled in this technique.) The velocity of the moving boundary is given in table 21. Setting the boundary nodes at the same pressure will simulate the 

periodic behavior of the damper. The boundary pressure is set at 0.0 psi.

5.7.3 Results

The Example 7 GFSSP input and output data files (ex7.dat and ex7.out) are included in 

appendix J (on CD inside back cover).

The pressure distribution predicted by GFSSP is shown in figure 98. The plot of pressure 

versus angle (i.e., node position) in figure 98 shows that the pressure is symmetric about the boundary pressure of zero psi. The model results are compared to experimental results in figure 99. 
128

Table 20. Branch dimensions of squeeze film damper.

Branch 

Number

Width 

(in)

Length 

(in)

Height 

(in)

12

23

34

45

56

67

78

89

910

1,011

1,112

1,213

1,314

1,415

1,516

1,617

1,718

1,819

1,920

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.94

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.82673

0.012578

0.017987

0.028221

0.042169

0.058320

0.074925

0.090183

0.102441

0.110370

0.113113

0.110370

0.102441

0.090183

0.074925

0.058320

0.042169

0.028221

0.017987

0.012578

Table 21. Moving boundary information of squeeze film damper.

Node

Normal Area 

(in2)

Velocity 

(ft/s)

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.777126

0.256180

0.484598

0.660503

0.764832

0.786280

0.722522

0.580468

0.375510

0.129861

–0.129861

–0.375510

–0.580468

–0.722522

–0.786280

–0.764832

–0.660503

–0.484598

–0.256180
129

In figure 99, the pressure is normalized with a characteristic pressure (CpRe) and the angle has 

been converted to compare to a dimensionless time (wt) for comparison with the experimental data. 

As seen in figure 99, the pressure profile of the GFSSP model compares favorably with the experimental results in shape and magnitude, the only major difference between the two results is a phase 

shift.

15

10

5

0

–5

–10

–15

Pressure (psi)

0 1 2 3 4 5 6 7

Angle (Radians)

Figure 98. Predicted circumferential pressure distributions in the squeeze 

film damper.

0.8

0.6

0.4

0.2

0

–0.2

–0.4

–0.6

–0.8

–1

p/CpRe

0 1 2 3 4 5 6 7

Initial 20 Node GFSSP Model

Experimental Data

ω t

Figure 99. Comparison of GFSSP model results with experimental data 

for squeeze film damper.
130

5.8 Example 8—Simulation of the Blowdown of a Pressurized Tank

5.8.1 Problem Considered

In the previous examples, the simulation of steady state flow in a given flow circuit was considered. In this example, the capabilities of the unsteady flow formulation of GFSSP to simulate 

the process of blowing down a pressurized tank will be employed.

Consider a tank with an internal volume of 10 ft3, containing air at a pressure and temperature of 100 psia and 80 oF, respectively. Air is discharged into the atmosphere through an orifice 

of 0.1-in diameter for a period of 200 s. GFSSP will be used to determine the pressure, mass flow 

rate, and temperature history of the isentropic blowdown process. These predicted values are then 

compared with the analytical solution. 

5.8.2 GFSSP Model

The physical schematic for example 8 is shown in figure 100(a) and a schematic of the corresponding GFSSP model is shown in figure 100(b). Figure 100(c) shows how the model looks 

in VTASC. The venting process can be modeled with two nodes and one branch. Node 1 is an 

internal node that represents the tank. For the unsteady formulation, the node volume and the 

initial conditions must be supplied for each internal node and a history file must be supplied for 

each boundary node. The history file contains the pressures, temperatures, and concentrations at 

discrete times. At a minimum, this file should include values at the process start time and at some 

time corresponding to the expected process stop time. Additional times can be included to account 

for nonlinear variation in the values if required. The code interpolates in the history file data to 

determine the values for a particular instant. Shown below is a listing of EX8HS2.DAT, which is 

the history file of example 8, used to provide the boundary conditions for node 2. The file listing 

has been annotated to explain the meaning of the entries. In the actual file, column headings and 

descriptions must not appear.

EX8HS2.DAT

2 - Number of data points

tau(sec) p(psia) T (oF) Concentration

0 14.700 80.00 1.00

1000 14.700 80.00 1.00

In addition to supplying the internal node volumes and history data files, the time step 

(DTAU), start time (TAUF), stop time (TAUL), and print interval (NPSTEP) must also be 

included within the model input data file when creating an unsteady flow (STEADY = .False.) 

model.

The initial pressure within the tank (node 1) was 100 psia. Resistance Option 22 was used 

for branch 12 with a flow coefficient of 1. Air is modeled using the ideal gas option that is available 

in VTASC. 
131

Fluid: Air

Pt=0 = 100 psia

Tt=0 = 80 °F

V = 17,280 in3

Ex8Hs2.dat

P = 14.7 psia

T = 80 °F

CL = 1 

A = 0.00785 in2

 i m = 0.0177 lbm/s

Atmosphere

p = 14.7 psia

d = 0.1 in

V = 10 ft3

Initial Conditions

Tank

pi

 = 100 psia

Ti

 = 80 °F

1 2

12

2

12

1 C O

(a)

(b)

(c)

Figure 100. Venting air tank schematics: (a) Physical schematic, (b) detailed 

model schematic, and (c) VTASC model.

As an interesting note, this example could be used to model an isothermal process by using 

the SORCEQ user subroutine. The temperature of the fluid remains constant in an isothermal process. In this example it is presumed that initially the air and tank wall are at the same temperature. 

During blowdown, the air temperature tends to drop. With heat transfer from the wall, temperature 

drop would be less compared to an isentropic process. For an isothermal process, there will be no 

change in temperature. This particular situation (isothermal) can be modeled by setting an infinite 

heat transfer coefficient between the wall and fluid in SORCEQ.

5.8.3 Results

The input and output files of this example are included in appendix K (see CD inside back 

cover) as ex8.dat and ex8.out. It may be noted that for each time step, solutions for each node and 

branch are printed in the output file.

5.8.3.1 Analytical Solution. The differential equation governing an isentropic blowdown 

process can be written as:
132

p

i p

⎛

⎝

⎜

⎞

⎠

⎟

(1−3γ )/2γ d p / pi ( )

dτ

 = γ A

i ρV γ c g i p i ρ

2

γ +1

⎛

⎝

⎜

⎞

⎠

⎟

(γ +1)/2(γ −1)

. (83)

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

τ = 0, p

i p = 1 . (84)

The analytical solution for p/pi is given by Moody as:29

p

i p = 1+ γ −1

2

⎛

⎝

⎜ ⎞

⎠

⎟ 2

γ + 1

⎛

⎝

⎜

⎞

⎠

⎟

(γ +1)/2(γ −1) γ c g i p

ρi

Aτ

V

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

−2γ /(γ −1)

. (85)

The analytical and GFSSP solutions are compared in figure 101. The figure shows a comparison between the GFSSP solution and the analytical solution of pressures. The difference in 

pressures is also shown plotted for three different time steps (1, 0.1, and 0.01 s). The discrepancies 

between analytical and numerical solutions are found to diminish with reduction in time step. This 

observation is in conformity with expectations.

5.9 Example 9 —A Reciprocating Piston-Cylinder

5.9.1 Problem Considered

This example further illustrates GFSSP’s capability to model complex unsteady flow. Figure 102 shows the piston-cylinder configuration considered by this example problem. The cylinder 

has a diameter of 3 in. Within the cylinder is nitrogen gas, sealed in by a piston moving at a rotational speed of 1,200 rpm and a stroke of 3 in. GFSSP will be used to predict the pressure and 

temperature within the system and the results will be compared with the isentropic solution.

5.9.2 GFSSP Model

In order to model this configuration, a coordinate transformation is utilized. In this new 

coordinate system, the endplate of the cylinder is modeled as another piston and the origin of the 

coordinate system is at the midpoint between the two ‘pistons.’ Figure 103 demonstrates the modified piston-cylinder arrangement.

The GFSSP model of the piston-cylinder arrangement consists of two internal nodes 

and one branch. (Note: the model does not have any boundary nodes.) Figure 104(a) shows the 

GFSSP piston-cylinder model. In order to model the motion of the piston, two special options are 

utilized—the moving boundary option and the variable geometry option. The moving boundary 

option is required to adequately model the work input by the motion of the pistons. The variable 

geometry option is required to model the variation of the geometry of branch 12. The initial condition of the nitrogen in the cylinder is 14.7 psia and 75 oF. Figure 104(b) shows how the model looks 

in VTASC.
133

120

100

80

60

40

20

0

0.1

0.09

0.08

0.07

0.06

0.05

0.04

0.03

0.02

0.01

0

Pressure (psia)

∆

P (psid)

Time (s)

0 20 40 60 80 100 120 140 160 180 200

P (GFSSP) (psi)

P (Moody) (psi)

(GFSSP Minus Moody, 1 s Time Step)

(GFSSP Minus Moody, 0.1 s Time Step)

(GFSSP Minus Moody, 0.01 s Time Step)

Figure 101. Comparison of the predicted pressure history by GFSSP 

and the analytical solution.

d = 3 in

L(t) = 4 + 3 cos ( t) in

x

Figure 102. Piston-cylinder configuration.
134

d = 3 in

L(t) = 4 + 3 cos ( t) in

x

Figure 103. Coordinate transformed piston-cylinder configuration.

1

1

12

12

2

2

L12 = 4 + 3 cos ( t) in

d12 = 3 in

p1i = p2i = 14.7 psia

Moving Boundary Term 

Associated With Nodes 1 and 2

Boundary Node

Internal Node

Branch (a)

(b)

Figure 104. GFSSP model of the piston-cylinder: (a) Detailed model schematic 

and (b) VTASC model.

5.9.3 Results

Appendix L (on CD inside back cover) contains the input, variable geometry history, 

and output files for this example. The results of the study are compared to an analytical solution 

(for constant ratio of specific heat, γ ). Equations (86) and (87) are used to obtain the analytical 

solution assuming an isentropic process: 

T(t) = T0 ( + 459.6) ρ(t)

ρ0

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

γ −1 ⎧

⎨

⎪

⎩

⎪

⎫

⎬

⎪

⎭

⎪

− 459.6 (86)

and
135

p(t) = p0

ρ(t)

ρ0

⎡

⎣

⎢

⎢

⎤

⎦

⎥

⎥

γ

, (87)

where T0, p0, and ρ0 are temperature, pressure, and density at time equal zero.

Figures 105 and 106 compare the results of the GFSSP piston-cylinder model with the 

analytical solution from equations (86) and (87). As these figures illustrate, the GFSSP model 

compares favorably to the analytical solution. It should be noted that the isentropic solution uses 

a ̨constant ratio of specific heats (g), whereas GFSSP accounts for the variation of specific heat 

ratios with changes in temperature and pressure.

800

700

600

500

400

300

200

100

0

Temperature (°F)

Time (s)

0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05

Isentropic Solution for =1.4

GFSSP Prediction

Figure 105. Predicted temperature history of piston-cylinder model.

5.10 Example 10 —Power Balancing of a Turbopump Assembly

5.10.1 Problem Considered

This example, Example 10, illustrates the modeling of the mechanical coupling between two 

flow components. In the turbopump assembly shown in figure 107, a coaxial shaft mechanically 

connects the pump and turbine. The power required by the pump must be transmitted from the 

turbine in order for the system to be in balance. The purpose of this example is to demonstrate this 

power balancing for a turbopump when used in a gas turbine cycle. The physical plausibility of the 

predicted results was demonstrated by performing parametric studies on shaft speed.
136

250

200

150

100

50

0

Temperature (°F)

Time (s)

0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04 0.045 0.05

Isentropic Solution for =1.4

GFSSP Prediction

Figure 106. Predicted pressure history of piston-cylinder model.

Gas Generator

Pump Turbine

m

w

Figure 107. Simplified turbopump assembly.

5.10.2 GFSSP Model

A model of the turbopump portion of a flow circuit is shown in figure 108(a). This model 

consists of an inlet from a hydrogen tank, a turbopump assembly (pump, turbine, and a connection between them (shaft)), two heat exchangers, a bypass dump outlet, and an outlet to the power 

turbine. The first of the heat exchangers (denoted on fig. 108(a) as the Regenerator) is used to heat 

a small portion of the main LH2 flow by using the ‘hotter’ hydrogen exiting the turbine, while the 
137

1

17

20

2

3

4 5 7 8

9

13 14 15 16

12 11 10

19 18

6

To Power 

LH2 Tank Turbine

12

45

1920

1011

1314

1112

46 68

78

89

1516 

Line

910 

Cold Line

1415 Hot Line

1213 

Turbine

1819 

Hot Line

1618

Q

1617

To 

Bypass

23 

Pump

57 Cold Line

Regenerator

Heat Exchanger

34 

Turbopump 

Assembly

(a)

(b)

20

17

19

12

12

11 9

5 7 8

6

4

13

1

2

3 14 15

1920

1112

23 1213

34

45 57 78

1314 1415 1516

89

46 68

1011 910

1819

1618

1617

18 Cv

Cv

P P

16

Cv

Cv

Cv 10

Figure 108. A flow circuit containing turbopump assembly: (a) Detailed model 

schematic and (b) VTASC model.

remainder of the LH2 flow bypasses this heat exchanger. The second heat exchanger is used to boil 

and superheat the hydrogen by means of external heat addition. The shaft speed for this model is 

set in the input file to 80,000 rpm. Figure 108(b) shows how this model looks in VTASC.
138

This model uses the following options:

• Branch resistance options:

(1) Pipe flow (option 1).

(2) Pump with pump efficiency (option 15).

(3) Valve with a given Cv (option 16).

• Special options:

(1) Heat exchangers.

(2) Turbopump assembly.

Each pipe flow branch has a length of 100 in, an inside diameter of 0.3927 in, and an 

absolute roughness of 0.00098175 in. Branches 89, 1112, and 1618 each have Cv=3.554 and 

A=0.19635 in2. Branch 12 has Cv=2.877 with A=0.19635 in2 while branch 1617 has Cv=0.00354 

with A=0.01 in2. The branch options chosen to represent the turbine (branch 1213) and the pump 

(branch 23) have no bearing on the model calculations except for the flow areas that are provided. 

For this case, option 15 was used for both branches with arbitrary inputs for pump horsepower 

and efficiency. The flow areas are 0.12112 in2 for branch 23 and 0.019635 in2 for Branch 1213. The 

turbopump characteristics are defined in the turbopump dialog shown in figure 109. The pump 

characteristic file is also shown below with annotations to explain the meanings of each value. The 

two heat exchanger dialogs are shown in figure 110(a) and (b).

EX11PMP23.DAT

18 – Number of lines of data

Flowrate/Speed Head/Speed2 Torque/(Density*Speed2)

0.000 8.680E-06 0.000

3.035E-05 8.971E-06 8.8724E-10 

6.071E-05 9.190E-06 9.7065E-10

9.106E-05 9.341E-06 1.0804E-09

1.214E-04 9.436E-06 1.2166E-09

1.518E-04 9.486E-06 1.3393E-09

1.821E-04 9.486E-06 1.4570E-09

2.125E-04 9.445E-06 1.5644E-09

2.428E-04 9.372E-06 1.6733E-09

2.732E-04 9.263E-06 1.7872E-09

3.035E-04 9.117E-06 1.9105E-09

3.339E-04 8.935E-06 2.0558E-09

3.643E-04 8.753E-06 2.2161E-09

3.718E-04 8.689E-06 2.2698E-09

3.749E-04 8.625E-06 2.2869E-09

3.794E-04 8.479E-06 2.3215E-09

3.807E-04 8.388E-06 2.3281E-09

3.810E-04 0.000E+00 0.000
139

Figure 109. Example 10 turbopump dialog.

(a) 

(b) 

Figure 110. Example 10 heat exchanger dialog: (a) Heat exchanger 1 and (b) heat exchanger 2.

5.10.3 Results

Appendix M (see CD inside back cover) contains the input, pump characteristics, and 

output files for this example (ex10.dat, ex10pmp23.dat, and ex10.out). The results of the study are 

illustrated in figure 111. 
140

20

17

19

12

12

11 9

5 7 8

6

4

13

1

2

3 14 15

1920

1112

23 1213

34

45 57 78

1314 1415 1516

89

46 68

1011 910

1819

1618

1617

18 Cv

Cv

P P

16

Cv

Cv

Cv 10

M = 0.224 1×10 lb/s

M = 0.145 1×10 lb/s

M = 0.787 1×10–1 lb/s M = 0.773 1×10–3 lb/s

M = 0.224 1×10 lb/s

M = 0.223 1×10 lb/s

P = 0.6000 1×102 psia

T = 0.4190 1×103 °F

P = 0.1780 1×104 psia

T = 0.4072 1×103 °F

P = 0.1777 1×104 psia

T = 0.1462 1×103 °F

P = 0.1645 1×104 psia

T = 0.1478 1×103 °F

P = 0.2539 1×103 psia

T = 0.2250 1×103 °F

P = 0.5766 1×103 psia

T = 0.2489 1×103 °F

P = 0.1765 1×104 psia

T = 0.3066 1×103 °F

P = 0.6272 1×103 psia

T = 0.3550 1×101 °F

W = 64.9 hp

 = 200 Btu/s Q

Figure 111. GFSSP RCS model results.

5.10.4 Parametric Study

In order to verify this complicated model, a parametric study on the shaft speed of the turbopump was conducted. Figures 112–114 illustrate the results of this model. Figure 112 illustrates 

the pressure differential across the turbopump for both the pump and the turbine as a function 

of the shaft speed. Figure 113 illustrates the hydrogen mass flow rate through the turbopump as 

a function of the shaft speed. Figure 114 illustrates the torque and horsepower transmitted in the 

turbopump as a function of the shaft speed. As each of these figures illustrates, a functional relationship is identifiable for each predicted variable as a function of shaft speed.
141

1,400

1,200

1,000

800

600

400

200

0

Pressure Differential (psi)

0 10,000 20,000 30,000 40,000 50,000 60,000 70,000

Shaft Speed (rpm)

Pressure Rise Across Pump

Pressure Drop Across Turbine

Figure 112. Parametric study results: turbopump pressure differential.

0.16

0.14

0.12

0.1

0.08

0.06

0.04

0.02

0

Mass Flow Rate (lbm/s)

0 10,000 20,000 30,000 40,000 50,000 60,000 70,000

Shaft Speed (rpm)

Figure 113. Parametric study results: turbopump hydrogen mass flow rate.
142

40

35

30

25

20

15

10

5

0

0 10,000 20,000 30,000 40,000 50,000 60,000 70,000

Shaft Speed (rpm)

Torque (in-lb)

Horsepower

Figure 114. Parametric study results: turbopump torque and horsepower.

5.11 Example 11—Steady State and Transient Conduction 

Through a Circular Rod, With Convection

5.11.1 Problem Considered

In the previous examples, focus was on GFSSP’s fluid modeling capabilities. In this example, 

GFSSP’s ability to model applications with conjugate heat transfer is introduced. The verification 

and validation of GFSSP’s conjugate heat transfer capability was performed by comparing with the 

known solution of a simple conduction-convection problem.4 The heat transfer in a homogenous 

circular rod between two walls was considered (fig. 115). The two walls are held at temperatures of 

32 oF and 212 oF, respectively. The 0.167-ft-diameter rod is 2 ft in length and is initially at a temperature of 70 oF. The heat transfer coefficient between the rod and the ambient air is 1.14 Btu/ft2-

hr-R and the thermal conductivity of the rod is 9.4 Btu/ft-hr-R. 

5.11.2 GFSSP Model

Figure 116(a) shows a GFSSP schematic of the system described in figure 115. The circular 

rod is represented by eight solid nodes and seven solid-solid conductors. Even though all of the 

material properties are not used for a steady state model, GFSSP still requires that placeholder 

values be input at each solid node. The thermal conductivity and specific heat temperature history 

files, which are shown below in an annotated form, are based on a user-defined material.
143

USER1K.PRP

2 – Number of lines of data

Temperature ® Thermal Conductivity (Btu/ft-sec-R)

0 0.002611

1000 0.002611

USER1CP.PRP

2 – Number of lines of data

Temperature ® Specific Heat (Btu/lbm-R)

0 0.1981

1000 0.1981

Two ambient nodes are used to model the hot and cold walls, and their interaction with the 

rod is modeled using two solid-ambient conductors. Because GFSSP is first and foremost a fluid 

analysis code, it is necessary to include a fluid flow path in any GFSSP model that is being developed. Therefore, a dummy flow circuit consisting of two boundary nodes, two internal nodes and 

three pipe flow branches was used to represent the ambient environment of figure 15. The details of 

the flow path were arbitrarily chosen with the sole intent of maintaining a constant temperature of 

70 oF at all points in the flow path to correctly simulate the ambient environment. The two internal 

flow nodes are connected to the solid rod by eight solid-fluid conductors that represent the heat 

transfer between the ambient and the rod in figure 115. Figure 116(b) shows how this model looks 

in VTASC.

Tamb = 70 °F

D = 0.167 ft

T = 32 °F T = 212 °F

L = 2 ft

Figure 115. Schematic of circular rod connected to walls at different temperatures.
144

11

2

2 3 4 5

3 4 5 6 7 8 9

12 13 14

125

1112 1213 1314

122

123

1 12

1 10

23 34 45 56 67 78 89 910

124

136 139

137 138

Water

P = 50 psia

T = 70 °F

P = 45 psia

T = 70 °F

L = 0.1 in

D = 1.73 in

T = 32 °F T = 212 °F

L = 12 in

D = 1.73 in

A = 18.85 in2

hc = 0.000317 Btu/ft2-s-R 

each branch

A = 3.14159 in2

hc = 0.022 Btu/ft2-s-R

Solid Node Properties 

Material: User 1

Ti

 = 70 °F

m = 1 lbm

Cp = 0.1981 Btu/lbm-R

k = 0.002611 Btu/ft-s-R

D = 3 in

A = 3.14159 in2

each branch

L = 12 in

D = 1.73 in

10

11 12 14

1112

12 23 34 45

122 123 124 125

1213 1314

HTC-R

FS

K

SS

K

SS

K

SS

56

K

SS

HTC-R

SA

910

HTC-R

SA

HTC-R

FS

HTC-R

FS

HTC-R

FS

6 7 8 9

67 78 89

136 137 138 139

HTC-R

FS

K

SS

K

SS

K

SS

HTC-R

FS

HTC-R

FS

HTC-R

FS

13

(a)

(b)

Figure 116. GFSSP model of circular rod for example 11: (a) Detailed model schematic 

and (b) VTASC model.

5.11.3 Results

The input and output files, including property files, of Example 11 have been attached 

in appendix N (see CD inside back cover). 

5.11.3.1 Analytical Solution

From the Thermal Analysis Workbook,4 the differential equation of energy transport 

is given by

d2T

dx2 − 4h

Dk

T −T ( ∞ ) = 0 (88)

with boundary conditions of T(0)=32 oF and T(L)=212 oF. The closed form of the solution is given 

by
145

T (x) = Tamb + 4.653e

4h

Dk

x

− 42.65e

− 4h

Dk

x

. (89)

Figure 117 compares GFSSP’s predicted temperature distribution along the rod with the 

closed-form solution. The comparison shows very close agreement between GFSSP and the closedform solution.

250

200

150

100

50

0

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Distance (ft)

Temperature (°F)

Closed Form Solution

GFSSP

Figure 117. Comparison of GFSSP temperature prediction and closed-form solution.

5.12 Example 12—Simulation of Fluid Transient Following Sudden Valve Closure

5.12.1 Problem Considered

This example takes advantage of GFSSP’s capability to model fluid transients. Fluid transients, also known as water hammer, can have a significant impact on the design and operation 

of both spacecraft and launch vehicle propulsion systems. These transients often occur at system 

activation and shut down, and they must be predicted accurately to ensure the structural integrity 

of the propulsion system fluid network.

Consider the system shown in figure 118. Liquid oxygen (LOX) at 500 psia and 200 oR flows 

through a 400-ft-long, 0.25-in inside diameter pipeline at a mass flow rate of 0.1 lbm/s. The corresponding downstream pressure is 450 psia. At time zero, a valve at the end of the pipe begins 

a 100-ms rapid closure. This example discusses how to predict the liquid’s response to the sudden 

valve closure, including the maximum expected surge pressure in the line.
146

5.12.2 GFSSP Model

The GFSSP model of the propellant tank and pipeline schematic of figure 118 is shown 

in figure 119. The system is represented by seven nodes and six branches. Node 1 represents the 

propellant tank as a boundary node. Node 7 represents the downstream pressure as a boundary 

node. Nodes 2 to 6 are internal nodes where pressure, temperature, and density are calculated. 

Branches 12 to 56 represent 80-ft-long, 0.25-in-diameter pipe segments. Branch 67 represents the 

valve as a flow through a restriction with a flow area of 0.0491 in2 and a flow coefficient of 0.6. 

In addition, GFSSP’s unsteady valve open/close option is used to model the rapid valve closure in 

branch 67. The valve open/close characteristics are defined in the valve open/close dialog shown in 

figure 120. The valve open/close characteristic file is also shown below with annotations to explain 

the meanings of each value. GFSSP’s restart option is used to set the initial conditions of the 

model. The model is run as a steady state model, and the solution is saved in files FNDEX12.DAT 

and FBREX12.DAT. This solution is then read in as the initial solution for the transient model 

run.

EX12VLV.DAT

7 – Number of lines of data

Time (sec) Flow Area (in2)

0.00 0.0491

0.02 0.0164

0.04 0.00545

0.06 0.00182

0.08 0.00061

0.1 1.E-16

100 1.E-16

Flow Rate = 0.1 lb/s

Valve Closes in 100 ms

450 psia

400 ft

D = 0.25 in

Liquid Oxygen 

at 500 psia 

and 200 R

Figure 118. Schematic of a propellant tank, pipeline, and valve.

1 7

56 67

5 6 R

45

4

34

3

23

2

12

Figure 119. GFSSP model of a propellant tank, pipeline, and valve.
147

Figure 120. Example 15 Valve Open/Close dialog.

Some consideration should be given to the time step for this model. In order to properly 

model the fluid transient, a time step must be chosen that is small enough to accommodate the 

model discretization. This is done by calculating the Courant number, which is the period of oscillation of one branch divided by the time step (eq. (90)). For the model to properly capture the fluid 

transient phenomena, the Courant number should be greater than unity.

Courant number = 4Lbranch

afluidΔτ

≥ 1 . (90)

The speed of sound (afluid) for LOX is 2,462 ft/s. Choosing a time step of 0.02 s gives 

a Courant number of 6.5, which satisfies the criteria for this model. 

5.12.3 Results

The input and output files including history and restart files of example 15 have been 

attached in appendix O (see CD inside back cover). Figure 121 compares GFSSP’s predicted pressure at node 6 with a Method of Characteristics (MOC) solution. Both solutions compare very well 

in the timing and character of the predicted pressure oscillations. The maximum surge pressure 

predicted by GFSSP is 624 psia compared to a MOC prediction of 636 psia. Additional studies 

have been performed using varying levels of discretization, different fluids, and even liquid-gas 

mixtures.30
148

700

600

500

400

300

200

100

0

Pressure (psia)

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

Time (s)

MOC

GFSSP

Figure 121. Comparison of GFSSP and MOC predicted pressure oscillations.

5.13 Publications

In addition to the examples described in this section, there are large numbers of GFSSP 

applications published in journals and presentations at conferences. A list of publications appears 

in appendix P (see CD inside back cover).
149

APPENDIX A—DERIVATION OF Kf FOR PIPE FLOW

It is assumed that there is a dynamic equilibrium that exists between the friction and the 

pressure forces. Therefore, the momentum conservation equation can be expressed as:

Pu − Pd = Kf m 2 , (91)

where Kf is a function of f, L, D, and r.

For a fully developed pipe flow, the momentum conservation equation can be written as:

τ π DL = (Pu − Pd )

πD2

4 . (92)

The Darcy friction factor, f, can be expressed as:

f = 8 τ cg

ρ 2 u . (93)

From the continuity equation,

u = 4m

ρ π 2 D . (94)

Substituting equations (93) and (94) into equation (92),

Pu − Pd = 8 fLm 2

cg ρ 2 π 5 D . (95)

Therefore,

Kf = 8 fL

cg ρ 2 π 5 D . (96)
150

 APPENDIX B—SUCCESSIVE SUBSTITUTION METHOD OF SOLVING COUPLED 

NONLINEAR SYSTEMS OF ALGEBRAIC EQUATIONS

The application of the successive substitution method involves the following steps:

(1) Develop the governing equations:

x1 = f

1 x1,x2,x3, ..., xn ( )

x2 = f

2 x1,x2,x3, ..., xn ( )

...

 xn = fn x1,x2,x3, ..., xn ( ) . (97)

If there are n number of unknown variables, there are n number of equations.

(2) Guess a solution for the equations:

Guess x1

*,x2*,x3*,...,xn* as an initial solution for the governing equations.

(3) Compute new values of x1,x2,x3,...,xn by substituting x1

*,x2*,x3*,...,xn* in the righthand side of equation (97).

(4) Under-relax the computed new value:

x = (1−α)x +αx, where a is the under-relaxation parameter.

(5) Replace x1

*,x2*,x3*,...,xn* with the computed value of x1,x2,x3,...,xn from step (4).

(6) Repeat steps (3)–(5) until convergence.
151

APPENDIX C—NEWTON-RAPHSON METHOD OF SOLVING COUPLED NONLINEAR 

SYSTEMS OF ALGEBRAIC EQUATIONS

The application of the Newton-Raphson method involves the following seven steps: 

(1) Develop the governing equations. The equations are expressed in the following form: 

f

1 x1,x2,x3, ..., xn ( ) = 0

f

2 x1,x2,x3, ..., xn ( ) = 0

...

 fn x1,x2,x3, ..., xn ( ) = 0 . (98)

If there are n number of unknown variables, there are n number of equations.

(2) Guess a solution for the equations. Guess x1

*,x2*,x3*, ..., xn* as an initial solution 

for the governing equations.

(3) Calculate the residuals of each equation. When the guessed solutions are substituted 

into equation (98), the right-hand side of the equation is not zero. The nonzero value is the 

residual. 

f

1 x1

*,x2*,x3*, ..., xn* ( ) = R1

f

2 x1

*,x2*,x3*, ..., xn* ( ) = R2

...

 fn x1

*,x2*,x3*, ..., xn* ( ) = Rn   . (99)

The intent of the solution scheme is to correct x1

*,x2*,x3*, ..., xn* with a set of corrections 

x′

1,x′2,x′3, ..., x′n such that R1,R2,R3, ...,Rn are zero.

(4) Develop a set of correction equations for all variables. First, construct the matrix 

of influence coefficients:
152

∂ f

1

∂x1

∂ f

1

∂x2

∂ f

1

∂x3

... ∂ f

1

∂xn

∂ f

2

∂x1

∂ f

2

∂x2

∂ f

2

∂x3

... ∂ f2

∂xn

... ... ... ... ...

∂ fn

∂x1

∂ fn

∂x2

∂ fn

∂x3

... ∂ f

n

∂xn

. (100)

Then, construct the set of simultaneous equations for corrections:

∂ f

1

∂x1

∂ f

1

∂x2

∂ f

1

∂x3

... ∂ f

1

∂xn

∂ f

2

∂x1

∂ f

2

∂x2

∂ f

2

∂x3

... ∂ f

2

∂xn

... ... ... ... ...

∂ fn

∂x1

∂ fn

∂x2

∂ fn

∂x3

... ∂ fn

∂xn

⎡

⎣

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎢

⎤

⎦

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

⎥

x′

1

x′

2

x′

n

⎡

⎣

⎢

⎢

⎢

⎢

⎢

⎤

⎦

⎥

⎥

⎥

⎥

⎥

=

R1

R2

Rn

⎡

⎣

⎢

⎢

⎢

⎢

⎢

⎤

⎦

⎥

⎥

⎥

⎥

⎥

 . (101)

(5) Solve for x′

1,x′2,x′3, ..., x′n by solving the simultaneous equations.

(6) Apply correction to each variable.

(7) Iterate until the corrections become very small.
153

REFERENCES

 1. Majumdar; A.K.; Bailey, J.W.; Schallhorn, P.A.; and Steadman, T.E.: Generalized Fluid 

System Simulation Program, U.S. Patent No. 6,748,349, June 8, 2004.

 2. Streeter, V.L.: Fluid Mechanics, 3rd ed., McGraw-Hill, 1962.

 3. Stoecker, W.F.: Design of Thermal Systems, 3rd ed., McGraw-Hill, 1989. 

 4. Owen, J.W.: “Thermal Analysis Workbook,” NASA TM–103568, Marshall Space Flight 

Center, AL, January 1992. 

 5. SINDA/G Thermal Analyzer, Network Analysis Inc., <http://www.sinda.com>, accessed 

January 2001.

 6. SINDA/FLUINT, Heat Transfer and Fluid Flow Design and Analysis Software, C&R Technologies, <http://www.crtech.com/sinda.html>, accessed January 2001.

 7. Seymore, D.C.: ROCETS Users Manual, Internal Report, NASA Marshall Space Flight 

Center, AL, October 1999. 

 8. Broyden, C.G.: “A Class of Methods for Solving Nonlinear Simultaneous Equations,” 

Mathematics of Computation, Vol. 19, No. 92, pp. 577–593, October 1965.

 9. Patankar, S.V.: Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corp., 

Washington, DC, 1980.

10. 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. 

11. 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, Glenn Research Center, 

OH, November 1973. 

12. Cryodata, Inc., User’s Guide to GASPAK, Version 3.20, November 1994. 

13. Winplot Version 4.54 User’s Manual, NASA Marshall Space Flight Center, AL, July 16, 2008. 

14. Holman, J.P.: Heat Transfer, 8th ed., McGraw Hill, 1997.
154

15. Miropolskii, Z.L.: “Heat Transfer in Film Boiling of a Steam-Water Mixture in Steam Generating Tubes,” Teploenergetika, Vol. 10, No. 5, pp. 49–52, 1963 (in Russian; translation Atomic 

Energy Commission, AEC-Tr-6252, 1964).

16. Colebrook, C.F.: ”Turbulent Flow in Pipes, with Particular Reference to the Transition 

Between the Smooth and Rough Pipe Laws,” J. Inst. Civil Eng., Vol. 11, pp. 133–156, 1938–

1939. 

17. White, F.W.: Viscous Fluid Flow, McGraw Hill, 1991.

18. Hooper, W.B.: “Calculate Head Loss Caused by Change in Pipe Size,” Chem. Engr., Vol. 95, 

No. 16, pp. 89–92, November 7, 1988. 

19. Yamada, Y.: “Resistance to a Flow through an Annulus with an Inner Rotating Cylinder,” 

Bulletin of the Japan Society of Mechanical Engineers, Vol. 5, pp. 302–310, 1962. 

20. Ito, H.; and Nanbu, K.: “Flow in Rotating Straight Pipes of Circular Cross Section,” Transactions of ASME, Vol. 93, Series D, No. 3, pp. 383–394, 1971.

21. Howell, G.W.; and Weather, T.M.: “Aerospace Fluid Component Designers Handbook,” 

Vol. 1, Rev. D, TRW Systems Group, RPL-TDR-64-25, February 1970.

22. Hooper, W.B.: “The Two-K Method Predicts Head Losses in Pipe Fittings,” Chem. Engr., 

Vol. 24, pp. 97–100, August 24, 1981.

23. Lebar, J.F.; and Cady, E.C.: “Viscojet Testing at NASA Lewis Research Center CCL-7 Test 

Facility,” Report No. A3-Y953-JFL-9G136, July 1990. 

24. Kays, W.M.; and London, A.L.: Compact Heat Exchangers, 3rd ed., McGraw-Hill, New York, 

1984.

25. Patankar, S.V.; and Spalding, D.B.: “A Calculation Procedure for Heat, Mass and Momentum 

Transfer in Three Dimensional Parabolic Flows,” Int. J. Heat Mass Transfer, Vol. 15, 

pp. 1787–1806, 1972. 

26. Hodge, B.K.: Analysis and Design of Energy Systems, 2nd ed., Prentice Hall, 1990. 

27. Schallhorn, P.A.; and Majumdar, A.K.: “Numerical Prediction of Pressure Distribution Along 

the Front and Back Face of a Rotating Disc With and Without Blades,” AIAA 97-3098, Presented at the 33rd Joint Propulsion Conference, Seattle, WA, July 6–9, 1997.

28. Schallhorn, P.A.; Elrod, D.A.; Goggin, D.G.; and Majumdar, A.K.: “A Novel Approach 

for Modeling Long Bearing Squeeze Film Damper Performance,” Paper No. AIAA 98-3683, 

Presented at the 34th Joint Propulsion Conference, Cleveland, OH, July 13–15, 1998.
155

29. Moody, F.J.: Introduction to Unsteady Thermofluid Mechanics, John Wiley, 1990.

30. Majumdar, A.K.; and Flachbart, R.H.: “Numerical Modeling of Fluid Transients by a Finite 

Volume Procedure for Rocket Propulsion Systems,” Paper No. FEDSM2003-45275, Proceedings of ASME FEDSM’03, 4th ASME/JSME Joint Fluids Engineering Conference, Honolulu, 

HI, July 6–10, 2003.

37. Majumdar, A.K.; and Van Hooser, K.: “A General Fluid System Simulation Program to 

Model Secondary Flows in Turbomachinery,” Paper No. AIAA 95-2969, 31st AIAA/ASME/

SAE/ASEE Joint Propulsion Conference and Exhibit, San Diego, CA, July 10–12, 1995.

38. Schallhorn, P.; Majumdar, A.; Van Hooser, K.; and Marsh, M.: “Flow Simulation in Secondary Flow Passages of a Rocket Engine Turbopump,” Paper No. AIAA 98-3684, 34th AIAA/

ASME/SAE/ASEE, Joint Propulsion Conference and Exhibit, Cleveland, OH, July 13–15, 

1998. 

39. Van Hooser, K.; Majumdar, A.; and Bailey, J.: “Numerical Prediction of Transient Axial 

Thrust and Internal Flows in a Rocket Engine Turbopump,” Paper No. AIAA 99-2189, 

35th AIAA/ASME/SAE/ASEE, Joint Propulsion Conference and Exhibit, Los Angeles, CA, 

June 21, 1999.

40. Reid, R.C.; Prausnitz, J.M.; and Poling, B.E.: The Properties of Gases and Liquids, 4th ed., 

McGraw Hill, 1987.

41. Epstein, M.; and Anderson, R.E.: “An Equation for the Prediction of Cryogenic Pressurant 

Requirements for Axisymmetric Propellant Tanks,” Adv. Cryog. Eng., Vol. 13, pp. 207–214, 

1968.

42. Van Dresar, N.T.: “Prediction of Pressurant Mass Requirements for Axisymmetric Liquid 

Hydrogen Tanks,” J. Propul. Power, Vol. 13, No. 6, pp. 796–799, 1997.
156

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining 

the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing 

this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operation and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid 

OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To)

4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S) 5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION

 REPORT NUMBER

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITOR’S ACRONYM(S)

11. SPONSORING/MONITORING REPORT NUMBER

12. DISTRIBUTION/AVAILABILITY STATEMENT

13. SUPPLEMENTARY NOTES

14. ABSTRACT

15. SUBJECT TERMS

16. SECURITY CLASSIFICATION OF:

a. REPORT b. ABSTRACT c. THIS PAGE

17. LIMITATION OF ABSTRACT 18. NUMBER OF

 PAGES

19a. NAME OF RESPONSIBLE PERSON

19b. TELEPHONE NUMBER (Include area code)

Standard Form 298 (Rev. 8-98)

Prescribed by ANSI Std. Z39-18

Generalized Fluid System Simulation Program,

Version 5.0—Educational

A.K. Majumdar

George C. Marshall Space Flight Center

Huntsville, AL 35812

National Aeronautics and Space Administration

Washington, DC 20546–0001

Unclassified-Unlimited

Subject Category 34

Availability: NASA CASI (443–757–5802)

Prepared by the Propulsion Systems Department, Engineering Directorate

Supplemental information is available on CD through CASI

M–1319

Technical Memorandum

NASA/TM—2011–216470

network flow, finite volume method, thermofluid system, two-phase, conjugate heat transfer

01–08–2011

UU 180

NASA

U U U

The Generalized Fluid System Simulation Program (GFSSP) is a finite-volume based general-purpose computer program for analyzing steady state and time-dependent flow rates, pressures, temperatures, and concentrations in a complex flow network. The 

program is capable of modeling real fluids with phase changes, compressibility, mixture thermodynamics, conjugate heat transfer between solid and fluid, fluid transients, pumps, compressors and external body forces such as gravity and centrifugal. The thermofluid 

system to be analyzed is discretized into nodes, branches, and conductors. The scalar properties such as pressure, temperature, and 

concentrations are calculated at nodes. Mass flow rates and heat transfer rates are computed in branches and conductors. The graphical user interface allows users to build their models using the ‘point, drag and click’ method; the users can also run their models and 

post-process the results in the same environment. The integrated fluid library supplies thermodynamic and thermo-physical properties of 36 fluids and 21 different resistance/source options are provided for modeling momentum sources or sinks in the branches. 

This Technical Memorandum illustrates the application and verification of the code through 12 demonstrated example problems. 

STI Help Desk at email: help@sti.nasa.gov

STI Help Desk at: 443–757–5802
The NASA STI Program...in Profile

Since its founding, NASA has been dedicated to the 

advancement of aeronautics and space science. The 

NASA Scientific and Technical Information (STI) 

Program Office plays a key part in helping NASA 

maintain this important role.

The NASA STI Program Office is operated by 

Langley Research Center, the lead center for 

NASA’s scientific and technical information. The 

NASA STI Program Office provides access to 

the NASA STI Database, the largest collection of 

aeronautical and space science STI in the world. 

The Program Office is also NASA’s institutional 

mechanism for disseminating the results of its 

research and development activities. These results 

are published by NASA in the NASA STI Report 

Series, which includes the following report types:

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant 

phase of research that present the results of 

NASA programs and include extensive data 

or theoretical analysis. Includes compilations 

of significant scientific and technical data 

and information deemed to be of continuing 

reference value. NASA’s counterpart of peerreviewed formal professional papers but has less 

stringent limitations on manuscript length and 

extent of graphic presentations.

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis.

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees.

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical conferences, 

symposia, seminars, or other meetings sponsored 

or cosponsored by NASA.

• SPECIAL PUBLICATION. Scientific, technical, 

or historical information from NASA programs, 

projects, and mission, often concerned with 

subjects having substantial public interest.

• TECHNICAL TRANSLATION. 

English-language translations of foreign 

scientific and technical material pertinent to 

NASA’s mission.

Specialized services that complement the STI 

Program Office’s diverse offerings include creating 

custom thesauri, building customized databases, 

organizing and publishing research results...even 

providing videos.

For more information about the NASA STI Program 

Office, see the following:

• Access the NASA STI program home page at 

<http://www.sti.nasa.gov>

• E-mail your question via the Internet to 

<help@sti.nasa.gov>

• Fax your question to the NASA STI Help Desk 

at 443 –757–5803

• Phone the NASA STI Help Desk at 

443 –757–5802

• Write to:

NASA STI Help Desk

NASA Center for AeroSpace Information

7115 Standard Drive

Hanover, MD 21076–1320
National Aeronautics and

Space Administration

IS20

George C. Marshall Space Flight Center

Huntsville, Alabama 35812