ITF 2000-1

DEVELOPMENT OF COMPUTER PROGRAM 

FOR FIRE SUPPRESSANT FLUID FLOW 

FINAL TECHNICAL REPORT 

VOLUME I 

MAIN REPORT

By 

Kemal Tuzla, Thomas Palmer and John C. Chen 

Lehigh University 

And 

Ramu K. Sundaram and Woon-Shing Yeung 

Duke Engineering and Services, Inc. 

September, 2000 

 Contractor: 

Institute of Thermo-Fluid Engineering & Science 

Lehigh University 

Bethlehem, PA 18015 USA 

Subcontractor: 

Duke Engineering and Services, Inc. 

Marlborough, MA 

1 
DISCLAIMER 

This report was prepared as a partial account of work performed for the National 

Institute of Standards and Technology by Lehigh University and Duke Engineering & 

Services, Inc. Neither Lehigh University nor Duke Engineering & Services, Inc., nor any of

their employees, makes any warranty, expressed or implied, or assumes any legal liability or 

responsibility for any third party’s use, or the results of such use, of any information, 

apparatus, product, or process disclosed in this report, or represents that its use by such third 

party would not infringe privately owned rights. 

2
TABLE OF CONTENTS 

DISCLAIMER 2

TABLE OF CONTENTS 3 

ACKNOWLEDGEMENTS 6

ABSTRACT 7

NOMENCLATURE 8

FIGURE CAPTIONS 11 

TABLE CAPTIONS 14 

1. DESCRIPTION OF PROJECT AND OBJECTIVES 15 

2. PREVIOUS WORKS AND SELECTION OF THE BASE CODE 17 

2.1. Experimental Data with Fire Suppressants 17 

2.2. Physical Property Data for Fire Suppressant Fluids 18 

2.3. Desired Feature of the New Code 18 

2.4. Two-Phase Codes Considered for this Work 19 

2.5. Selection of the Base Code 20 

3. EXPERIMENTAL WORK 21 

3.1. Experimental Test Facility 21 

3.1.1. Pressure Sensors 24 

3.1.2. Flow Rate Sensor 25 

3.1.3. Temperature Measurement Sensors 26 

3.1.4. Void Fraction Sensor 28 

3.2. Test Matrix Performed In Study 29 

3.3. A Typical Experimental Run 32 

3.4. Observations and Conclusions 33 

4. CODE DEVELOPMENT WORK 63 

4.1. INTRODUCTION 63 

4.1.1. SCOPE OF PROGRAM 63 

4.1.2. BACKGROUND 64 

4.2. CODE THEORY 66 

4.2.1. GOVERNING EQUATIONS 67 

4.2.2. NUMERICAL SOLUTION 70 

4.2.3. STATE RELATIONSHIPS AND FLUID PROPERTIES 71 

3
4.2.3.1. Basic Thermodynamic Property Tables 71 

4.2.3.2. Generation of Thermodynamic Property Tables 73 

4.2.3.3. Implementation Of Thermodynamic Property Tables 73 

4.2.3.4. Sample Binary Data File Generation for Halon 1301 74 

4.2.3.5. Critical Region Properties 77 

4.2.3.6. Transport Properties 78 

4.2.4. CONSTITUTIVE MODELS 80 

4.2.4.1. Wall Friction 81 

4.2.4.2. Interphase Friction 82 

4.2.4.3. Interphase Heat and Mass Transfer 83 

4.2.4.4. Wall Heat Transfer 83 

4.2.4.5. Coefficient of Virtual Mass 85 

4.2.4.6. Release of Dissolved Gas 85 

4.2.5. SPECIAL PROCESS MODELS 89 

4.2.5.1. Choked Flow 89 

4.2.5.2. Horizontal Stratification and Entrainment 89 

4.2.5.3. Abrupt Area Change 90 

4.2.5.4. User-Specified Form Losses 90 

4.2.5.5. Cross Flow Junction 91 

4.2.6. COMPONENT MODELS AND SYSTEMS 91 

4.2.6.1. Branch 91 

4.2.6.2. Valves 92 

4.2.6.3. Trip Systems 93 

4.2.6.4. Control Systems 93 

4.3. CODE ARCHITECTURE 102 

4.3.1. COMPUTER ADAPTIBILITY 102 

4.3.2. TOP LEVEL ORGANIZATION 102 

4.3.3. INPUT PROCESSING 103 

4.3.4. TRANSIENT PROCESSING 104 

4.4. CODE INSTALLATION AND EXECUTION 109 

4.4.1. SYSTEM REQUIREMENTS 110 

4.4.2. INSTALLATION 110 

4
4.4.3. FSP VERSION 1.0 FILES 112 

4.4.4. CODE EXECUTION 113 

4.4.4.1. Editing a File 113 

4.4.4.2. FSP Execution 114 

4.4.4.3. Rstrip Execution 115 

4.4.4.4. Ace Execution 116 

4.5. INPUT PREPARATION AND SAMPLE PROBLEM 119 

4.5.1. System Configuration 119 

4.5.2. Development of FSP Model 119 

4.5.3. Input File Listing 125 

4.5.4. Results & Data Comparison 130 

4.6. ASSESSMENT AGAINST EXPERIMENTS 135 

4.6.1. FSP ASSESSMENT AGAINST NIST EXPERIMENTS 135 

4.6.1.1. NIST Experimental Apparatus 135 

4.6.1.2. FSP Simulations of NIST Experiments 136 

5. Comparison of Code Predictions to Lehigh Experimental Data 145 

5.1. Comparison with Experimental Data Collected at Lehigh Test Facility 145 

6. SUMMARY 168 

7. REFERENCES 172

5
ACKNOWLEDGEMENTS 

This research was supported in part by the U.S. Department of Defense, through the

Strategic Environmental Research and Development Program (SERDP). Dr. Richard Gann, 

Director of Next Generation Fire Suppression Technology Program (NGP) and Dr. William

Grosshandler, Program Director, provided valuable guidance throughout the duration of the 

project. Dr. Jiann Yang of NIST provided technical expertise in the experimental work and 

also needed information on fluid properties. Their support is gratefully acknowledged. 

It is noted that the present computer program was developed from a base computer 

program, RELAPS/MOD3, developed originally at the Idaho National Engineering 

Laboratory for nuclear reactor safety analysis under the sponsorship of the U.S. Nuclear 

Regulatory Commission. 

6
ABSTRACT 

A computer program has been developed for predicting the hydraulic performance of 

delivery systems of fire-suppressant fluids. The program was developed under the Next 

Generation Fire Suppression Technology Program being managed by the National Institute 

of Standards and Technology. 

The program is based on a one-dimensional, two-fluid model of two-phase flow. In this 

model, separate conservation equations are written for the liquid and gas phases for mass, 

momentum and energy. Constitutive relationships are specified for interphase transport of

mass, momentum and energy. Heat transfer between the fluid and passive structures such as 

pipe walls can be modeled. The program also contains built-in models for wall friction and 

two-phase critical flow. The transport of non-condensable gas in the system, namely the 

Nitrogen fill gas, as well as the Nitrogen released from solution during agent discharge, is 

modeled via separate mass conservation equations, with constitutive relations to specify the 

rate of gas release. The conservation equations are solved using a semi-implicit numerical 

method, with user-supplied boundary and initial conditions. 

The program was deliberately made flexible in terms of types of fluids and piping 

layout. The current version of the program allows the user to select any one of five fluids – 

water, Halon-1301, CO2, HFC-227ea or HFC-125. Modules are available in the program with 

which the user can model a delivery system, including one or more supply tanks and a 

combination of piping networks. The user can also model valves in the system, with specified

valve opening times if needed. 

The suppressant discharge is a highly transient process, generally lasting less than a 

second to few seconds. The program was benchmarked against transient experimental data 

available in the literature, as well as experiments conducted as a part of this project, on the 

discharge of HFC-227ea and HFC-125 in a specially prepared discharge loop. The present 

experiments lasted between 1.5 to 6 seconds. In addition to transient system pressure at 

various points, these experiments also measured critical parameters, such as the transient 

mass discharge history, fluid temperature and the void fraction near the exit. These are first 

dynamic measurements of mass flow rate, fluid temperatures and void fraction during 

suppression discharge. This allowed a more comprehensive assessment of the computer

program than possible with previously available experimental data. The results of the 

assessment showed that the program is capable of predicting the performance of various 

delivery systems with several fluids. 

This report is prepared in two volumes. First volume provides a detailed description of 

the experimental work and discussion of the results. It also includes a brief description of the

theory and numerical solution method, instructions for installation on a personal computer 

with the WINDOWS operating system, and instructions on preparation of the input needed to 

describe the system being analyzed. Volume II of the report consists of three appendices: a)

an input manual for the computer program, b) information on the experimental data, and c) 

comparison of the present data to predictions of the computer program. 

7
NOMENCLATURE

A area, cross-sectional area (m2

) 

At throat area (m2

) 

a interfacial area per unit volume (m-1) 

B drag coefficient, momentum conductance 

Bx body force in x coordinate direction (m/s2

) 

C coefficient of virtual mass 

Cp specific heat at constant pressure (J/kg.K) 

Cv specific heat at constant volume (J/kg.K), valve flow coefficient 

D diameter (m) 

DISS energy dissipation function (W/m3

) 

FIF, FIG interphase drag coefficients (liquid, vapor) (s -1) 

FI interphase drag coefficient (m3/kg.s) 

FWF, FWG wall drag coefficients (liquid, vapor) (s -1) 

fF friction factor of pipe 

G mass flux (kg/m2

 -s) 

g gravitational constant (m/s2

 ) 

H volumetric heat transfer coefficient (W/K.m3

) 

h specific enthalpy (J/kg), heat transfer coefficient (W/m2

 .K), energy transfer 

 coefficient for Γg 

K energy form loss coefficient 

k thermal conductivity (W/m.K) 

L length 

M molecular weight, mass (kg) •

m instantaneous rate of mass flow (kg/s) 

P pressure (Pa), channel perimeter (m) 

Pi initial nitrogen pressure in supply vessel (Pa) 

PN partial pressure of Nitrogen in gas phase 

Po atmospheric pressure (Pa) 

P1,P2,P3 pipe position #1, #2, #3 pressure (k/Pa) 

p wetted perimeter (m) 

Pr Prandtl number 

PT1,PT2. . . pressure taps for pressure drop measurements 

Q volumetric heat addition rate (W/m3

) 

q heat transfer rate (W), heat flux (W/m2

) 

R radius (m) 

Rn, Rs gas constants (noncondensable, steam) (N.m/kg.K) 

rc critical gas bubble radius (m) 

s specific entropy (J/kg.K) 

T temperature (K), trip 

Tc critical temperature (K) 

TR reduced temperature (K) 

T1,T2,T3 pipe position #1, #2, #3 temperature 

t time (s) 

U specific internal energy (J/kg) 

 volume (m3

), specific volume (m3/kg) 

8
v mixture velocity (m/s), phasic velocity (m/s) 

vc choking velocity (m/s) 

X flow quality, static quality, mass fraction 

x spatial coordinate (m) 

xN dissolved gas concentration (mol/mol)

∆Z height of volume (m) 

z elevation change coordinate (m) 

Symbols 

α void fraction 

β coefficient of isobaric thermal expansion (K-1) 

Γ volumetric mass exchange rate (kg/m3

 .s) 

∆P pressure drop (k/Pa) 

∆PNUC pressure difference at initiation of gas release (Pa) 

∆t increment in time variable (s) 

ε wall vapor generation/condensation flag 

κ coefficient of isothermal compressibility (Pa-1) 

μ viscosity (kg/m.s) 

ν kinematic viscosity (m2

/s) 

ρ density (kg/m3

), 

σ surface tension (J/m2) 

υ specific volume (m3

/kg) 

Subscripts

CHF value at critical heat flux condition 

cr, crit critical property or condition 

cv collection vessel 

f liquid phase 

fg phasic difference (i.e., vapor term-liquid term) 

g vapor phase 

I interface 

i interface 

j, j+1, j-1 spatial noding indices for junctions 

K spatial noding index for volumes 

k numbering index for volumes 

L spatial noding index for volume 

M spatial noding index for volume 

m mixture property 

min minimum value

ref reference value 

s steam component of vapor phase, superheated, solute 

sat saturated quality 

T turbulent 

v vapor quantity 

w wall 

wall wall 

9
wg, wf wall to vapor, wall to liquid 

μ viscosity 

 

Superscripts

s saturation property 

10
FIGURE CAPTIONS 

Figure 3.1: Schematic of Test Facility: Source Vessel, Discharge Piping, and Collection 

 Vessel. 

Figure 3.2: Schematic of Source Vessel Showing Mass Inventory Measurement. 

Figure 3.3: Pressure versus Temperature Saturation Curve for HFC-227ea & HFC-125. 

Figure 3.4: Marotta Valve. 

Figure 3.5: Top View of Test Facility, schematic of Discharge Piping System. 

Figure 3.6: Schematic of Pressure Fittings for Temperature and Pressure Ports. 

Figure 3.7: Schematic of Rack Used to Hold the Source Vessel. 

Figure 3.8: Picture of Test Facility: Source Vessel, Discharge Piping, Collection Vessel, 

 Support Structure and Instrumentation. 

Figure 3.9: Pressure Histogram for Experimental Run #5. 

Figure 3.10: Plot of Source Vessel Pressure Showing the Repeatability of the 

Experimental Runs. 

Figure 3.11: Schematic of Alternate Piping Configuration #1. 

Figure 3.12: Schematic of Alternate Piping Configuration #2. 

Figure 3.13: Pressure Drop Histogram Across a 90° Elbow from Run #9. 

Figure 3.14: Diagram Transducer Setup for Mass Inventory Measurement. 

Figure 3.15: Mass Inventory Histogram of Run #5 

Figure 3.16: Schematic of Thermocouple Construction. 

Figure 3.17: Temperature Histogram from Film Thermocouples of Run #5. 

Figure 3.18: Temperature Histogram from Shielded Thermocouples of Run #5. 

Figure 3.19: Schematic of Capacitance Probe for Bench Setup. 

Figure 3.20: Bench Setup to Test Capacitance Sensor. 

Figure 3.21: Void Fraction from Bench Tests of Capacitance Sensor. 

Figure 3.22: Schematic of Capacitance Probe for Test Facility. 

Figure 3.23: Cross-Section of Capacitance Setup Showing Electrical Field. 

Figure 3.24: Void Fraction Histogram of Run #5. 

Figure 3.25: Comparison of Experimental Run Length with Straight Pipe (#1) and a Pipe 

with Fittings (#8). 

Figure 3.26: Comparison of Pressure Drop Across the Top and Bottom Branches of a Tee 

Split From Experimental Run #13. 

Figure 3.27: Comparison of Collection Vessel Temperature with HFC-227ea (#5) and 

HFC-125 (#16) as a Test Fluid. 

Figure 3.28: Cumulative Mass for Experimental Run #18 Using Water as a Test Fluid. 

Figure 3.29: Temperature and Void Fraction Histograms from Experimental Run #19 

Using Water as a Test Fluid. 

Figure 4.2.1: Staggered mesh numerical method. 

Figure 4.2.2: Flow Chart of Generation Program. 

Figure 4.2.3: Curve Fit for Henry Constant. 

Figure 4.2.4: Effect of Release Coefficient on Calculated Pressure Response. 

Figure 4.2.5: Comparison of Calculated and Measured Bottle Pressure for Test 146. 

Figure 4.2.6: Void Fraction Responses in Lower Portion of Bottle. 

Figure 4.2.7a: 90-degree Tee model with cross flow junction. 

Figure 4.2.7b: Tee model with normal junctions. 

11
Figure 4.3.1: Top Level Structure. 

Figure 4.3.2: Modular Structure of Transient Calculations. 

Figure 4.3.3: Transient (Steady-State) Block Structure. 

Figure 4.5.1: Schematic of Lehigh Test Apparatus. 

Figure 4.5.2: FSP Nodalization Diagram for Lehigh Tests. 

Figure 4.5.3: Comparison of Calculated and Measured Bottle Pressure for Sample 

Problem. 

Figure 4.6.1: NIST Test Apparatus (adapted from Gann et al, 1995). 

Figure 4.6.2: Piping Configurations for NIST Experiments. 

Figure 4.6.3: FSP Input Model. 

Figure 4.6.4: Comparison of Calculated and Measured Bottle Pressure for B-59. 

Figure 4.6.5: Comparison of Calculated and Measured Bottle Pressure for B-68. 

Figure 4.6.6: Comparison of Calculated and Measured Bottle Pressure for B-61. 

Figure 4.6.7: Comparison of Calculated and Measured Bottle Pressure for B-60. 

Figure 5.1: Comparison of Experimental and Predictive Source Vessel Pressure Data in 

Sensitivity Study on Critical Radius Value for Conditions of Run #4. 

Figure 5.2: Comparison of Experimental and Predictive Source Vessel Pressure Data in 

Sensitivity Study on Gas Release Rate Value for Conditions of Run #4. 

Figure 5.3: Comparison of Experimental and Predictive Void Fraction Data in 

Sensitivity Study on Gas Release Rate Value for Conditions of Run #4. 

Figure 5.4: Comparison of Experimental and Predictive Pressure Data using Final 

Values for Code Operating Parameters for Conditions of Run #4. 

Figure 5.5: Comparison of Experimental and Predictive Void Fraction Data using Final 

Values for Code Operating Parameters for Conditions of Run #4. 

Figure 5.6: Comparison of Cumulative Mass for Experimental Run #5 and Predictive 

Computer Code Data. 

Figure 5.7: Comparison of Void Fraction for Experimental Run #5 and Predictive 

Computer Code Data. 

Figure 5.8: Comparison of Source Vessel Pressure for Experimental Run #5 and 

Predictive Computer Code Data. 

Figure 5.9: Comparison of Pipe Position #1 Pressure for Experimental Run #5 and 

Predictive Computer Code Data. 

Figure 5.10: Comparison of Pipe Position #2 Pressure for Experimental Run #5 and 

Predictive Computer Code Data. 

Figure 5.11: Comparison of Pipe Position #3 Pressure for Experimental Run #5 and 

Predictive Computer Code Data. 

Figure 5.12: Comparison of Collection Vessel Pressure for Experimental Run #5 and 

Predictive Computer Code Data. 

Figure 5.13: Comparison of Pressure Drop Across a Capped Tee for Experimental Run 

#9 and Predictive Computer Code Data. 

Figure 5.14: Comparison of Pressure Drop Across a 90° Elbow for Experimental Run #9 

and Predictive Computer Code Data. 

Figure 5.15: Comparison of Pressure Drop Across a Union for Experimental Run #9 and 

Predictive Computer Code Data. 

Figure 5.16: Comparison of Source Vessel Fluid Temperature for Experimental Run #16 

and Predictive Computer Code Data. 

Figure 5.17: Comparison of Pipe Position #1 Fluid Temperature for Experimental Run 

#16 and Predictive Code Data. 

12
Figure 5.18: Comparison of Collection Vessel Gas Temperature for Experimental Run 

#16 and Predictive Computer Code Data. 

13
TABLE CAPTIONS 

Table 3.2: Matrix of Experimental Runs Conducted at the Lehigh Test Facility. 

Table 4.1: PROFISSY Results. 

Table 4.4.1: Contents of the WinZip Archive fspv1.zip 

14
1. DESCRIPTION OF PROJECT AND OBJECTIVES

The objective is to develop a computer program (or code) capable of predicting the 

delivery rates of a fire suppressant fluid to a designated location, via a piping network 

system. 

The current delivery systems for Halon generally consist of a suppressant vessel 

connected to the delivery location by a piping network. The suppressant vessel is pressurized 

by an inert gas, usually nitrogen. The suppressant fluid is maintained in a liquid state by the 

pressurization of the system, usually up to several MPa. The system is activated by a quick 

opening valve in response to either an automatic or manual trigger. This causes the gas 

driver to expel the suppressant fluid through the piping system. The delivery system has to

be able to begin delivering the fluid to the discharge location very quickly, generally within a 

few seconds. Complete discharge of the system is required to occur within 0.1 to 10 seconds. 

Since Halon and its potential replacement fluids are in a vapor state at standard 

pressure and temperature, at the discharge location the fluid would be a two-phase mixture of 

superheated liquid and vapor. Thus significant thermodynamic and thermal non-equilibrium

can be expected between the phases. Also, due to the large pressure difference between 

upstream source vessel and downstream exit, continuous flashing is anticipated as the fluid 

travels through the piping, and two-phase critical (choked) flow can occur at various 

locations. 

Hence, to accurately predict performance of the delivery system, the computer code 

should be a transient two-phase code, which allows for phase non-equilibrium. The short 

delivery times require fairly high flow rates, which promote homogeneous two-phase flow, 

i.e. little slip. However, in imbalanced piping networks, with side tees and other fittings,

some separated flow could occur (stratified flows, slug/plug flow, etc.). Hence the code

should be able to predict slip between phases and the corresponding effect on pressure drop. 

This can be especially important when separated flow encounters directional change, such as 

at a side tee. The ability to predict the transport of non-condensable gas is also important. 

The fluid, pressurized with a driver gas, is initially saturated with the gas. During the 

delivery, as the system depressurizes, the driver gas comes out of solution and expands. This 

gas evolution phenomenon needs to be accounted for. One additional requirement is that the 

code should be useful for estimation of the transient hydrodynamic loads in the piping 

network. Thus the momentum equations need to be sufficiently detailed to estimate the

unbalanced force in piping sections between locations of elbows and other fittings. 

There has been some work conducted on modeling Halon 1301 flows. In a 1988 

summary of the state-of-the-art, DiNenno and Budnick (1988) indicate that most of the flow 

calculations are performed by vendors using proprietary methods. These methods are 

verified against standard NFPA procedures for calculations and by discharge to testing in 

sealed tests. However, there are significant uncertainties in the calculations and also in

extrapolation of test results to actual configurations. 

15
Therefore, the objective of the present work was to develop a computer code capable 

of predicting single and two-phase hydrodynamic behavior of fire suppressant fluids during 

transport through delivery piping systems. 

A key technical approach in the present program was utilization of advancements 

made in other applications that deal with multi-phase flows. In particular, the highly 

sophisticated computer codes that have been developed for thermal-hydraulic analysis of 

nuclear power systems have all the characteristics required for analysis of fire suppressant 

systems. These include models that account for relative slip between liquid and vapor 

phases, thermodynamic non-equilibrium between the phases, changes in two-phase flow 

regimes, critical choke flows, and transport of non-condensable gases. Such codes are also 

structured for numerical analysis of fast transients, well capable of the transients anticipated

for suppressant systems. 

Furthermore, since the fluids being considered for use in fire suppressant systems is 

different from the water/steam system of nuclear power systems, physical property packages 

need to be modified and expanded. Here again, we have leveraged existing technology by 

borrowing from existing property packages. Specifically the REFPROP package, developed

by the National Institute of Science and Technology (NIST), covers many refrigerant fluids

including those that are potential candidates for the fire suppressant systems. 

The present project contained an experimental task, to obtain data needed to assess 

the code. The experimental program envisioned a discharge loop using new fire suppressants. 

In this experimental program, we successfully measured major flow parameters, which have 

not been measured heretofore. These include measurements of instantaneous discharge flow 

rate, fluid temperature, and void fraction at various locations along the discharge pipe. 

This program has been carried out as a team effort; the Institute of Thermo-Fluid

Engineering and Science at Lehigh University was the main contractor. Duke Engineering 

and Services is subcontractor to Lehigh University and responsible for code development 

work. 

16
2. PREVIOUS WORKS AND SELECTION OF THE BASE CODE 

2.1 Experimental Data with Fire Suppressants 

To fully test the computer code developed in this research, comparisons of code 

predictions against experimental data will be made. A literature survey has been conducted

to find possible sources for experimental data for transient two-phase flows. One widely 

referenced study is the "The Marviken Tests” conducted using water/steam for the nuclear 

reactor industry, Marviken Power Station (1982) and EPRI Report NP-2370 (1982). These 

tests used a liquid/vapor filled vessel under high pressure and measured the vertical discharge 

characteristics as this system was released to atmospheric conditions. A similar test was 

carried out by Edwards and O'Brien (1979), but focused on the effects of a horizontal 

discharge. Even though the Marviken and Edwards tests were conducted with water/steam

systems, they still provide some means for assessing portions of the code, which use 

mechanistic models for predicting non-homogeneous, non-equilibrium two-phase flows. 

These mechanistic models are fluid-independent and are invoked for all fluids.

Three studies, which more closely reflect the intentions of the current study, have 

used fire suppressants as their test fluids, rather than water and steam. The work undertaken 

by Elliot et al. in 1984, include both theoretical and experimental investigations. In this 

work, Halon 1301 (widely used in both military and commercial applications) was used as

the suppressant fluid. The tests discharged Halon under nitrogen gas pre-pressure, resulting 

in flow of a liquid/vapor/non-condensable-gas mixture through a configuration of piping and 

an end nozzle. This study provides data on pressure decay as a function of time, allowing for 

the estimation of discharge flow rates, which are useful for accessing code predictions. This 

study also provides observations of experimental phenomena as the discharge proceeded. An 

example is the phenomenon of evolution of the dissolved nitrogen from the liquid Halon, 

which causes an increase in the pressure of the discharge vessel and an increase in the driving 

force propelling the discharge of the fluid. 

The most comprehensive tests conducted to this date were launched in the mid 1990's

by the National Institute of Standards and Technology, Pitts et al. (1994) and Yang et al. 

(1995). Their tests not only utilized Halon 1301 as a test fluid, but also utilized HFC-125, 

HFC-227ea, and CF3I. They conducted several experimental runs with each fluid, providing 

extensive data for various pipe configurations and operating conditions. Experiments with 

different initial fill volumes, pressures, and temperatures were performed to determine their 

effects on the two-phase system. Their findings provide important information about the 

characteristics of the two-phase transient flow as it discharges from the source vessel and 

through the piping. The pressure histograms from this earlier program will be used as a 

benchmark in gauging the effectiveness of the code being developed in the current program. 

These prior studies provide the beginning of an experimental database. However, 

there are no data on several key parameters important to code assessment. Primary among 

these is the instantaneous mass discharge rate and void fraction, at the exit of the source 

vessel and at the exit of the distribution pipe. Such parameters reflect the mechanisms of 

phase change (vaporization) and interfacial momentum transfer (slip velocity) that must be 

correctly modeled in the computer code. Without specific experimental measurements of 

such key parameters, it is impossible to fully assess the validity of any prediction. This is the 

17
reason why the experimental program at Lehigh has been structured to obtain a number of 

additional measurements, above and beyond the usual measurement of pressure decay. 

2.2 Physical Property Data for Fire Suppressant Fluids 

Two-phase flow calculations require thermodynamic and transport properties of the 

fluid. These properties are required for single-phase liquid, for single-phase vapor, and for 

saturated vapor/liquid mixtures. The following thermodynamic properties are needed as

functions of absolute pressure: 

• Saturation temperature 

• Specific heat for gas and liquid phases 

• Thermal expansion coefficient for gas and liquid phases 

• Isothermal compressibility for gas and liquid phases 

• Specific volume for gas and liquid phases 

• Specific entropy for gas and liquid phases 

In addition, the following transport properties are also required for two-phase flow 

calculations: 

• Dynamic viscosity for gas and liquid phases 

• Thermal conductivity for gas and liquid phases 

• Surface tension 

Most fire suppressants are also used as refrigerants. The refrigerant property package 

REFPROP developed at NIST by Gallagher et al. (1993) covers most of the fire suppressants. 

For example, it contains properties of Halon 1301, and present fire suppressant candidate FC227ea. After evaluation, we find that REFPROP is a suitable package for determining the 

properties required by the new code for two-phase flow calculations with suppressant fluids. 

2.3 Desired Feature of the New Code 

As discussed in the original proposal, Chen, Tuzla, Sundaram (1997), the base code

should be developed with some specific features in terms of its technical ability to model the 

basic phenomena of two-phase flows. In addition, user convenience and ease of

development are also desirable features. The present development seeks to build the 

following features into the new code: 

Required Features 

Desirable Features 

1. Non-proprietary, publicly available 

2. Portable to various computer and operating systems 

3. Ease of development (well-documented, flexible architecture) 

1. Non-homogeneous two-phase flow accounting for relative slip between phases 

2. Non-equilibrium two-phase flow, allowing for metastable conditions during flashing 

18
3. Choking under single and two-phase conditions 

4. Non-condensable (dissolved gas) transport 

5. Multi-component mixtures (liquids, gases, powders) 

6. Flexible system modeling capability (choice of piping components and layout) 

7. Robust set of equations, for simulation of fast transient two-phase flows (for example,

momentum equation to include frictional losses, energy equation to include

dissipation) 

8. Demonstrated ability to model representative systems (assessment against 

experiments) 

4. Ease of application (user-friendly inputs) 

There are other intangible features that have also been considered, such as the existence of an 

active user group, which would facilitate future modifications and developments.

2.4 Two-Phase Codes Considered for this Work 

The relevant codes are those used in the nuclear industry for loss-of-coolant accident 

analyses (LOCA). These codes address the major phenomena of interest (non-homogeneous, 

non-equilibrium two-phase flows, choking conditions, fast transients). The survey was

limited to codes developed and available in the United States and consisted of RELAP5, 

RETRAN, TRAC-P, TRAC-B and GOTHIC. Codes developed in other countries that could 

be used (ATHLET, CATHARE) were not considered because easy access to these codes may

not be available. 

All these codes have many of the required and desirable features for the present 

application. They are all based on a two-fluid model of two-phase flow, with constitutive 

equations to specify the interfacial transport terms. Although some of them have multidimensional modeling capability, they are predominantly used in a one-dimensional mode. 

A brief history of these codes is presented below, with some relevant information: 

1. RELAP5 

The latest version of this code is RELAP5/MOD3, version 3.2, available since 1995, 

RELAP5/MOD3 Code Manual (1995). This code was developed at the Idaho National 

Engineering and Environmental Laboratories (INEEL) under sponsorship of the U.S. Nuclear 

Regulatory Commission (USNRC). It is used for a variety of reactor accident analyses, 

including LOCA analyses, for a variety of types of reactors. It is also used in a variety of 

general thermal-hydraulic analyses in nuclear and non-nuclear systems. It is non-proprietary, 

and is publicly available, for a transmittal fee, from the USNRC. An active user group, 

including U.S. and international organizations, meets once a year. Information on user 

problems, error correction, etc. is transmitted to user group members once every three 

months. 

2. RETRAN 

The latest version of this code is RETRAN-03, RETRAN-03 (1992). This code was

developed by Computer Simulations and Analysis, Inc., under sponsorship of the Electric 

Power Research Institute (EPRI). It is used primarily for analyses of reactor systems, with 

19
focus on the primary cooling system, and usually for non-LOCA scenarios. It is proprietary, 

and can be obtained from EPRI, with license fees and user group membership fees. The user 

group, including U.S. and international organizations, is active and meets regularly. 

3. TRAC-P 

The latest version of this code is TRAC-PF1/MOD2, TRAC-PF1/MOD2 Code (1992). It 

was developed by the Los Alamos National Laboratory (LANL) under USNRC sponsorship. 

It is used primarily for LOCA analysis of pressurized water reactors. It is non-proprietary 

and can be obtained, for a transmittal fee, from the USNRC. There are relatively small 

number of users and no regular user group meetings. 

4. TRAC-B 

The latest version of this code is TRAC-BF1/MOD1, TRAC-BF1/MOD1 Code (1992) and 

TRAC-BF1/MOD1 Models (1992). It was developed by INEL under USNRC sponsorship. 

It is used for accident analyses of boiling water reactors. It is non-proprietary and can be

obtained from the USNRC for a transmittal fee. The user group is active and meets

regularly. 

5. GOTHIC

The latest version of this code is GOTHIC version 5.0, GOTHIC Containment Analysis 

Package Version 5.0 (1995). It was developed by Numerical Applications, Inc., under EPRI 

sponsorship. It is used primarily for analysis of reactor containments. It was developed from 

the COBRA series of codes originally developed at the Pacific Northwest Laboratory under 

USNRC sponsorship. The original COBRA series of codes are non-proprietary, but were not 

maintained by the USNRC after the mid-1980's, and are now difficult to obtain. The 

GOTHIC code is proprietary and can be obtained from EPRI for license fees and user

membership fees. The user group is active and meets regularly. 

2.5 Selection of the Base Code 

The documentation for each of the candidate codes was reviewed. It was concluded 

that all the codes listed above have many of the required and desirable features, but none 

have all the features. The conclusion was that the RELAP5 code (MOD3, version 3.2) has 

more of the required and desirable features than the other codes. In addition, the code has 

been adapted previously for operation with a fluid other than water, and the code architecture

to accomplish this is available in the latest version of the code. This is especially helpful for

insertion of new property packages for suppressant fluids. Besides, the code development 

team at Duke Engineering and Services has accumulated significant operating experience 

with RELAP5, thus providing confidence that the program objectives could be met 

successfully with this code. Based on these characteristics, RELAP5/MOD3 version 3.2 has 

been selected as the base code platform for the present development. 

20
3. EXPERIMENTAL WORK 

The objective for the experimental work was to collect experimental data using new 

fire suppressant agents as test fluids. This data will then be used in the assessment of the 

product code for simulating possible suppressant delivery systems. Descriptions of pertinent 

dimensions and design features of the test facility are provided below. 

3.1 Experimental Test Facility 

The test facility consists of a source vessel for the suppressant fluid, a quick opening

discharge valve, a holding rack, a piping network, and a collection tank. This arrangement

allows the system to be operated as a closed loop, with recycle of the test fluid. A simple

schematic has been provided in figure 3.1, showing the layout of the facility. The 

dimensions of this test system adhere as closely as possible to the standard sizes used in

present-day fire suppression systems. 

The source vessel, one of the main components of the facility, is shown in figure 3.2. 

The source vessel has a total internal volume of 3.81 liters and was initially charged with

liquid to about 60% of that volume. The source vessel was constructed from a 2 meter long 

piece of 2-inch (5.08 cm) nominal diameter stainless steel pipe with an internal diameter of 

1.94 inches (4.93 cm) and a wall thickness of 0.22 inches (0.55 cm). Using a vessel of large 

length-diameter ratio improves the accuracy of measuring liquid inventory. This will be 

explained in more detail in Section 3.1.2, but simply speaking, we will use a pressure 

difference associated with the static head of the liquid in the vessel to measure the liquid

inventory in the vessel at any time during the experimental run. The long length, or height, 

of the source vessel provides a greater range of static head measurements, thus increasing the 

accuracy of determining the transient liquid inventory. This pipe will be positioned vertically 

and capped at both ends by stainless steel pipe flanges approximately 1 inch (2.54 cm) thick. 

Two slip-on flanges with o-ring grooves for the pressure seals were welded to either end of

the pipe so that the top and bottom plates could be bolted in place. The top plate has an 

identical outer diameter as the flanges, at 6.5 inches (16.51 cm), and was tapped to allow for 

a mixing line inlet and a fill/relief line. The bottom plate has a larger diameter of 11 inches 

(27.94 cm) for the purpose of anchoring the vessel to the holding rack minimizing any 

vibration or movement during operation. This plate has also been tapped to allow the 

discharge valve to be screwed into place and sealed with a crushable, tempered aluminum 

gasket. Other features incorporated into the source vessel include pressure and temperature

ports, a liquid level gauge, a mixing line, and the differential pressure transducer with fluid 

transmission line for measurement of liquid inventory. 

Since tests will require pressures up to 40-50 atmospheres in the source vessel,

pressurization by a nitrogen cover gas will be employed. As shown in figure 3.3, for FC227ea, the equilibrium vapor pressure of suppressant fluids will only reach about 5 

atmospheres at room temperature, insufficient for the desired pressure without the additional 

partial pressure of an inert cover gas such as nitrogen. A potential problem with this method 

is that the total pressure in the vessel will not reach equilibrium until the liquid has become

saturated with dissolved nitrogen. In order to expedite this process, a mixing system has

been incorporated into the design, using a pump to recycle liquid suppressant from the

21
bottom of the source vessel to a spray nozzle located in the vapor/gas space. By spraying

small droplets of liquid through the nitrogen gas, the surface area of liquid in contact with 

nitrogen will increase, therefore increasing the rate at which nitrogen can be dissolved into 

the liquid. 

Additional ports in the source vessel include an opening near the top, which will 

serve as the connection to the nitrogen supply and as one leg of the transmission line for the 

differential pressure transducer. A liquid level gauge made of armored glass has also been 

attached to the middle section of the vessel to visually determine the liquid fill level in the

vessel at the beginning of each experimental run. Also, there are two thermocouple ports for 

the vapor and liquid spaces and a pressure port in the vapor space. The last port is in the top 

plate, providing access for filling of the vessel and also serving as a connection for a safety

pressure-relief device. 

The next pertinent item in our design is the release mechanism of the source vessel, 

for initiating discharge of the test agent. A quick-opening plunger valve with an inlet 

diameter of 1.25 inches (3.18 cm) and an exit diameter of 1.75 inches (4.45 cm) was used for

this purpose. This valve, shown in figure 3.4, is a Model MV121KJ-2 valve made by 

Marotta Scientific Controls Inc. A feature of this valve is the use of the fluid pressure to 

drive the plunger upon command; thereby eliminating the need for springs or motors. To set 

or close the valve, the plunger is pushed up into its armed position and held in place by a 

small locking latch. Two O-rings are attached to the plunger to form a pressure seal along 

the valve casing. Pressure must then be provided on the inlet side of the valve for the valve 

to operate properly. To open the valve, an electrical charge provided by 20-volt DC power 

source releases the latch, allowing the fluid pressure to rapidly drive open the plunger. This 

release is much faster than that of traditional solenoid valves, and is similar to that of a burst 

diaphragm or explosive-charge release, which is used in present-day suppression systems.

The valve is connected to the 1/2-inch discharge piping in the loop by a tapered transition 

piece. This transition incorporates a flange connection to permit access to the O-ring seals on

the plunger for lubrication. 

The piping in the test facility consists of 1/2-inch nominal diameter stainless steel 

pipe with an internal diameter of 0.546 inches (1.39 cm) and a wall thickness of 0.15 inches 

(0.38 cm). There are five other pieces including the adapter piece described above, which 

make up the discharge piping system. Two of these pieces are approximately 12 inches (30.5 

cm) in length with flanges at either end to provide for easy disassembly. All the flanges are 

machined with an O-ring groove to provide an adequate seal against the high pressures 

experienced in the tests. These 12-inch sections will be used as access points at the beginning 

and end of the piping system, as indicated in figure 3.5. They will also be the first and last

locations of the pressure and temperature measurement ports. Ports for the pressure 

transducers and film thermocouples used for these measurements required special weld 

fittings in order to insure a complete pressure seal in the loop. These fitting were specially

designed to cause minimal disturbance to the flow, while allowing for the sensor to be as 

close to the flow as possible. See figure 3.6 for a schematic diagram of these fittings. 

The other two pieces in the piping system include a special section used for the void 

fraction measurement and a long middle section of pipe. The section used in the

measurement of void fractions will be described in section 3.1.4. The long middle section of 

22
pipe represents the distribution spool piece often found in fire-suppressant piping systems. 

Both ends are flanged and are capped by valves, allowing for the isolation of this section 

when dismantling. The middle piece will be tapped in the middle for temperature and 

pressure ports as described previously. This middle piece of long discharge piping allows for

easy reconfiguration to different piping schemes (i.e. inclusion of tees or elbows). This whole 

piping system is then connected by a flange to the collection vessel, providing a closed flow 

path from the source vessel, through the piping system, to the collection vessel. 

The collection vessel is designed to serve two purposes: to catch the test fluid for 

recycle and to enable monitoring of the transient discharge pressure during test runs. The 

collection vessel has an internal volume of approximately 30 liters, chosen to limit the final 

pressure after discharge to less than 10-12 atmospheres. This vessel is made of a 10 inches 

diameter stainless steel pipe with an internal diameter of 9.56 inches (24.3 cm) and a wall 

thickness of 0.59 inches (1.50 cm). The bottom of the collection vessel is a hemispheric cap 

with a 1/2-inch port, to be used for emptying the contents. At the top end is a pair of 1.75 

inches thick flanges sealed by an O-ring. The outside flange has access holes to allow for a

cooling coil, a vacuum line, and a temperature port. The cooling coil will carry either cold

water or liquid nitrogen for the purpose of condensing the spent suppressant fluid for recycle. 

The vacuum line is connected to the top of the vessel for the purpose of relieving the excess 

pressure before the next run. The line is connected to a knock-back condenser to prevent a 

minimal amount of suppressant from escaping. Additional ports in the sides of the vessel 

provide connections for the piping system and a fast-response pressure transducer. 

In anticipation of potentially large reaction forces during fluid discharge, a rugged 

framework was designed to support the test facility. As shown in figure 3.7, an 8 inch Ibeam, 10 feet high, was chosen to give a rigid backing for the source vessel while being 

heavy enough to dampen any vibrations exhibited during the test runs. This rack I-beam is 

directly bolted to the building structure for rigidity. The rack I-beam has three welded 

shelves for use in supporting the source vessel. The first is positioned just below the middle 

of the beam and is used to hold the source vessel from the bottom. Since this area

experiences most of the forces during the discharge tests, the shelf in this area is bolted to the

large flange at the bottom of the source vessel. This shelf is equipped with support legs for 

added stability. The next shelf is located near the top of the rack and anchors the top section 

of the source vessel by clamping to the 2-inch pipe. The last shelf, positioned lower down on 

the beam, is used to hold the mixing pump. The I-beam is also equipped with boltholes

along its face in order to add support structures for the level gauges and mixing lines. As for

the piping system, it is supported by a metal framework, which has provisions for clamping 

the pipe and preventing movement in both horizontal and vertical directions. The collection 

vessel is attached to a movable dolly to allow for opening of the piping system when access

is necessary. To account for substantial forces associated with the fluid discharging into the 

vessel, the collection vessel is held in place by support braces attached to a building I-beam. 

A picture of this entire system is shown in figure 3.8, where the source vessel, straight piping 

system, collection vessel, support structure and instrumentation can be seen. 

23
3.1.1 Pressure Sensors 

The test facility is equipped to measure instantaneous pressure readings in the source

vessel, along the piping system, and in the collection vessel by means of fast response 

pressure transducers. Two types of transducers were used in the experiments conducted, the 

Entran Model EPX-V01 with a range of 0-1000 PSIG (0-68 atmospheres) and the Validyne 

Model DP15 with varying ranges from all within 0-1000 PSIG. The Entran transducers were 

primarily used in the early runs when temperature drop during operation was not a concern. 

While the Validyne transducers were used for both pressure drop measurements and absolute

pressure measurements. The Validyne transducers were used with a transmission fluid for 

speed of transmitting the pressure wave in the pressure drop readings and also for insulation

when the temperature drop was considerable as in the HFC-125 tests. The Entran 

transducers have a diaphragm, which deflects according to the difference of the pressures, 

experienced at its front and backsides and provides an output through a strain gauge bridge 

located on the backside of the diaphragm. While the Validyne transducers used a

magnetically permeable diaphragm which when deflected caused a magnetic reluctance,

which was then interpreted by the inductance value of coils in the transducers body. In order

to install these transducers at measurement locations, special fittings had to be designed and 

fabricated. As explained previously in section 3.1, it is desirable to locate the transducers as 

flush to the inside surface of the pipe as possible, in order to minimize disturbances to the 

flow. This was accomplished with the specially designed fittings shown in figure 3.6. With 

this fitting welded into place on the pipe, the transducer can be screwed into the tapped hole,

allowing an O-ring to seal against the high fluid pressures. Bench top experiments were 

conducted to determine whether the transducers had appropriate response times for the fast

transient tests conducted. Both transducers were found to be at an acceptable level of 140 

milliseconds for pressurization and 190 milliseconds for depressurization when operating at 

20-30 atmospheres. A detailed analysis of these tests was provided in the first annual report. 

Although, the best test for the instrumentation would be to view their performance in actual 

testing situations. Figure 3.9 shows a typical output (run #5) from the pressure transducers 

installed in the test facility. This run used HFC-227ea as a test fluid and was initially 

charged to a source vessel pressure of 4.22 MPa. As one can see in the figure the pressures 

in the source vessel, pipe and collection vessel remained constant until the moment when the

valve was actuated (time = 0 milliseconds). Then the pressure in the source vessel dropped 

quickly while the pressures in the rest of the facility rose. This occurred until the pipe was 

filled with fluid (time = ~200 msec) and then there was a steady decrease in all the pressures

until all of the fluid was discharged from the source vessel (time = ~1500 msec). At this 

point a small pressure increase or “bump” was prevalent in the piping pressure traces due to 

the exit of the remainder gas contents of the source vessel. We also carried out a set of runs 

with similar initial conditions to check the repeatability of the measurements and the 

reliability of the instrumentation. The operating conditions are shown in the Table 3.1.1.1. 

below. 

Table 3.1.1.1: Operating Condition For Repeatability Tests.

Run Fill Volume Ambient Temperature Source Vessel Pressure Downstream Pressure 

#3 2358 32 4180 720

#4 2370 29 4195 525

#5 2388 31 4220 720

24
All the tests used HFC-227ea, as a test fluid the repeatability was very good as shown by 

figure 3.10, which compares the source vessel pressure histogram of the three experimental

runs. The slight discrepancies are most likely caused in the differences in the initial and

environmental conditions and not due to the instrumentation. With the example of the output 

from the entire facility and the fact that we had good repeatability we feel satisfied that the 

transducers did a good job of catching the pressure fluctuations during our experimental runs. 

There was also a necessity to measure the pressure drop across particular sections of

the piping system. While the absolute measurements at different locations allowed us to

measure the variation of absolute pressure across the discharge pipe, there was also an

interest to measure pressure drop across fittings in order to determine frictional losses across

these piping elements. Figures 3.11 and 3.12 show the two other piping configurations we

tested in this study besides the straight pipe configuration. Validyne pressure transducers 

were connected on either side of the fittings to measure the distinct pressure drop associated

with these fittings. This can be seen by the ports specifically designated to the fittings shown 

on the diagrams (i.e. PT2-PT3 across a 90° elbow). A typical output from this measurement 

is shown in figure 3.13, where pressure drop versus time is displayed. This plot shows the

pressure drop across a 90° elbow from run #9 using HFC-227ea as a test fluid and an initial 

pressure of 3915 kPa. There is a notable pressure drop instantly as the fluid front passes 

through the elbow, which deteriorates as time passes. This continues until approximately 

2300 milliseconds when the last liquid is discharged and the remaining gas leaves the system, 

denoted by the sudden drop in pressure attributed to a reduction in the frictional component 

of pressure drop associated with gas flow over liquid flow. Similar plots were made for the 

pressure drop across capped tees, through tees and unions. Also, a comparison was done

between the vertical branches of a through tee to study the effect of phasic composition on 

the separation of the flow. With these measurements, it was possible to successfully analyze 

the effects of different fittings on the pressure losses in the discharge piping. 

3.1.2 Flow Rate Sensor 

As mentioned previously, a key parameter that we wish to measure is the

instantaneous mass flow rate of fluid during transient discharge from the source vessel. To 

the best of our knowledge, this has not been successfully measured in prior experiments due 

to its inherent difficulty. If fluid inventory in the source vessel could be measured as a

function of time, then the rate of discharge would be given by the gradient of inventory 

versus time. Because of possible void generation in the liquid (due to gas evolution) and 

associated level swelling, observation of liquid level does not provide indication of fluid 

inventory. Our approach is to measure the static head of fluid in the source vessel, as a 

dynamic function of time. Since frictional and kinetic terms are negligible in the overall 

momentum balance for conditions in the vessel, we expect to obtain the fluid inventory in the 

vessel directly from static head measurements: 

 P

g

A M = • ∆ (3.1.2-1) 

where M is the fluid inventory, A is the cross-sectional area of the vessel, g is gravitational 

acceleration, and ∆P is the pressure difference from the bottom to the top of the vessel (static

25
head). The instantaneous rate of mass flow out of the vessel is then obtained from the time

derivative, 

 

 dt

d P

g

A

dt

dM

m

∆ = = • •

 (3.1.2-2)

 

To obtain the desired accuracy for the rate of mass flow, a precise measurement of ∆P 

is necessary. We used a Validyne Model DP15 differential pressure transducer to obtain this

static head ∆P. As shown in figures 3.2 and 3.14, this transducer was connected to the vapor 

space and the liquid volume at top and bottom of the source vessel, respectively. The 

transmission line connecting the vapor space to the liquid space was filled with a Meriam

Red fluid with a specific gravity of 2.95 to maximize response time. The theory behind this 

was that pressure signals will be transmitted to the transducer at the speed of sound in liquid, 

sufficiently fast to obtain dynamic data during suppressant discharge and much quicker than 

through a compressible gas volume. Therefore, the connections to the liquid and vapor

spaces were made at the wall of the Marotta valve and the wall of the source vessel, 

respectively. This allowed for an almost flush-mounted differential pressure transducer with 

an excellent response time for measuring liquid inventory. Figure 3.14 shows a diagram of 

the set up for this differential pressure transducer. Notice that we attempted to fill the 

transmission line as close to the source vessel as possible, therefore reducing any pockets of 

fire suppressant vapor and making a more accurate measurement.

We were pleased to find that our thinking worked very well, as shown in figure 3.15, 

a plot of the cumulative mass leaving the source vessel versus time. This example shows the 

same experimental run described in the previous section (3.1.1). It shows that at the point 

when the valve was actuated a quick release of test fluid left the source vessel and continued 

gradually until the source vessel was empty, at approximately 3.37 kg. This corresponded 

with the pressure drop in the source vessel shown in figure 3.9 and with the end of the liquid 

flow and beginning of the gas flow at a time of ~1500 milliseconds. This is what we would 

expect since we did not have much fluctuation in the pressure histogram leading us to believe 

that there was a steady flow of fluid through the facility. Once again, we had excellent 

repeatability and felt that this method for measuring the instantaneous liquid inventory gave 

an acceptable representation. 

3.1.3 Temperature Measurement Sensors 

Fluid temperatures along the discharge pipe are to be measured using film 

thermocouples. It has been experimentally shown at Lehigh, Cokmez-Tuzla, Tuzla, Chen 

(1993), that this type of thermocouple can provide response times in the order of 

milliseconds. A schematic of the film thermocouple is shown in Figure 3.16. Since the hot 

junction of the thermocouple is comprised of the metallic film, it is essential to maintain a

continuous electric conductive-path through that film. Two issues are involved: compatibility 

of the film material to the test fluid, and the erosive effects of the flow on the film. During 

the past experiments at Lehigh, a metallo-organic paint was used to make the film at the tip 

of the thermocouple. This film worked well in experiments with high temperature steam and 

water. After a search for alternatives, we decided once again to use this same material for the 

junction film. A bench setup was prepared to test the thermocouple with this film junction. 

26
After forming the film, the thermocouple assembly was exposed to test liquids for duration of

up to five minutes. Tests with the primary test fluid, HFC-227ea, indicate that the junction 

film holds up well to such exposure. No bench test for erosion of the film was attempted 

prior to the running of the test facility, due to the significant effort involved. We found that 

the film held up for anywhere from 5-10 experimental runs before deterioration in the film

deemed it necessary to dismantle the thermocouple from the facility and reapply the metallic

film. We determined the need for reapplication of the metallic film by an increased 

fluctuation in the output signal of the thermocouple or in some cases a loss in the signal 

altogether. Many of the failures in the thermocouple were not attributed to the metallic film

at all, but to the disintegration of the inner wiring of the thermocouple. This was most likely 

due to the vibration encountered in the test facility when the experimental runs were 

conducted and the tight seal on the thermocouple casing needed to hold the thermocouple in 

place during testing. 

An example of a typical output from the four film thermocouples located throughout

the test facility is shown in figure 3.17. The top set of data is the signal from the film

thermocouple located in the vapor space of the source vessel. The next three sets are from

the film thermocouples in pipe positions #3, #2, and #1, respectively from top to bottom.

This plot shows no significant temperature change in the vapor space, which is due to the fact

that there was only a gas expansion involved in this region. While an adiabatic gas 

expansion from a volume of 1.5 L to 34 L should have resulted in a large decrease in 

temperature, in our case this was not an adiabatic expansion due to the availability of a large

heat sink in the stainless steel walls. With the considerable warming from the facility walls

and the thermocouple’s location at the top most portion of the entire facility, the only 

temperature drop experienced would be minimal. As for the other set of thermocouple 

traces, we believe the most significant reason for the difference in temperature drop as we 

move along the pipe is also due to the warming effects. 

There were also two standard, 1/32” diameter, shielded thermocouples placed at

different locations in the facility, one each in the source and collection vessels. Figure 3.18 

shows a standard temperature trace from these thermocouples, again from run #5. The two 

traces shown in this diagram are significantly different from each other. The collection

vessel temperature does a quick dip and then warms back up to room temperature. This is 

due to the filling of the vessel with cold fluid, which is probably still going through the 

process of dissolution of nitrogen from liquid suppressant. Also, there is an energy loss 

associated with the phase change of the suppressant fluid as it goes from liquid in the pipe to 

gas in the collection vessel. This goes on until the collection vessel reaches the saturation

pressure of the suppressant fluid and then the vessel begins to heat up due to compression of 

the gas in the volume. On the other hand the thermocouple located at the bottom of the 

source vessel initially does not change much since it is submersed in liquid suppressant. 

Then as the liquid in the source vessel runs out, it experiences the vaporization of the 

suppressant fluid on its surface and the expansion of the remaining gas in the system. To 

fully understand the complete picture of the changes of phase taking place in the 

experimental runs, we need a representation of the void in the facility, which is described

below.

27
3.1.4 Void Fraction Sensor 

One of the major flow variables for two-phase flow is the void fraction. Many flow 

parameters, such as thermodynamic quality, vapor and liquid mass fluxes, and transport 

properties of vapor and liquid affect the void fraction. In turn, the void fraction strongly 

affects pressure drop, holdup inventory, and flow rate for any given operating condition.

Thus, it is desirable to experimentally measure the void fraction during the transient 

discharge of suppressant. Once again, this is a difficult measurement in two-phase flows and 

was not attempted in prior experiments of suppressant flow. In the present work we 

attempted to measure the transient void fraction at the end of the discharge pipe. The 

desirable features of a sensor for void fraction are: 

• No disturbance on flow 

• Integrate the void fraction over the cross-section 

• Fast response 

Based on previous experience at Lehigh University, we selected to use a capacitance 

sensing method for this measurement, in expectation that this method will satisfy all of the 

above desirable features. 

A number of different capacitance probes were designed, fabricated, and bench 

tested. The final probe configuration is shown schematically in figure 3.19. As indicated the 

probe consists of a positive electrode and a ground electrode, both placed on the outside

surface of the discharge pipe. The section of the discharge pipe at this location needs to be 

electrically nonconductive to allow capacitance measurements. For the purpose of these 

bench tests, a glass tube is being used as the pipe material. The positive electrode is a strip of

copper in contact with the outer surface of the pipe. The electrode span is as wide as the ID 

of the pipe. The ground electrode is placed on the opposing side of the pipe. When powered 

by a high-frequency voltage, this probe would measure an overall capacitance between the

two electrodes, which includes the whole cross-section of the tube. This capacitance can

then be calibrated to give the void fraction of the vapor/liquid mixture in the pipe. There is a

guard electrode around the positive electrode to prevent/reduce stray capacitance. 

This probe was tested on a bench loop with water and HFC-134a. The schematic of 

the bench setup is shown in figure 3.20. The test started with the glass tube filled with water. 

During the test, the tube was emptied from the bottom so that the liquid level receded with a 

constant speed. The recorded signal from the capacitance sensor was compiled as a voltage

signal from the capacitance meter and then reduced into void fraction with a calibration using 

the signals associated with all liquid and all vapor. An example of the output from this test is

shown in figure 3.21. It is seen that initially the probe senses 100% liquid in the tube. As the 

liquid level approaches the top edge of the sensor, the signal slowly starts to show some void 

in the sensing volume (edge effect). When the liquid level is within the sensing volume, the 

variation of measured void fraction is linear with liquid level. Again, there is an edge effect

when the liquid level is close to the bottom edge of the sensing volume. Later, when the 

liquid level falls clear of the sensing volume, the probe senses all air in the tube. Despite the

fact that we guard the sensing area, there are some edge effects in the axial direction. The 

above test was repeated with HFC-134a under pressure and similar results were observed. 

Our biggest problem with changing this design to fit the conditions experienced in our 

28
experimental test facility would be containing the higher pressures. In the previous bench 

examples we dealt with only low-pressure systems and in the present application the system

pressure is as high as 1250 kPa at the location of the void probe. The other difficulty was 

that we must not have an electrically conductive surface between the probe and the tested 

medium. Due to these limitations we developed a design that would be constructed on the 

inner surface of the pipe as shown in figure 3.22. There were three main pieces to this 

design; the first was the electrical plug welded to the metal pipe, which served as a pressure 

seal and as an electrical connection to the capacitance meter. The next two pieces were the

guard and the probe separated by a small layer of insulation. To insure that these pieces were 

not washed away in the flow during the runs, they were covered in a thin layer of epoxy. 

This capacitance probe used the pipe wall as its ground creating an electrical field similar to 

the one shown in figure 3.23. While the orientation of this design was horizontal as 

compared to vertical in the bench test setup, we believed that it would give us similar results. 

We conducted an insitu calibration by filling the area with all liquid suppressant and then 

letting it run out so that only suppressant vapor remained. An example of the final output of 

an actual experimental run is provided in figure 3.24 where void fraction versus time of run 

#5 is shown. The signal initially starts out at a void fraction of 1.0 (all vapor) and then 

shoots down to approximately 0.4 as the area is flooded with the transient flow of two-phase 

fluid. The signal jumps around as the mixture of the two phases changes depending on the 

speed and location of the dissolution of the nitrogen from the liquid suppressant. As the flow 

slows down, so does the fluctuation in the capacitance signal, until the liquid runs out (time = 

1700 msec) and the remaining gas is vented into the collection vessel causing the signal to 

abruptly go back to the 1.0 region. Our largest problem with the capacitance signal resulted 

from the static charge build up of the fast moving dispersed flow. It caused the signal to

jump out of range and required a significant filtering effort to be conducted. This is also the 

reason for the offset in the final value at the end of the run. Despite these slight problems, 

we feel that the capacitance probe design gave an excellent representation of the void fraction 

values experienced during the experimental runs.

3.2 Test Matrix Performed In Study 

After the completion of the test facility, we needed to choose the appropriate initial 

conditions in which to conduct the experimental runs. Our main goal was to encompass a 

wide variety of test parameters while keeping the total number of runs at a minimum to 

lessen the wear on our instrumentation. Table 3.1 shows the experimental runs conducted in 

this study. We used HFC-227ea, also known as HFC-227ea as our primary test fluid, since it 

was recommended at the time as the alternative fire suppressant with the best fire retarding 

ability. We also used HFC-125 or R125 since it had been studied recently and there was 

some test data that we could use for future comparison. In the preliminary experiments, 

water and carbon dioxide were used to test the new instrumentation after installation. For

example, run #19 used water as a test fluid and was undertaken to test the capacitance sensor 

and film thermocouple without exposing them to high pressure or a corrosive fluid. Run #18 

also used water and had the intent of testing the mass flow sensor response and all of the

thermocouples in the facility. In a similar manner, runs #20 & #21, which used CO2, were 

used to measure the response of the film thermocouples to a sudden drop of temperature as 

low as -40° Celsius, while also allowing for us to gauge the error encountered by the pressure 

transducers at this temperature. These four runs utilized the straight pipe configuration in the 

test facility shown previously in figure 3.5. The HFC-125 runs (#14-#17) were conducted 

29
with only one objective, to test the influence of changing the source vessel pressure. The 

four runs were carried out at source vessel pressures varying from 2050 to 5000 kPa. It must

also be noted that any runs using source vessel pressures with initial conditions using low

pressures were conducted at fill levels between 85-90%. This was done in an attempt to 

facilitate as much vaporization during the flow discharge as possible since the downstream 

pressure would remain at a lower pressure with the smaller amount of nitrogen gas available 

to expand. 

The majority of our experimental runs used HFC-227ea as a test fluid and were

broken into three main piping configuration groups. First, there were the tests using the 

straight pipe configuration. These tests had three different goals: (1) to test the repeatability 

of the instrumentation, (2) to observe the effect of changing the source vessel pressure, and 

(3) to determine the effect of a heat source or “hotspot” along the piping system. The

repeatability runs (#3-#5) were conducted at an initial source vessel pressure of about 42 

atmospheres. The result of the source vessel pressure histogram is shown in figure 3.10. 

While the “hotspot” run (#7), was conducted at a similar starting pressure with a foot long 

section at 55° Celsius between the 2nd and 3rd instrumentation ports on the pipe. Only one 

such run was needed since it was the conclusion of this group that a minimal effect resulted 

from the heat addition under such rapid transient conditions. The third goal was to observe 

the effect of varying the initial source vessel pressure from 2125 to 4910 kPa (#1-#6). The 

second main piping configuration group used the piping configuration shown in figure 3.11, 

where a vertical “U” turn was built into the piping system. It allowed us to measure the effect 

of pipe fittings on the pressure loss in the system. The pressure drops across a capped tee, 

90° elbow and union were measured with varying initial source vessel pressures from 2120-

4910 kPa in runs #8-#10. The third piping configuration group used a complete vertical split 

as shown in figure 3.12, where the tee was no longer capped but used to join two symmetric 

piping sections. The two sections were identical to the one used in the previous grouping 

although this time the pressure drop across the split was measured. A sensor was placed

before the tee and after the 90° elbows on either arm of the split to test whether there was an 

effect from the splitting of the flow. These test were also conducted at three varying initial 

source vessel pressures from 2140-4910 kPa in runs #11-#13. These tests gave us a good 

look at the effects of changing the initial source vessel pressures, changing the test fluids, 

changing the piping configurations and the repeatability of our instrumentation. 

30
Table 3.2: Matrix of Experimental Runs Conducted at the Lehigh Test Facility 

Run # Test Fluid Piping Configuration Fill Volume [mL] Psvo [kPa] Pcvo [kPa] 

1 HFC-227ea Straight Pipe 3571 2125 10 2 HFC-227ea Straight Pipe 2340 3235 635 3 HFC-227ea Straight Pipe 2358 4180 720 4 HFC-227ea Straight Pipe 2370 4195 525 5 HFC-227ea Straight Pipe 2388 4220 720 6 HFC-227ea Straight Pipe 2150 4910 715 7 HFC-227ea Straight Pipe

w/ Hotspot @ 55°C 

2200 4230 710 

8 HFC-227ea Capped Tee (Config. #1) 3420 2120 125 9 HFC-227ea Capped Tee (Config. #1) 2500 3915 575 

10 HFC-227ea Capped Tee (Config. #1) 2450 4910 615 11 HFC-227ea Through Tee (Config. #2) 3507 2140 600 12 HFC-227ea Through Tee (Config. #2) 2400 3935 655 13 HFC-227ea Through Tee (Config. #2) 2460 4910 615 14 HFC-125 Straight Pipe 2900 2050 110 15 HFC-125 Straight Pipe 2363 3635 175 16 HFC-125 Straight Pipe 2400 4245 165 17 HFC-125 Straight Pipe 2467 5000 125 18 Water Straight Pipe 3585 2875 120 19 Water Straight Pipe 3100 4150 100 20 CO2 Straight Pipe ~1000 5215 100 21 CO2 Straight Pipe ~1000 5360 100 

31 
3.3 A Typical Experimental Run 

Since we have explained our objectives behind the building and operating of the test 

facility, the operating procedure and output should be explained in more detail. First we set 

the plunger in the Marotta valve into the closed position using the arming tool. Once the 

valve is closed, the next task is to recycle any fluid remaining in the collection vessel back to 

the source vessel using a magnetically driven centrifugal pump. If there is not enough test 

fluid in the collection vessel to reach the desired initial fill level for the next run, then the 

recycle line is closed and a line to a fresh-supply tank above the test facility is opened, 

allowing the remaining liquid needed for the run to be pumped into the source vessel. Once 

this is accomplished, the supply line is closed and a circulation line is opened allowing the 

liquid at the bottom of the source vessel to be pumped to a spray nozzle at the top of the 

vessel. The spray nozzle then creates a fine mist of liquid to float through the vapor space of 

the vessel before joining the rest of the liquid. At this point we begin to slowly add nitrogen 

gas to the source vessel through a connection at the top of the vessel. We continue to run the 

circulation pump until the liquid is saturated with nitrogen. This is recognized by the constant 

pressure in the source vessel, which usually took about 120 minutes. The circulation pump is

then turned off and the contents of the source vessel are allowed to settle and cool down to 

room temperature. The data acquisition system is initiated and run for 500 milliseconds

before the Marotta valve is actuated. This allows for the collection of the initial conditions of

the system. Once the valve is actuated, the contents of the source vessel discharge through 

the Marotta valve into the piping system and finally collect in the collection vessel. To fully 

understand the output from a single run, a typical experimental run (#5) will be examined in

more detail. The pressure histograms can be seen in Figure 3.9 where the traces go from the 

source vessel at the top to the collection vessel at the bottom and the pipe ports in the middle. 

In specific, the black line (diamond) is the source vessel trace and the green (circle), red 

(square) and blue (triangle) traces are the first, second and third pipe positions, respectively, 

with the orange (cross) line being the collection vessel. Initially the source vessel is at 4220 

kPa and the pressure downstream of the closed valve is at 720 kPa. After the valve is 

opened, the pressure in the source vessel quickly drops as its contents fill the piping system, 

explaining the sudden rise in pressure there. At the point when the pipe is completely filled, 

~200 milliseconds, there is a consistent drop in pressure as the fluid empties into the 

collection vessel. At ~1600 milliseconds the fluid begins to run out and the flow in the 

system is predominately vapor. This quickly causes the entire system to equilibrate to the

same ending pressure, ~1100 kPa. 

Not only can we follow this discharge with pressure histograms, but we can also 

follow the temperature change throughout the system. Figure 3.17 shows the output of the 

highly sensitive film thermocouples placed at the top of the source vessel in the vapor space 

and along the bottom of the piping system at the three port locations. The same color and 

symbol scheme is used as in the pressure histogram with the source vessel being on top, 

followed by the third, second, and first pipe positions, respectively. There was little 

temperature change experienced in the fluid flow portion of the discharge since there was 

little flashing of the suppressant. The main temperature drop occurs at the point when the 

vapor discharge begins (~1600 ms). This is due to the sudden expansion of the remaining 

nitrogen gas and falls as low as 10° Celsius at the first pipe position. There are lesser

temperature drops as it travel farther down the pipe due to the increased distance from the 

initial point of expansion and the warming of the vapor by the system walls. There were also

32 
two 1/32” shielded thermocouples placed in both the source vessel and collection vessel. 

The outputs of these are shown in Figure 3.18, where the black (diamond) symbols are for 

the bottom, fluid space of the source vessel, and the orange (cross) symbols are for the top,

vapor space of the collection vessel. There is a small dip in the temperature of the collection 

vessel as the fluid initially reaches the large volume of the vessel and then it warms back up

to room temperature rather quickly. As for the source vessel temperature, this thermocouple 

experiences the temperature drop due to the expansion of the gas that the film thermocouple 

in the vapor space did not. This is most likely due to the fact that the film thermocouple at 

the top of the vessel is in contact with stagnant gas and did not have enough time to sense its 

temperature while the lower thermocouple experienced the shock of the quick expansion 

from the fast flowing gas. 

The next insight into this study was provided by the differential pressure transducer, 

which measured the mass flow rate of the discharge. As seen in figure 3.15, this

instrumentation allowed for the measurement of the cumulative mass of liquid leaving the

system. This gives us a good idea of how much liquid has left the source vessel and might

possibly land on the desired destination, in this case, into the collection vessel. Ideally, we

would have liked to locate such a sensor near the desired destination, but since we could not 

develop an accurate, quick response method of measuring the collected fluid, this seemed the 

most suitable alternative. Figure 3.15 verifies what we have already observed from the

previous measurements of pressure and temperature. At first there is a primarily liquid flow, 

starting rapidly and slowly dying out as the inventory is relinquished to the piping system, 

until no liquid is left at ~1600 ms. Along with this additional measurement to backup our 

other assessments, we also measured the void fraction leaving the pipe by use of a

capacitance sensor located at the end of the pipe. This instrument obtained an electrical 

signal, which pertained to a capacitance measurement of the flow passing through this 

section. This signal was then converted into a measurement of void fraction through a

previous calibration where a value of 0.0 denotes all liquid flow and a value of 1.0 denotes 

all vapor flow with everything in the middle being a mixture of the two. Figure 3.24 shows 

the output from this sensor during the experimental run. Once again we notice that the main 

liquid flow seems to end in the 1600-1800 millisecond range. The slight discrepancy

between this measurement and the mass inventory measurement is most likely caused by the 

fact that they are on either ends of the pipe. This trace shows us that at the beginning there is

a flow of predominately liquid, void fraction of 0.4-0.2, which slowly rises to a more

consistent value of 0.5 until the liquid runs out and the value jumps back up to an all vapor 

level. There are many peaks in the trace due to the fact that the flow is highly turbulent at 

this period in the flow, with the phasic composition changing with the rate of nitrogen 

coming out of solution. The flow slows down a little near the end and the trace smooth out, 

until the end when there is no more flow. These four different methods of taking 

measurements allows us to believe that we have a good grasp on what is going on during the 

rapid transient phase of the experimental run. For the detailed output of all of the

experimental runs conducted in this study see Appendix B. 

3.4 Observations and Conclusions 

After analyzing the full range of experiments conducted in this study we found many 

findings and observations. The most exciting developments in this study pertained to the 

previously unsuccessful measurements of mass flow rate, void fraction and local 

33
temperature. Our transient fluid inventory method gave excellent results providing us with 

not only another check for the computer code developed in this work, but also an innovative 

technique currently unavailable in two-phase flow experiments. The void fraction 

measurement using a capacitance sensor has been used in other applications, but never in 

such a rapid-transient setting. In addition, most methods used in previous studies required

the probe to extend out into the measured flow, which causes obstruction to the flow and 

obvious changes in the resulting flow characteristics. Our method provides a quick, accurate 

signal with minimal disturbance to the flow, resulting in the measurement of more realistic 

flow characteristics. Last, the film thermocouples used in this study also allowed for the

measurements to be made without disturbing the fluid flow. In previous studies 

thermocouples were always extended into the middle of the pipe in order to get the most 

accurate representation of the local temperatures causing them to break and deflect the fluid

as it flowed. These techniques had another shortcoming of not allowing a fast response due

to the thermocouple sheathing. In our case, we have only a paper-thin layer shielding the 

thermocouple junction from the flow while having no effect on the flow due to their flushmounted design. These techniques allowed this study to provide a more comprehensive 

comparison for the computer code simulations to be tested against. 

The experiments in this study primarily focused on the use of HFC-227ea as a test 

fluid. The first set of experimental runs examined the effect of changing the initial source

vessel pressure on the flow characteristics through a straight, 1/2” nominal diameter pipe, see 

Figure 3.5. The main difference between these runs was the speed at which the source vessel 

discharged its contents into the piping system and hence into the collection vessel. All the

runs exhibited the same basic phenomena as seen in Figure 3.9; (1) a quick pressure drop as 

the pipe filled, (2) a continuous drop in the pressure as the source vessel discharged its 

remaining fluid, and (3) a slight pressure “bump” in the pressure measurements in the pipe 

due to the venting of the remaining gas and vapor in the source vessel. Similar observations 

can be made by examining the void fraction and mass flow measurements. The lowest void 

fraction can be seen at the beginning of the runs, Figure 3.24, when the liquid content in the

flow is the greatest. Near the end of the run, we see a marked increase in the void fraction as 

more nitrogen comes out of solution causing it to rise. The corresponding mass flow 

measurement, Figure 3.15, shows the fastest discharge at the beginning of the run when the 

source vessel pressure is at its highest value providing the greatest driving force. Once the 

pressure begins to wane, so does the flow rate shown in the measurement, until it ends at the

end of the liquid flow for the run. 

The second set of experiments were conducted using an alternative piping system

which contained tees, elbows, and a union as shown in Figure 3.11. The main goal for these 

runs were to measure the pressure drop associated by adding specific pipe fittings into the 

piping network. As we would have expected, the flow was significantly reduced by the pipe

fittings, as shown in Figure 3.25. The comparison shows that the run (#1) with the straight 

pipe finishes discharging its contents at approximately 4600 milliseconds compared with the 

run (#8) that used the piping layout containing two capped tees, a 90° elbow, a union and two 

45° elbows, which finished at about 7500 milliseconds. To study these pressure drop

measurements, a comparison to the literature values for non-compressible, turbulent flow

through pipe fittings using the friction loss factor method was conducted. This method is

based on the procedure where pipe fittings are set equal to equivalent lengths of pipe 

depending on the anticipated pressure drop associated with the fitting. In our experimental

34
setup, the pressure ports used to measure the pressure drop across the fittings included not 

only the fitting but also some piping in between, see figure 3.11. In each case the additional 

piping plus the fitting between the pressure ports was included to remain as accurate as 

possible. For example the 90° elbows was made up of the pressure port then 2 inches of 

piping, then the elbow, then 3 more inches of piping and the final pressure port. Therefore, 

what we denote as the “90° elbow” is actually the elbow plus 5 inches of piping. The 

literature values were obtained by using the literature value for the friction loss factor, K, for

the particular fitting and adding it to the calculated value for the additional piping, found by 

using equation 3.4-1. 

 D

L K f = 4 f (3.4-1) 

Where ff is the friction factor of the pipe found by using a chart in the literature with the 

diameter and estimated roughness of the pipe, L and D are the length and diameter of the 

pipe, respectively. Once a total K is found for the fitting and additional piping it can be used 

in equation 3.4-2 to find the literature value for ∆P. 

 2

2 K v P ρ ∆ = (3.4-2) 

Where ρ is the mixture density calculated using the experimental void fraction 

measurements, and v is the flow mixture velocity calculated from the experimental values of 

mass flow through the pipe and the mixture density. 

Table 3.4.1: Comparison of Experimental Pressure Drop Values to Literature Values*. 

Run #8 Run #9 Run #10 

Initial Source Vessel Pressure [kPa] 2120 3915 4910 

Void Fraction 0.649 0.507 0.485

Liquid Density [kg/m3] 1403.5 1392.8 1389.0

Mixture Density [kg/m3] 514.43 704.73 733.18

Velocity [m/sec] 3.446 9.007 7.630

90o Elbow 

Experimental ∆P [kPa] 9900 27800 36000

Literature Value ∆P (mixture ρ) [kPa] 21594 (118%) 106073 (282%) 72775 (102%)

Literature Value ∆P (liquid ρ) [kPa] 7915 (20%) 53673 (93%) 38414 (7%)

Capped Tee 

Experimental ∆P [kPa] 11300 40400 62600

Literature Value ∆P (mixture ρ) [kPa] 29778 (164%) 195398 (384%) 134059 (114%)

Literature Value ∆P (liquid ρ) [kPa] 14580 (29%) 98871 (145%) 70763 (13%)

Union 

Experimental ∆P [kPa] 11100 31500 50400

Literature Value ∆P (mixture ρ) [kPa] 17048 (54%) 83742 (166%) 57454 (14%)

Literature Value ∆P (liquid ρ) [kPa] 6248 (44%) 42373 (35%) 30327 (40%)

* From Fundamentals of Momentum, Heat, and Mass Transfer by Welty, Wicks, and Wilson. 

35
Since our experiments took place under two-phase conditions, we used the mixture density in

order to get as close to the actual conditions in the pipe as possible. The one drawback of 

this method is that the void fraction sensor was not located at the local position of the fitting, 

but at the end of the pipe so it may not be a precise representation of the void conditions at 

the fitting. For this reason and because the literature values for K were developed using liquid 

turbulent flow, we also included a comparison to values using the liquid density and liquid 

velocity in place of the mixture density and mixture velocity. See table 3.4.1 for the detailed

values of this comparison. As the comparison shows, there are some discrepancies between 

the literature and the experimental values, with the calculations using liquid density being in 

reasonable agreement. In order to compare whether the mixture or the liquid density is best 

used when calculating the literature value for ∆P, the error has been provided in parentheses. 

This error was determined using the experimental value as the actual value. The 

experimental values for pressure, void fraction, density and velocity were taken in the middle 

of the runs, with run #8 being from 1500 to 3500 milliseconds and runs #9 & #10 being from 

500 to 1500 milliseconds. Some of the possible reasons for the discrepancies in the 

comparison between the calculations using the liquid density and those using the mixture 

density may be: 

1) Because of the phasic composition of the flow at the particular time during the run. For

example, at the beginning of the run there is a large amount of liquid flowing but there is 

also a large quantity of nitrogen degassing, so it is hard to tell how this will effect our 

calculations.

2) The capped tee is near the beginning of the piping layout while the void sensor that we

are using for the calculation of the mixture density is at the end. 

3) The fact that the literature value for the union is actually considered as if the fitting is just 

a section of straight pipe. In theory this should be the case since the union when 

assembled should not interfere with the flow, but there are still ridges from the male and

female sections connecting on the inside of the union possibly causing disturbances to the 

flow. 

4) The role of the two-phase flow that we have in the experiments. An assumption made

when modeling two-phase flow through pipe fittings is that the differences with singlephase liquid flow are marginal, but this may not be the case. The void fraction may play

a large part in the pressure drop experienced through a pipe fitting since it will influence

the flow conditions through the fitting. 

Overall, the liquid density yielded a better calculation than the mixture assumption with an 

average error of 47.3% compared with 155.3%. These findings lead us to believe that the 

phasic composition of the flow through the fitting plays a major role in the pressure drop

experienced at any instantaneous time during the discharge and local flow characteristics will 

play a large role in what magnitude of pressure drop occurs. 

The third set of experiments also used an alternative-piping configuration, which 

included a vertical split of the flow, see Figure 3.12. The goal of these runs was to examine 

whether there would be any difference in pressure drop between the two legs of the split. At 

slower flow velocities, separation of the flow may occur based on the flow regime of the

discharge. For example, if the flow were stratified then all the liquid may flow through the 

bottom leg of the split and the entire vapor may travel through the top leg. This would be 

36
evident from the difference in pressure drop between the two legs since the vapor would have 

a much lower pressure drop than the liquid. In our experimental runs we examined no such 

trend, as can be seen in figure 3.26, which depicts the pressure drop across both legs after a 

vertical tee split. The plots for both the top and bottom branch are almost identical which 

was the case for the entire tested pressure range of 2140 kPa to 4910 kPa. This proves our 

theory that at these high velocities the flow is in the bubbly flow regime causing an even split 

between the two branches. 

The fourth set of experiments used the same piping configuration as the first set, the

straight pipe, but we changed the test fluid to HFC-125. The main difference observed in 

these runs over the HFC-227ea runs was the drop in temperature as the liquid vaporized. If a 

comparison is made between run #5 (HFC-227ea) and run #16 (HFC-125), which have initial 

source vessel pressure within 5 kPa of each other, there is an undeniable difference in the 

temperature traces in the collection vessel. As Figure 3.27 shows, the temperature in run #5 

only drops about 3 degrees compared to a temperature drop of almost 35 degrees in run #16. 

While the other temperature traces between these runs are similar, this may be because the 

flow passes through the pipe at such a fast rate with a high convective heat transfer 

coefficient, which causes the fluid to warm up. However, in the collection vessel the

volume/surface ratio is much larger and convective heat transfer to surface area is smaller. 

Although, it would have been beneficial to measure the temperature at the end of the pipe to 

justify this theory we were unable to take this measurement due to a thermocouple failure. 

Once the liquid reaches the collection tank it must vaporize to reach the saturation conditions

in the collection vessel causing a drop in temperature. These conditions did not significantly 

effect the other measurements near the end of the pipe for either fluid. For instance, for 

either fluid the void fraction followed a similar trend as expected from the previous runs, it 

dropped off as the initial flow passed the sensor and then rose slowly until the liquid ran out, 

then it returned to the all vapor value. These plots can be examined in more detail in

Appendix B as the forth graphs under runs #5 and #16. During runs with both of these fluids 

we did not observe a noticeable pressure recovery due to the degassing of the nitrogen, which 

was seen in previous studies. Although, when examined closely, we see a sudden slow down 

in the pressure decay after the initial pressure drop, which may indicate a slight recovery. 

This theory is based on the fact that in the HFC-227ea runs there was a much smoother 

transition shown on the source vessel pressure histogram between the initial pressure drop 

and the filling of the pipe. This can be seen best by once again comparing runs #5 & #16, 

where there is a slight “bump” in the source vessel pressure trace at about 400 milliseconds

of run #16. The absence of a more pronounced pressure recovery due to the degassing of 

nitrogen can be accounted for by either of two reasons. The speed of the flow is too great for 

a noticeable increase in pressure from degassing to be observed at any one point or the 

degassing took place gradually over the run preventing one noticeable increase to be 

possible. We do know that the degassing occurred at some point since the void fraction 

validated that a two-phase flow existed and we are certain it was not due to vaporization 

since there was no temperature drop to make us believe otherwise. 

The last sets of run were conducted during preliminary runs to test the

instrumentation after installation in the system. The first test using water, run #18, allowed 

for a sampling of the mass flow rate sensor. As can be seen in figure 3.28, this test worked 

out very well and gave an excellent representation of the instantaneous mass inventory in the 

source vessel. This test also showed us that our thermocouples could pick up quick drops in 

37
temperature, as shown in the pipe position 1, trace of the film thermocouple which picked up 

the nitrogen expansion at the end of the run, see Appendix B, run #18, graph 2. The second 

test, run #19, gave us an excellent sample of the capabilities of the capacitance probe for 

measuring void fraction. Since there would be no vaporization or dissolution of gas in this

run, the void fraction should go from a value of 1.0 to 0.0 as the liquid flowed. This is 

exactly what was found at the beginning of the run until it stopped at precisely the same point 

as the temperature dropped signifying the start of the nitrogen expansion through the system,

as can be seen in Figure 3.29. The final set of experimental runs used carbon dioxide as a 

test fluid and were conducted with an open collection vessel to prevent a high pressure at the

end of the run. The goal of these tests was to get a grasp on the response time of the film

thermocouples and what range of temperature drop they could collect. The thermocouples 

were able to pick up a quick drop in temperature of up to –40° Celsius. Detailed graphs of

the measurements of all 21 experimental runs can be found in Appendix B. 

38
2

0

0

c

m(

6

.

5

6

'

)

Tsvs

Psv

Tsvu

4.93 cm (1.94") ID

DISCHARGE 

VALVE

SOURCE TANK

V = 3.8 Liters

Tcv

COLLECTION VESSEL

V = 30 Liters

PIPING SYSTEM

Pcv

P1

T1

P2

T2

P3

T3

CIRCULATION LINE

SUPPRESSANT RECYCLE LINE

Rapid Discharge Facility (R

Lehigh University

CAPACITANCE SENSOR

∆P

Figure 3.1: Schematic of Test Facility: Source Vessel, Discharge Piping, and Collection Ve

39
2

0

0

c

m

(

6

5

.

6

'

)

Tsvs

Psv

Tsvu

4.93 cm (1.94") ID

∆P

Figure 3.2: Schematic of Source Vessel Showing Mass Inventory Measurement. 

40
0

500

1000150020002500300035004000

-60 -40 -20 0 20 40 60 80 

Temperature [degrees Celsius]

Pressure [kPa]

HFC-125

HFC-227ea

Figure 3.3: Pressure versus Temperature Saturation Curve for HFC-227ea & HFC-125

41
1.25"

Inlet

Outlet

Electric 

Switch

1

7

.

5

"

Figure 3.4: Marotta Valve. 

42
10.4' (317 cm)

5.5" (14.0 cm)

P Port

7.5" (19.1 cm)

T Port

42.25" (107.3 cm)

44.25" (112.4 cm)

P Port

T Port

5.75

P 

Specially Designed Section for

Void Fraction Measurement

Source Vessel & Holding Rack Collection 

Figure 3.5: Top View of Test Facility, schematic of Discharge Piping System. 

43
0.15"

1/8"

11/32"

1/8" NPT

T Port Fitting

Thermocouple

Fitting

1/2" Pipe

11/32"

1/2" Pipe

P Port Fitting

Pressure

Transducer

PRESSURE TRANSDUCER PORT

THERMOCOUPLE PORT

Figure 3.6: Schematic of Pressure Fittings for Temperature and Pressure Ports. 

44
1

5

"

4

2

.

5

"

1

0

.

1

'(

3

0

8

c

m

)

2

0

"

8"

SIDE VIEW FRONT VIEW

Figure 3.7: Schematic of Rack Used to Hold the Source Vessel. 

45
Figure 3.8: Picture of Test Facility: Source Vessel, Discharge Piping, Collection Vessel, Support Structure 

46
0

500

10001500200025003000350040004500

-500 0 500 1000 1500 2000 25

Time [msec]

Pressure [kPa]

Source Vessel

Pipe Position #1

Pipe Position #2

Pipe Position #3

Collection Vessel

Figure 3.9: Pressure Histogram for Experimental Run #5. 

47
0

500

10001500200025003000350040004500

-500 0 500 1000 1500 2000 25

Time [msec]

Pressure [kPa]

RRR

Figure 3.10: Plot of Source Vessel Pressure Showing the Repeatability of the Experimental 

48
Capped Tees

Discharge

Valve

Union

2 - 45o Elbow

90o Elbow

PT1

PT3 PT4

PT2

Figure 3.11: Schematic of Alternate Piping Configuration #1. 

49
Discharge

Valve

Tees

Union

2 - 45o Elbow

90o Elbow

2 - 45o Elbow

Union

90o Elbow

PT1

PT2PT3

Figure 3.12: Schematic of Alternate Piping Configuration #2. 

50
0

20406080

100120140160

-500 0 500 1000 1500 2000 2500 

Time [msec]

Pressure Drop [kPa]

Figure 3.13: Pressure Drop Histogram Across a 90° Elbow from Run #9.

51
+

Diaphragm

Source Vessel

Side of Valve

Fill Level

Differential Pressure Transducer

Figure 3.14: Diagram Transducer Setup for Mass Inventory Measurement. 

52
0.0

0.51.01.52.02.53.03.5

-500 0 500 1000 1500 2000 25

Time [msec]

Cumulative Mass [kg]

Figure 3.15: Mass Inventory Histogram of Run #5 

53
- +

Negative Lead

Positive Lead Insulators

Conductive Paint (Metallic Film) 

Figure 3.16: Schematic of Thermocouple Construction. 

54
0

5

101520253035

-500 0 500 1000 1500 2000 25

Time [msec]

Temperature [C]

Source Vessel Pipe Position 

Pipe Positi

Pipe Position #1

Figure 3.17: Temperature Histogram from Film Thermocouples of Run #5. 

55
0

5

101520253035

0 500 1000 1500 2000 2500 30

Time [msec]

Temperature [C]

Collection Vessel

Source Vessel

Figure 3.18: Temperature Histogram from Shielded Thermocouples of Run #5. 

56
Section B-B

B

Ground electrode (-)

Positive Probe

Insulation

B

45 mm

Glass 

BNC connector

Guard electrode (+)

 6 mm

Figure 3.19: Schematic of Capacitance Probe for Bench Setup. 

57
level 0

Ground

electrode

capacitance 

probe

20 mm

Positive

electrode

Water

Sensing

volume

Figure 3.20: Bench Setup to Test Capacitance Sensor. 

58
-0.200

0.0000.2000.4000.6000.8001.0001.200

0 10 20 30 

Water height [mm]

void fraction 

α

Sensing

All Air volume

All Water

Figure 3.21: Void Fraction from Bench Tests of Capacitance Sensor. 

59
Brass mad

to screw in

BNC Conn

guard elect

Teflonferrule

1.2"

positive electrode

probe tip

Flow

Epoxy

Ground e

Probe tip

19mm

16mm

Insulation between guard and pipe

& between probe and guard

probe tip side view

Figure 3.22: Schematic of Capacitance Probe for Test Facility. 

60
positive electrode

Ground electrode

Figure 3.23: Cross-Section of Capacitance Setup Showing Electrical Field. 

61
0.0

0.20.40.60.81.01.2

-500 0 500 1000 1500 2000 25

Time [msec]

Void Fraction

Figure 3.24: Void Fraction Histogram of Run #5. 

62
4. CODE DEVELOPMENT WORK 

4.1 INTRODUCTION 

4.1 .1 SCOPE OF PROGRAM 

This computer program (or computer code), named FSP, was developed for the 

analysis of fire-suppressant fluid delivery systems. It was developed because there was a 

need for a flexible computer program in the public domain that would allow modeling of

delivery systems for a variety of suppressant fluids. It was developed as a part of the Next 

Generation Fire Suppression Program (NGP), managed by the National Institute of Standards 

and Technology (NIST). The primary objective of the NGP was to develop fire-suppressant 

fluids to replace Halon-1301. Hence the main focus of this computer program was to be able 

to model typical delivery systems with several fluids considered as viable alternatives to

Halon-1301. 

Typical delivery systems in current use consist of a suppressant storage vessel 

connected via a piping network to one or more delivery locations. The storage vessel is under 

pressure of an inert cover gas (usually nitrogen). The suppressant fluid is maintained in a 

liquid state by the pressurization of the storage vessel by the cover gas, usually to several 

MPa. The piping network is initially empty and exposed to the environment. The system is 

activated by a quick opening valve or a rupture disc located between the storage vessel and 

the piping network, in response to an automatic or manual trigger. This causes the gas driver

to expel the suppressant fluid through the piping network. The typical delivery time of the

suppressant fluid to the delivery location is in the range of 0.1 to 10 seconds. 

Since most suppressant fluids are in a vapor state at standard pressure and 

temperature, as the suppressant is released from the storage vessel, it will flash inside the

piping network. At the discharge location, the fluid will be a mixture of superheated vapor 

and liquid. Thus significant thermodynamic and thermal non-equilibrium can be expected 

between the phases. Because of the large pressure difference between the upstream source 

vessel and downstream exit, continuous flashing can be anticipated as the fluid travels down 

the piping system, and two-phase critical flow (choked flow) can be expected inside the pipe. 

Depending on the mass flux in the piping system and its geometrical characteristics, the 

vapor and liquid phases can be traveling at different velocities, under non-homogeneous 

conditions. Further, since the storage vessel could be maintained under pressure of the cover 

gas for significant periods of time prior to discharge of the suppressant, the suppressant 

liquid will be saturated with the cover gas. As the liquid travels down the pipe and 

depressurizes, this dissolved cover gas will come out of solution and expand. The effect of 

this on the suppressant flow would need to be accounted for. 

Hence, to accurately predict performance of the delivery system, the computer code

should be able to predict transient two-phase flows and allow for non-homogeneous, nonequilibrium conditions. The code should be able to predict two-phase choking. It should 

account for the release of dissolved gas from the liquid phase. The code should allow users to 

select one of several fluids as the suppressant of choice. Also, it should allow the ability to 

model various delivery systems, with the user being able to select the geometrical 

63 
characteristics. These would typically include the size of the storage vessel, the layout of the 

piping network, and choice of pipe sizes. And finally, the code should be non-proprietary, 

available in the public domain. These were the objectives of the computer program described 

here. 

4.1.2 BACKGROUND 

There are several available methods for calculating discharge rates of fire-suppressant

fluid delivery systems. In a summary of the state-of-the-art, DiNenno and Budnick (1988)

indicate that most flow calculations are performed by vendors using proprietary methods.

These methods are verified against standard NFPA procedures for calculations and by 

discharge testing in scaled tests. It is noted that there are significant uncertainties in the 

calculations and in extrapolating test results to actual configurations. 

Pitts et al (1994) describe a mathematical model for calculating the discharge of firesuppressant fluids from a pressurized discharge tank across an orifice at the exit and showed 

that the transient pressure history in the discharge vessel calculated by the model compared 

reasonably well to data from experiments. Elliott et al (1984) developed a computer program

called HFLOW for predicting the discharge of Halon-1301 from a discharge vessel through a 

piping system. The results of the program compared favorably to results of experiments 

conducted with Halon-1301 in several piping arrangements. The theoretical model was based 

on a homogeneous model of two-phase flow. Elliott et al (1984) also modified an existing

computer program called SOLA-LOOP (Hirt et al, 1979), which was originally written for 

steam-water flows, to be able to model Halon-1301. This computer program was a nonequilibrium, non-homogeneous model of two-phase flow, using the drift flux approach to 

calculate slip between the phases. It was shown that this program, modified to account for 

Halon-1301 properties, was able to capture the physics of the agent discharge more

accurately. 

The methods described by Pitts et al (1994) and Elliott et al (1984) appear to be the

only methods available in the public domain for calculating performance of fire-suppressant 

fluid delivery systems. One way to develop a general, flexible computer program that would 

be able to model a variety of piping networks and handle several candidate fluids would be to 

extend this work by modifying these methods. Another approach is to take advantage of the 

advances in two-phase flow modeling that have taken place in the nuclear power industry, 

which have resulted in several computer codes available in the public domain. The 

development of the computer code described here followed the latter approach. 

Based on a survey of computer codes available in the nuclear power industry (Tuzla 

et al, 1998), the RELAP5/MOD3.2 computer code (RELAP5, 1995) was chosen for further 

development for modeling fire suppressant fluid systems. This code contains most of the

basic models required for modeling fire suppressant fluid systems. It is based on a two-fluid 

model of two-phase flow and can model non-equilibrium, non-homogeneous conditions. It 

contains a model for transport of non-condensable gas through the system. It was written 

originally for water-steam systems, but the code architecture allows for adding alternate 

fluids. The code has been tested extensively against experiments for its ability to model lossof-coolant accidents in nuclear power systems and is suitable for modeling rapid two-phase 

64
transient flows. The thermal-hydraulic phenomena in the early phase of these postulated 

accidents are similar to the phenomena expected during fire suppressant fluid discharge from 

a pressurized tank through a piping network. Further, the code input is modular, allowing the 

user to choose among several available modules to optimally describe the initial and

boundary conditions required for the calculation, and the geometrical characteristics of the 

system. 

This manual provides a summary description of the FSP computer code, which was

derived by modifying the RELAP5/MOD3.2 computer code, which is referred to in this

manual as the base code or, sometimes, as RELAP5. This manual does not contain a detailed 

description of the code, particularly the numerical solution method. This is because the base 

code has been documented extensively (RELAP5, 1995) and most of the documentation 

applies to the FSP code. Hence users interested in the details of the solution method should 

refer to that documentation. 

The main modifications to the base code in developing the FSP code were: 

1) Fluid Properties: The base code was based on water-steam systems. This was 

supplemented in the FSP code by using the REFPROP fluid properties package 

(Gallagher et al, 1993). The REFPROP package, developed by the National Institute of 

Standards and Technology (NIST), covers many refrigerant fluids including those that are 

potential candidates for fire suppressant systems. In this version of FSP, four fluids have 

been incorporated, Halon-1301, CO2, HFC-227ea, and HFC-125. Any of these alternate 

fluids can be selected by the user. Water properties, which were the default option in the 

base code, have not been removed, and hence water is also one of the fluid options in 

FSP. 

2) Constitutive Relationships: The code requires specification of constitutive relationships,

such as for the calculation of wall drag and interphase drag, to provide closure for the 

basic two-phase flow conservation equations used by the two-fluid model. In the base 

code, several of these relationships were based on experimental data from water-based 

systems and were uniquely applicable to these systems. These relationships were 

modified in FSP to be fluid-independent. The specific areas modified were the interphase 

drag models in vertical flow, and wall heat transfer models for critical heat flux and 

transition boiling. 

3) Release of Dissolved Gas: The base code contained transport equations for noncondensable gas in the gas phase and also for a solute (boron) in the liquid phase. In FSP,

these transport equations were modified. The noncondensable gas transport equation was

modified to include a source term to account for release of dissolved gas. The boron 

transport equation was altered to reflect transport of dissolved gas as the solute instead of 

boron. A sink term was added to this equation for conserving the total amount of 

noncondensable gas in the system. This approach allowed modeling the release of gas

anywhere in the system based on local conditions. New constitutive relations were added 

to represent the gas release rate. This requires specification of two additional parameters, 

a critical bubble radius for initiating the gas release, and a rate constant. These parameters 

have been implemented in FSP as user input variables, and are further described in later 

sections of this manual. 

65
The above modifications were the major modifications made to the base code in 

developing the FSP code. There were some additional minor modifications, primarily for 

making the code usable on PC-based computer platforms. This Version of FSP runs on a 

WINDOWS 95/98 operating system. Installation instructions are shown in a section of this

manual. 

4.2 CODE THEORY 

This computer program, FSP, was developed to predict the transient hydrodynamic

behavior of a fire suppressant fluid as it discharges from a storage vessel through a piping 

network. It was developed from a base code called RELAP5/MOD3.2 (RELAP5, 1995), 

which was developed primarily for nuclear plant transient and accident analyses, and

therefore uses water as the working fluid. It is based on a one-dimensional, two-fluid model 

of two-phase flow and allows prediction of non-homogeneous, non-equilibrium conditions. It 

is written in Standard FORTRAN-77, and was originally developed to run on the UNIX 

Operating System, but has been adapted to other operating environments. The version of 

RELAP5 chosen for further development is generally considered to be machine-independent 

and portable to a variety of commonly used operating systems. 

To run this program, a user generates an input file. This input file is essentially the 

description of the problem to be executed. It contains the geometrical description of the 

system under analysis, the necessary boundary and initial conditions, and other constraints 

imposed by the user, such as partition of the system into control volumes (nodalization), 

selection of particular code options and control of the time step used in the numerical 

solution. The code documentation contains a user manual in Appendix A describing how to 

generate an input file. 

In executing the program, FSP generates two files, an output file that can be read by 

the user, and a restart file written in binary format. The output file provides information at 

time intervals requested by the user in the input file. The restart file contains all the

information the code needs such that the user can restart the code to run for a longer time 

period, with or without additional changes to the problem description, and is also used for 

generating graphical output. 

The RELAP5 code also contains many user-convenient features and special models, 

which were developed primarily for nuclear plant analyses, such as the definition of specific

components in nuclear systems, and code models to capture the transient behavior of nuclear 

fuel rods. Some of these are useful for the intended application to fire suppressant fluids and 

fluid systems, and some are not. In the development of the FSP code, most of the existing 

features of RELAP5 were retained, although not necessarily useful in modeling fire 

suppressant delivery systems. 

The description of the code presented here is a brief summary, intended to provide 

basic information for the user. The base code has been documented in great detail (RELAP5,

1995) and the user is referred to that documentation for additional information, particularly 

the numerical solution method. 

66
4.2.1 GOVERNING EQUATIONS 

The thermal-hydraulic model solves eight field equations for eight primary variables 

as a function of time (t) and distance (x). The eight primary dependent variables are pressure 

(P), phasic specific internal energies (Ug,Uf), vapor volume fraction or void fraction (αg), 

phasic velocities (vg,vf), noncondensable quality (Xn), and solute density (ρb). Note that the 

noncondensable quality is defined as the mass fraction of noncondensables in the gas phase,

the remaining fraction being the mass fraction of steam in the gas phase. Of these seven

variables, six are state variables (P,Ug,Uf,αg,Xn,ρb) and two are velocities (vg,vf). The

corresponding eight field equations are the conservation equations for mass, momentum and 

energy for each phase, and the conservation of mass for the noncondensables in the gas and 

liquid phases. The equations are described in detail in RELAP5-1, 1995. The simple forms of 

these equations are shown below, as one-dimensional area-averaged equations, on a per unit 

volume basis. 

Conservation of Mass: 

Gas Phase: 

g g g g g A g v

t A x = Γ∂

∂

+

∂

∂ ( ) 1 (α ρ ) α ρ + ΓN (4.2.1-1)

Liquid Phase: 

f f f f f A f v

t A x = Γ∂

∂

+

∂

∂ ( ) 1 (α ρ ) α ρ (4.2.1-2)

The terms on the right hand side of equation 4.2.1-1 are the interfacial vapor generation per 

unit volume, and the noncondensable gas released from solution per unit volume. The terms 

on the right hand side of equation 4.2.1-2 are the corresponding terms for the liquid phase. 

For continuity across the interface, Γg + Γf = 0. The source term for noncondensable gas, ΓN, 

which is the rate of dissolved gas release per unit volume, is treated similar to an external 

source term. Thus, there is no corresponding term in the liquid phase mass conservation 

equation. Essentially, this assumes that the liquid density is not affected by the presence of 

dissolved gas. 

Conservation of Momentum: 

Gas Phase: 

67




 





∂

∂ − ∂∂

+

∂

∂ − −

+ Γ − − Γ − −

+ − ∂

∂ = − ∂

∂

+

∂

∂

x

v

v

x

v

v

t

v v

C A

A v v Av A FIG v v

B A A FWG v

x

P A

x

v

A

t

v

A

f

g

g

f

g f

g f m

g gI g N g g g g f

g g g x g g g

g

g g

g

g g

( )

( ) ( ) ( )

( ) ( ) 2

1 2

α α ρ

α ρ

α ρ α ρ α α ρ α ρ

 (4.2.1-3) 

Liquid Phase: 





 





∂

∂ − ∂

∂

+

∂

∂ − −

+ Γ − − −

+ − ∂

∂ = − ∂

∂

+

∂

∂

x

v

v

x

v

v

t

v v

C A

A v v A FIF v v

B A A FWF v

x

P A

x

v

A

t

v

A

g

f

f

g

f g

f g m

f fI f f f f g

f f f x f f f

f

f f

f

f f

( )

( ) ( ) ( )

( ) ( ) 2

1 2

α α ρ

α ρ

α ρ α ρ α α ρ α ρ

 (4.2.1-4) 

On the right hand side of equation 4.2.1-3, the second term is the body force and the 

subsequent terms are the wall friction, momentum transfer due to mass transfer (from vapor

generation as well as dissolved gas release), interfacial frictional drag and force due to virtual

mass. 

Conservation of Thermal Energy: 

Gas Phase: 

wg ig g g N N gf g

g g

g

g g g g g g g

Q Q h h Q DISS

v A

A x

P

t

U v A P

A x

U

t

+ + Γ + Γ − +

∂

∂ − ∂

∂ = − ∂

∂

+

∂

∂ ( ) ( ) 1 ( ) α

α

α ρ α ρ (4.2.1-5) 

Liquid Phase: 

wf if f f gf f

f f

f

f f f f f f f

Q Q h Q DISS

v A

A x

P

t

U v A P

A x

U

t

+ + Γ + +

∂

∂ − ∂

∂ = − ∂

∂

+

∂

∂ ( ) ( ) 1 ( ) α

α

α ρ α ρ (4.2.1-6) 

68
On the right hand side of equation 4.2.1-5, the third to last terms are the wall heat transfer, 

interfacial heat transfer, energy transfer due to mass transfer (vapor generation as well as

release of dissolved gas), sensible heat due to presence of noncondensables, and dissipation. 

Conservation of Mass for Noncondensables in Gas Phase: 

g g n g N

g g n X v A

t A x

X = Γ

∂

∂

+

∂

∂

( ) 1

α ρ

α ρ (4.2.1-7)

where 

Xn= mass fraction of noncondensible gas in the gaseous phase 

ΓN= noncondensible gas generation rate 

This equation essentially conserves mass of noncondensables in the gas phase, and allows 

transport/distribution of the noncondensables along with the gas phase. It is assumed that the 

noncondensables are transported at the same velocity as the steam in the gas phase, that they 

are in thermal equilibrium with the vapor, and that the properties of the gas phase are mixture 

properties of the vapor/noncondensable mixture. The sensible heat transfer terms due to the 

presence of noncondensables in the energy equations represent heat transfer at the 

noncondensable gas-liquid interface. This is necessary because the interfacial terms use 

saturation temperature based on the local vapor partial pressure. 

Conservation of Mass for Noncondensables in liquid phase (dissolved gas): 

s f s

s v A

t A x = Γ∂

∂

+

∂

∂ ( ) 1

ρ ρ (4.2.1-8)

where 

ρs = density of the dissolved nitrogen (solute), i.e., mass of dissolved nitrogen per unit 

volume 

Γs = negative of mass rate of nitrogen released per unit volume, or solute generation rate 

This equation essentially provides the means to model the release of dissolved gas at any 

location based on local conditions. The gas release model itself is described in a subsequent 

section, along with the assumptions used in its implementation in the code. 

In the base code, it is assumed that at the gas-liquid interface, there is no storage of

mass, momentum or energy. In FSP, the dissolved gas release term is modeled similar to an 

external source. It is assumed separately that Γg + Γf = 0 (similar to the base code), and that

ΓN + Γs = 0. Summing equations 4.2.1-3 and 4.2.1-4 provides a relationship between the 

interfacial drag terms FIG and FIF. Summing equations 4.2.1-5 and 4.2.1-6 provides a 

relationship between the interfacial mass transfer and interfacial heat transfer. These are the 

essentially the interface jump conditions needed for closure of the conservation equations. 

Note that in FSP, because of the way in which the dissolved gas release term has been 

69
treated, the conservation of mass, momentum and energy at the gas-liquid interface is only

approximately satisfied. This is an inherent assumption in the code. The results of executing

the code for typical FSP applications have not shown significant effect due to this 

approximation. If new applications show significant effect, this assumption would have to be

re-visited. 

To complete the equation set, several terms on the right hand side of the conservation 

equations need to be specified, such as FWG, FIG, Qwg and Qig. These are the constitutive 

relations, described in a subsequent section. Typically, these are based on mechanistic 

models and empirical correlations, and are written in terms of velocities, thermodynamic and 

transport properties, and the temperature of the liquid and gas phases. Thus, additional 

relationships have to be defined to relate these to the primary dependent variables of the 

code. 

4.2.2 NUMERICAL SOLUTION 

To solve the above set of equations, the differential conservation equations are cast 

into a numerically convenient set of finite difference equations, using a staggered spatial 

mesh arrangement. In this method, mass and energy are conserved over a control volume

while momentum is conserved over a cell bounded by the mass and energy cell centers. A

semi-implicit method is followed for time advancement where the implicit terms are

formulated to be linear in the dependent variables at the new time. Essentially, this requires

specification of the partial derivatives of fluid density with respect to the primary dependent 

variables. Most of the source terms on the right hand side of the equations are also treated 

semi-implicitly. The noncondensable dissolved gas release term is treated explicitly. 

Figure 4.2.1 shows the staggered mesh concept used in the numerical solution. In the

method, scalar variables (pressure, internal energies and void fraction) are defined at the cell

centers, at K, L, and so on. Vector quantities (velocities) are defined at the cell boundaries, at 

the junctions connecting the cells, at j, j+1, and so on. In this method, to define the flux terms

at the boundaries of the mass and energy conservation cells, the void fraction and fluid 

properties are ‘denoted’ from the upstream cells. Similarly, to estimate friction losses needed

in the momentum conservation, ‘volume velocities’ are estimated at cell centers using an 

averaging procedure of the velocities at the junctions connecting to each cell. The 

documentation of the base code (RELAP5, 1995) provides a detailed description of the 

method. Here, it is noted that the user can access ‘junction properties’ as well as ‘volume

properties’ in developing the input model and examining the output of the code. This is 

further described in the FSP user manual shown in Appendix A. 

The governing equations, together with the interface jump conditions, constitutive 

relations, thermodynamic and transport properties, and property derivatives, as a function of 

the primary dependent variables, forms a closed set of equation. The code solves this set for 

the primary dependent variables as a function of time and spatial location. Details of

derivation of the finite difference equations and the numerical solution method are described 

in the base code documentation (RELAP5, 1995). 

70
4.2.3 STATE RELATIONSHIPS AND FLUID PROPERTIES 

As discussed previously, the code uses five primary dependent state variables. These 

five state variables are 

P: total pressure

α: gas phase volume fraction 

Ug: gas phase specific internal energy 

Uf: liquid phase specific internal energy 

Xn: mass fraction of noncondensable gas in the gas phase 

All thermodynamic variables are expressed in terms of the above five state variables. 

In addition, several state derivatives are needed for the numerical scheme. The next section

presents the detailed requirements of the model fluid property package. 

4.2.3.1 Basic Thermodynamic Property Tables

Consider first the case without noncondensables. This is referred to as the SingleComponent, Two-Phase Case. 

A. Single-Component, Two-Phase Case 

Three sets of properties are required for this case.

Liquid Phase

The following are needed in the basic thermodynamic property package: 

Temperature: Tf = Tf (P, Uf) 

Isobaric Specific Heat cpf = cpf (P, Uf) 

Isobaric Thermal Expansion Coefficient: βf = βf (P, Uf) 

Isothermal Compressibility: κf= κf (P, Uf) 

Specific Volume υf = υf (P, Uf) 

Specific Entropy sf = sf (P, Uf) 

The definitions for isobaric specific heat, isobaric thermal expansion coefficient and

isothermal compressibility can be found from standard thermodynamic text books (e.g., Van 

Wylen and Sonntag, 1985). From these basic properties, the following derivatives are

evaluated in the proposed code: density with respect to internal energy, temperature with 

respect to internal energy, density with respect to pressure, and temperature with respect to 

pressure. 

Vapor Phase

Similarly, the following are needed in the basic thermodynamic property package: 

71
Temperature: Tg = Tg (P, Ug) 

Isobaric Specific Heat cpg = cpg (P, Ug) 

Isobaric Thermal Expansion Coefficient: βg = βg (P, Ug) 

Isothermal Compressibility: κg= κg (P, Ug) 

Specific Volume υg = υg (P, Ug) 

Specific Entropy sg = sg (P, Ug) 

and the same derivatives for the liquid phase are evaluated for the vapor phase. 

Saturation Properties

In addition to the phasic properties and derivatives, the saturation properties are needed for

interface evaluation and extrapolation to metastable states (i.e., superheated liquid or 

subcooled vapor). Two tables are needed, one based on temperature and the other based on 

pressure, as shown below for the liquid phase, 

a. Give saturation temperature, Ts

: 

Saturation pressure P = P(Ts

) 

Isobaric Specific Heat cs

pf = cspf (Ts

) 

Isobaric Thermal Expansion Coefficient: βf = βs

f (Ts

) 

Isothermal Compressibility: κf= κs

f (Ts

) 

Specific Volume υf = υs

f (Ts

) 

Specific Entropy sf = ss

f (Ts

) 

where the superscript s denotes saturation state. 

b. Give saturation temperature, P: 

Saturation Temperature Ts

 = Ts (P) 

Isobaric Specific Heat cs

pf = cspf (P) 

Isobaric Thermal Expansion Coefficient: βf = βs

f (P) 

Isothermal Compressibility: κf= κs

f (P) 

Specific Volume υf = υs

f (P) 

Specific Entropy sf = ss

f (P) 

In addition, the derivative of the saturation temperature with respect to pressure is required. 

This is evaluated using the Clapeyron equation (Van Wylen and Sonntag, 1985). Similar 

saturation properties are needed for the vapor phase. 

B. Two-Component, Two-Phase Case 

In the presence of noncondensable gases, the gaseous phase (mixture of vapor and 

noncondensable gases) is assumed to be a Gibbs-Dalton mixture of vapor (either treated as an 

ideal gas or from property tables) and ideal gases. 

72
An additional variable is calculated in the code for this case. It is the noncondensable 

quality, Xn, and is defined as the mass fraction of the noncondensable gases in the gaseous 

mixture, i.e., 

Xn = Mass of Noncondensable Gases / Mass of Gaseous Phase = Mn / (Mn +Mv ) 

where Mv is the mass of the vapor phase. The gaseous phase pressure is the sum of the

partial pressure of each gaseous component, as usual. The equation of state for each 

noncondensable component is simply the perfect gas equation of state. 

The liquid properties and derivatives are calculated in the same manner as described in Part 

A above. For the gaseous phase, the independent state properties are taken as P, Ug and Xn. 

All properties of the gaseous phase are expressed in terms of these three variables. Note that 

Ug is the specific internal energy of the gaseous phase, not just the vapor phase. The

derivatives required are vapor temperature with respect to pressure, internal energy, and

noncondensable quality, as well as interface temperature with respect to pressure, internal 

energy and noncondensable quality. The interface temperature, TI, is based on the partial

pressure of the vapor, 

TI = TI (Ps) 

The evaluation procedure for the required derivatives is somewhat involved. A detailed 

discussion is given in the base code documentation (RELAP5, 1995) and not described 

further here.

4.2.3.2 Generation of Thermodynamic Property Tables

As presented in Section 4.2.3.1 above, data are required for single phase properties as 

well as saturation properties in the solution method of the code. These data need to be in 

tabular form and need to be generated for the fluid specified by the user. Once generated, 

these property tables (in binary format) are used to compute properties at any given pressure 

and internal energy through iteration, interpolation or both. The code contains the necessary

coding to perform the required interpolation calculations for a variety of options needed by 

the user. 

4.2.3.3 Implementation Of Thermodynamic Property Tables 

Various fluids have previously been implemented into the base code, RELAP5.

Examples are blood, ammonia, liquid sodium and lithium. As a result, RELAP5 was 

structured to facilitate the implementation of a new fluid. The procedure for incorporating a 

new fluid is summarized as follows:

1. RELAP5 contains a property generation program. A set of pressure and temperature 

data points that will be used for the property tables can be supplied to this program 

following a prescribed format. 

73
2. The required tables are generated via the program, which contains the necessary 

formulae, empirical correlations, canonical functions, etc. The output from this 

program contains: 

a. saturation properties as a function of pressure 

b. saturation properties as a function of temperature 

c. single phase properties as a function of pressure and temperature 

For light water substance, for instance, this particular program calls other subroutines

that calculates the Gibbs and Helmholtz functions as given in the 1967 ASME steam

tables (Meyer et al, 1997). For a new fluid, all information can be included in the 

generating program and hence avoiding the need for additional subroutines. 

3. The output from the generating program are the required property tables. The format

of the output file is as follows: 

a. Echo back the user input temperature and pressure points. 

b. Tabulate the statistics of the table including number of temperature points, 

number of pressure points, number of saturation temperature points, number 

of saturation pressure points, and number of words in the binary file. 

c. Saturation properties versus temperature. 

d. Saturation properties versus pressure. 

e. Thermodynamic properties versus pressure and temperature. For each 

pressure in the user supplied pressure range, the properties for each 

temperature.

The above information is written to a binary file. Once this binary file is available and 

linked to the code executable, the code is able to calculate properties as called for by

the various interpolation subroutines. Following this procedure, for the new fluids 

implemented into the FSP code, a generation program was written which would

directly access the REFPROP code (Gallagher et al, 1993) which contains methods

for calculating thermodynamic properties of the fluids. The generation program’s 

output binary file is then linked to the FSP code. The built-in interpolation 

subroutines in FSP then take care of generating the required properties by using the 

information in this binary file. This procedure is transparent to the user of FSP and is 

being described here for information for users who may want to implement other

fluids into the code in the future. 

4.2.3.4 Sample Binary Data File Generation for Halon 1301 

As an example, this section presents a sample generation program for Halon 1301. 

The property relations are accessed from the REFPROP V. 4.0 software package (Gallagher 

et al, 1993). REFPROP is a calculation package for a wide variety of refrigerants including 

Halon 1301 and the other fluids included in the FSP code. Halon 1301 was the test case for 

checking out the implementation procedure because there was significant experimental data 

readily available at the time the development of FSP was initiated. 

74
The major undertaking was to interface the appropriate subroutines in REFPROP 

with the coding in the generation program. FSP expects three property tables: (i) saturation 

properties given temperature, (ii) saturation properties given pressure, and (iii) single phase 

properties versus temperature at a given pressure value (i.e., isobar). The values of pressure 

and temperature are supplied by the user and are input to the generation program. The 

needed properties include specific volume (m3

/kg), internal energy (J/kg), thermal expansion 

coefficient (1/K), isothermal compressibility coefficient (1/Pa), constant pressure specific 

heat (J/kg-K) and entropy (J/kg-K). 

Figure 4.2.2 shows the flow chart of the generation program. The major steps of the 

program are described below: 

1. The critical states (pressure, temperature and specific volume) as well as the 

minimum and maximum allowable values were supplied via data statements. 

2. The user supplied temperature and pressure values are read. The input data are

examined to ensure that they are within the minimum and maximum limits specified 

in Step 1 and are in ascending order. 

3. Various pointers and table statistics are calculated. 

4. Call Subroutine table3 to generate the saturation table as a function of temperatures. 

5. Call Subroutine table4 to generate the saturation table as a function of pressures. 

6. Call Subroutine table5 to generate the single phase tables as a function of temperature 

for a given pressure. 

7. Write output to an ASCII file and a binary file. 

The interface with REFPROP are accomplished with the following steps: 

1. In the main program, initialize the REFPROP subroutines. 

2. In Subroutine table3: 

a. Calculates the saturation pressure, specific volume, enthalpy, constant 

pressure specific heat, and specific entropy of each phase, by calling

appropriate subroutines in the REFPROP package. The internal energy is then 

calculated from (h - pv), where h is the enthalpy, p is the pressure and v is the 

specific volume. Finally, the thermal expansion and isothermal 

compressibility coefficients are estimated by numerical differentiation. Thus,

for the thermal expansion coefficient, the following are calculated 

βf = (1/vf ) [vf - vf (T-dT , p)]/dT along an isobar, and 

βg = (1/vg ) [vg (T+dT , p) - vg]/dT along an isobar. 

For the isothermal compressibility coefficient, the following are calculated

κf = -(1/vf ) [vf (T, p+dp) - vf]/dp along an isotherm, and 

κg = -(1/vg ) [vg - vg (T, p-dp)]/dp along an isotherm

75
Currently, dT and dp are set at 1 o

K and 1 kPa respectively. The single phase 

specific volumes: vf (T-dT , p), vg (T+dT , p), vf (T, p+dp) and vg (T, p-dp) are 

also calculated by calling appropriate subroutines in the REFPROP package. 

3. In Subroutine table4: 

The process described in Subroutine table3 above is repeated here, except that the 

saturation temperature is calculated given the pressure value. 

4. In Subroutine table5: 

In this subroutine, the fluid state is first determined by comparing the given 

temperature to the saturation temperature to determine the fluid state. The 

corresponding single phase properties (the specific volume, enthalpy and hence 

internal energy, constant pressure specific heat, and entropy) are then calculated by

calling appropriate subroutines in the REFPROP package. Finally, the thermal 

expansion and isothermal compressibility coefficients are estimated using the 

procedure as discussed in Item 2 above. 

The exact calling sequence of each of the REFPROP subroutines is described in the Users

Manual of REFPROP. 

In addition to generating the binary property data file, several modifications are needed to 

enable RELAP5 to execute on a fluid different than water. These modifications include: 

1. Increasing the number of fluid types 

2. Adding the name of the fluid (e.g., “halon”) to the fluid type information 

3. Modifying the thermal and transport property subroutines 

With the above modifications, and the data file, the code could be used to simulate systems 

containing Halon 1301. To validate the modifications and to ensure the data file is interfaced 

properly with FSP, the following validation steps were taken: 

1. Verify that the code modifications perform as designed and do not affect the code 

inadvertently. This was accomplished by running a test case with Halon 1301 as the 

fluid type and comparing the code calculated results with the property data file. To 

verify that the new coding did not affect other parts of the code, a test case with water

as the fluid type was simulated. The results are compared to those of the same case 

simulated with the unmodified base code version. The two sets of results were 

identical as they should be. 

2. Verify the correctness of the interpolation subroutines for fluid types other than 

water. This is partially accomplished in Part 1 above. A comprehensive assessment 

of these interpolation subroutines was performed by developing a driver program to 

execute each of the interpolation subroutines interactively, and comparing the results 

to the results obtained by running the FSP code. 

76
3. Verify the correctness of the property data file. This was accomplished by comparing 

the results in the ASCII output data file with the results directly from REFPROP. 

This validation was conducted using Halon 1301 as the test case. The results of the 

validation have confirmed that the modifications and data file were correctly implemented, 

and that the various interpolation subroutines in the code performed as designed. This 

procedure was repeated for the other fluids. The property generation program is included in 

the code transmittal package for the FSP code. 

4.2.3.5 Critical Region Properties 

The REFPROP package imposes certain limits on pressure and temperature in the 

calculation of saturation properties. In particular, the temperature cannot be more than 0.97 

times the critical temperature, and pressure cannot be more than 0.9 times the critical 

pressure. These limits are set in subroutines BUBLT and BUBLP respectively. Attempt to 

remove these limits in REFPROP can produce inconsistent trends in the vicinity of the

critical region. For instance, if subroutines BUBLT and BUBLP are executed at the critical 

point, the resulting liquid and vapor values do not equal to each other, as they should be at 

the critical point. Attempt to interpolate from these inconsistent property tables can lead to 

hard failure (e.g., IEEE overflow or underflow) in the subsequent FSP calculation using these 

tables. Hence special consideration must be given in the critical region and at the critical 

point. 

The following observations can be made, regarding the critical region, which have 

been utilized in getting around the problems encountered with the use of REFPROP: 

1. At the critical point, both β and κ are infinity (inflexion point). For computation, 

some arbitrary values are assigned. 

2. At the critical point, both the vapor and liquid specific volume are equal to the critical 

volume, vcrit, which is known. 

3. For all other properties (u, s, cp), the critical values are estimated by visually

extrapolating the liquid and vapor values below the critical region. 

4. The input values for pressure and temperature satisfy the REFPROP internal limits of 

0.97 Tcrit and 0.9 Pcrit as discussed above. 

Using the above as guidelines, the following critical properties are developed for the 

candidate fluids of interest: halon-1301, carbon dioxide, HFC-227ea and HFC-125. 

Halon 1301 

vcrit = 1.34309 e-3 m3

βcrit = 0.03 K-1

κcrit = 6.0 e-8 Pa-1

ucrit = 2.6e5 J/kg 

cp,crit = 1200 J/kg-K 

77
scrit = 1.2 J/kg-K 

Carbon Dioxide 

vcrit = 2.17451 e-3 m3

βcrit = 0.023 K-1

κcrit = 6.0 e-8 Pa-1

ucrit = 3.0e5 J/kg 

cp,crit = 5000.0 J/kg-K 

scrit = 1.5 J/kg-K 

HFC-227ea 

vcrit = 1.68783 e-3 m3

βcrit = 0.03 K-1

κcrit = 6.0 e-8 Pa-1

ucrit = 3.5e5 J/kg 

Test 146 of Elliott et al was chosen for study of how this model performs in the code. 

This experiment was a discharge test of Halon-301 from a supply bottle fitted at the bottom

with a nozzle. The test measured the bottle pressure as a function of time when the flow was 

suddenly discharged from the bottle through the exit nozzle. There was no piping network 

connected to the bottle in this test. 

cp,crit = 1400 J/kg-K 

scrit = 1.3 J/kg-K 

R125 

vcrit = 1.74843 e-3 m3

βcrit = 50.0 K-1

κcrit = 1.0e-6 Pa-1

ucrit = 3.1e5 J/kg 

cp,crit = 8000 J/kg-K 

4.2.3.6 Transport Properties

Besides thermodynamic properties, a number of transport properties are needed of the 

candidate fluid. These include: 

1. Dynamic viscosity μf or g = μf or g (T,P) 

2. Thermal Conductivity kf or g = kf or g (T,P) 

3. Surface Tension σ = σ(T,P) 

These are implemented in the code as functions.

78
These properties are available from various sources, including REFPROP. In generating the 

functions needed, REFPROP was used to generate a table of data points and the data fitted

with correlations. 

An example for HFC-125 is provided below. 

This is R125 in the REFPROP package. Its chemical name is Pentafluroethane (C2HF5). 

The following information are from REFPROP. 

Tcrit = 66.18o

C 

Pcrit = 3629 kPa 

vcrit = 1.74843e-3 m3

Mol. Wt. = 120.03 kg/kg-mol 

From REFPROP, the following correlations were developed. For liquid properties, the effect 

of pressure is minimal and both viscosity and conductivity are mainly functions of 

temperatures. 

Liquid viscosity: 

μf = 2051.0 – 23.45 T + 0.2873 T2

 – 0.003754 T3

where μf is in micropoise, and T in o

C. 

Liquid conductivity 

kf = 0.0678 – 0.0004564 T 

where kf is in W/m-K and T in o

C. 

For vapor properties, viscosity and thermal conductivity are dependent on temperature and 

pressure. At atmospheric pressure, 

Vapor viscosity: 

μg = 120.52 + 0.4144 T at atmospheric pressure 

Vapor conductivity 

kg = 0.0123 + 7.93e-5 T at atmospheric pressure 

The surface tension data are as follows: 

T (o

C) Surface Tension (mJ/m2) 

-80 18.67 

-40 12.48 

0.0 6.88 

40 2.14 

79
66.18 0.0 

This was correlated as 

σ = 6.867 – 0.123 T + 0.000293 T2

 

where σ is in mJ/m2

 and T in oC., and the second coefficient was adjusted slightly (from 

0.125 to 0.123) to give zero value for σ at the critical temperature. 

4.2.4 CONSTITUTIVE MODELS 

The constitutive relations include models for defining two-phase flow regimes and 

flow-regime-related models for interphase drag and shear, the coefficient of virtual mass, 

wall friction, wall heat transfer, and interphase heat and mass transfer. Heat transfer regimes 

are defined and used for wall heat transfer. For the virtual mass, a formula based on the void 

fraction is used. 

Note that while the FSP code contains models for wall heat transfer, these models 

have not been exercised or assessed with the various new fluids that have been implemented 

into the code. Also, wall heat transfer is not expected to play a significant role in the analysis

of typical suppressant delivery systems because of the short time period of the expected

transient. Hence, the user input manual described in Appendix A does not contain the 

information needed to model wall heat transfer. If a user encounters a special application of 

FSP where wall heat transfer may need to be included in the model, the base code 

documentation (RELAP5, 1995) contains the needed information. 

The constitutive relations needed are for the terms shown in equations 4.2.1-1 through 

4.2.1-8: Γg, Γf, ΓN, FWG, FWF, FIG, FIF, C, Qwg, Qwf, Qig, Qif, Qgf, DISSg and DISSf. In 

general, these relationships are dependent on the thermodynamic state, void fraction, fluid 

properties and velocities. As shown in Figure 4.2.1, the thermodynamic state, void fraction 

and fluid properties are defined at volume centers, while the velocities are defined at the 

junctions. Fluid properties and the void fraction are ‘donored’ to junctions and volume

velocities are estimated by momentum flux weighting of the velocities at the junctions that

connect to each volume. It is convenient to use flow regimes to characterize the two-phase 

flow. This allows specifying a more mechanistic form to the various constitutive relations,

which is desirable since the code is structured to handle various different fluids. 

Four flow regime maps are included in the base code, but only two of these are 

relevant for applications of FSP. The other flow regime maps will not be encountered by the 

user and not discussed here. The flow regime maps in FSP address flow regimes and

transitions between flow regimes for horizontal flow and vertical flow. 

80

The flow regime maps are based on the work of Taitel et al, 1980, and Mishima and 

Ishii, 1984, who have developed flow regime classifications and mechanistic models to 

describe flow regime transitions. However, some of their transition criteria are quite 

complex, and further simplification has been carried out to efficiently apply these criteria in

the code. When wall heat transfer is not modeled, the flow regimes and the transition criteria 

are dependent on hydrodynamic conditions only. When wall heat transfer is modeled, there is
a possibility that the wall may encounter the condition of critical heat flux (CHF). In postCHF condition, the wall will mostly be not wetted by the liquid. This alters the flow regime. 

To capture this effect, the code recognizes the effect of wall heat transfer explicitly, and if

post-CHF conditions are indicated, a separate flow regime map is used to characterize the 

two-phase flow patterns. Again, note that the post-CHF regimes would be invoked only if the 

user models wall heat transfer, and the local fluid and wall conditions indicate post-CHF 

conditions. As discussed above, wall heat transfer is not expected to be modeled in most

applications of the FSP code. Hence, these flow regimes are not described here. The base 

code documentation (RELAP5, 1995) contains a detailed description, if needed. 

The vertical flow regime map (for both up and down flow) is for volumes whose

elevation angle φ is such that 45 < |φ| < 90 degrees. This map is modeled as nine 

regimes--four for pre-CHF heat transfer, four for post-CHF heat transfer, and one for vertical 

stratification. As noted previously the post-CHF regimes will not be encountered if wall heat

transfer is not modeled by the user. The pre-CHF regimes are the same flow regimes as

would exist if there was no wall heat transfer. For pre-CHF heat transfer, the regimes 

modeled are the bubbly, slug, annular mist, and mist-pre-CHF regimes. 

The horizontal flow regime map is for volumes whose elevation angle φ is such that 0 

< φ < 45 degrees. This map does not distinguish between pre- and post-CHF regimes. 

Otherwise it is similar to the vertical flow regime map. The horizontal flow regime map 

consists of horizontally stratified, bubbly, slug, annular and mist-pre-CHF regimes. 

In reviewing the base code, RELAP5/MOD3.2, it was found that most of the 

constitutive relations were structured to be reasonable for the fluids included in FSP. The 

only areas where correlations specific to water were used were in the computation of 

interphase drag in vertical flow, and in wall heat transfer. Thus only these were modified. In 

the description presented below, the emphasis is on providing the user with basic information 

on what is contained in FSP. Details of the various models can be found in the base code 

documentation (RELAP5, 1995). Additional information is provided below in the specific 

areas where the base code was modified. The dissolved gas release model was developed for 

FSP, and so it is described in more detail than the other models. 

4.2.4.1 Wall Friction 

The relevant wall friction terms are FWG and FWF in the phasic momentum

equations. Note that the wall friction force terms include only wall shear effects. Losses due

to abrupt area change are calculated using mechanistic form loss models. Other losses due to 

elbows or complicated flow passage geometry are modeled using energy loss coefficients

supplied by the user in the input, as shown in Appendix A, under the cards for junction 

information. These terms are also directly used in computing the energy dissipation terms,

DISSg and DISSf in the phasic thermal energy equations. 

The method is to first calculate an overall two-phase frictional pressure drop. This is

done by a standard two-phase multiplier method (as shown in Wallis, 1972). The two-phase 

friction is then partitioned between the phases (Chisolm, 1967) to get the phasic wall friction 

components. A further adjustment is made to incorporate flow regime effects, using a wetted 

wall fraction. This provides a more physical representation for certain regimes where it is 

81
clear that one phase would be more influenced by wall friction than the other, as in annular 

flow, where the wall is fully wetted by the liquid phase. 

The model uses the hydraulic diameter and surface roughness data provided by the 

user for each volume. There is an option for the user to turn off wall friction altogether, as

shown in Appendix A, but this is not recommended unless as a sensitivity study to debug 

results from a specific analysis. 

4.2.4.2 Interphase Friction 

The interphase friction (or drag) per unit volume in the phasic momentum equations

is expressed in terms of a phasic interphase friction coefficient as: 

Fig = αgρgFIG(vg – vf) 

and 

Fif = αfρfFIF(vg – vf) 

The magnitude of the interphase drag force per unit volume on the liquid is assumed to be 

equal to that on the gas phase. This provides a relationship between FIG and FIF, 

αgρgFIG = αfρfFIF 

The interphase drag force per unit volume represents the interphase momentum transfer rate

and is due to the difference between the phasic velocities. Thus, analogous to the transfer of 

heat or mass, it can be expressed as: 

Fi,gf = Fi (vg – vf) 

Where, Fi is a momentum conductance term, and can be expressed as: 

Fi = AgfBgf

Where Agf is the interfacial area per unit volume, and Bgf is a momentum transfer coefficient. 

The code uses the flow regime characterization to derive relationships for Agf and Bgf

in terms of fluid properties, flow geometry, phasic velocities and the void fraction. These are 

then related back to the terms FIG and FIF needed in the phasic momentum equations. 

The base code, RELAP5/MOD3.2, uses mechanistic models for the above 

calculations, except for vertical flow in the pre-CHF regimes. In this case, the base code uses

a correlation developed empirically from water data. In FSP, to make this independent of the 

fluid type, the empirical correlation (called the Drift Flux correlation or the EPRI correlation) 

has been replaced with the mechanistic models. The models are similar to what exists in the 

base code for the various flow regimes and therefore are not further described here. It is

noted that for the calculation of interphase drag in vertical flow, in FSP, the vertically 

82
stratified regime is not explicitly used. Vertical stratification is recognized primarily for the 

calculation of interphase heat and mass transfer.

In building the input file, under the volume control flag information, there is a

parameter, b, that allows the user to select an optional interphase drag modeling approach for

rod bundles. In FSP, in addition to replacing the interphase drag models for vertical flow, this 

option has been essentially de-activated. Thus choosing a value of 0 or 1 for this parameter 

will not have an effect. It is recommended that the user choose the default value, which is 

zero. 

4.2.4.3 Interphase Heat and Mass Transfer 

The interphase heat and mass transfer rates are related by the assumption that energy

is not stored at the interface. It should be noted that the code separately recognizes interphase 

mass transfer in the bulk fluid and at an interface near the wall. This allows the code to 

mechanistically treat some special cases, such as subcooled boiling. In this case, there is a

radial temperature gradient within the liquid, and boiling occurs near the wall while there is 

condensation in the bulk fluid. Thus, the general treatment is: 

Γg = Γig + Γw

where Γw is computed by the wall heat transfer models. When wall heat transfer is not 

modeled by the user, this term is zero. Γig is computed from the heat transfer at the bulk 

interface. The interface is assumed to be at the saturation temperature. The heat transfer at 

the bulk interface is then: 

Qig = Hig (Ts

 – Tg) 

Qif = Hif (Ts

 – Tf) 

Qig + Qif + Γig(hg – hf) = 0 

Where hg and hf can be at saturation or at the value corresponding to the local pressure and

temperature, depending on whether the process is one of boiling or condensation. 

The terms Hig and Hif are the interphase heat conductance terms, and are conveniently 

modeled as a product of the interphase area per unit volume and the interphase heat transfer 

coefficient. These are in turn determined as a function of the flow regime, and modeled as a 

function of fluid properties, flow geometry, velocities and the void fraction. 

4.2.4.4 Wall Heat Transfer 

The total wall heat flux, q”, is modeled as a sum of heat fluxes to the vapor and liquid 

phases: 

q” = hg (Tw – Tref,g) + hf (Tw – Tref,f) 

83
where hg and hf are heat transfer coefficients to the gas and liquid phases respectively (per 

unit wall surface area), and Tref,g and Tref,f are gas and liquid reference temperatures. The 

reference temperatures can be the local gas or liquid temperature or the saturation 

temperature, depending on the correlation being used for the heat transfer coefficient.

A heat transfer surface concept is used, based on the heat flux, wall temperature and 

saturation temperature, to decide the mode of heat transfer (such as convection to single 

phase liquid, bulk nucleate boiling, condensation). Correlations are used for the heat transfer 

coefficients in the various heat transfer regimes to compute the heat flux. The total heat flux

is ensured to be continuous over the heat transfer-wall temperature surface, to avoid 

numerical difficulties as the heat transfer regimes go through transitions. 

Also calculated is the portion of the wall heat flux (Qw

ig and Qwif) that contributes

directly to mass transfer. Then the interphase mass transfer due to wall heat transfer is

calculated from:

Qw

ig + Qwif + Γw(hg – hf) = 0 

where, as before, hg and hf can be either at saturation conditions or at the local pressure and 

temperature, depending on whether the process is boiling or condensation. 

The base code, RELAP5/MOD3.2 contains generic models for the heat transfer 

coefficients in most of the heat transfer regimes, that can be considered to be reasonably 

fluid-independent. In two areas, the correlations are more empirical and also based solely on 

water data. Hence, these have been modified in the FSP code. One area is the critical heat 

flux (CHF). The base code uses a table look-up procedure for calculating CHF, using data 

obtained for water at a variety of conditions. This was replaced in FSP by a simpler model, 

based on a modification to the Zuber correleation(Zuber, 1961). 

qCHF = 1.31 (0.96 - αg) ρg hfg [σg (ρf-ρg)/ρg

2

]

1/4

Zuber’s correlation has been shown to match CHF data for a variety of fluids under pool 

boiling conditions. The modification to Zuber’s correlation is the term (0.96 - αg) and has 

shown to match flow data with water and Freon. Note that this condition of CHF is normally 

one that is associated with hydrodynamic dryout of a heated surface. This is the type of CHF 

one may expect when wall heat transfer is modeled with only passive heat structures, such as

typical pipe walls in fire-suppressant delivery systems. In implementing the above model in 

FSP, a lower limit of 1000 W/m2

 was imposed to prevent negative values, at high void 

fractions. 

The other area modified is in the correlation for wetted area fraction used for 

computing post-CHF heat transfer. This correlation utilizes an exponential decay form to 

capture the decrease of wetted area fraction as the wall temperature increases. The decay

constant was derived as a function of flow rate and void fraction, from water data. In FSP, 

the decay constant has been set to a fixed value. This is considered reasonable because this 

heat transfer regime is rarely encountered, and will provide continuity for the heat transfer

surface if encountered. 

84
Wall heat transfer is not expected to be significant for most FSP applications, because 

of the very short transient times expected. Hence, the wall heat transfer models were not 

tested for the various fluids implemented in FSP. The input manual, shown in Appendix A, 

does not describe how the user can model wall heat transfer. Nevertheless, the models do

exist in the code. To access this feature, the user is currently directed to the base code 

documentation (RELAP5, 1995). 

4.2.4.5 Coefficient of Virtual Mass 

The last term in the phasic momentum equations is the dynamic drag term due to 

virtual mass effects, also called the inertial drag force per unit volume. In the base code, in

the numerical solution procedure, the spatial dependence of the velocities in this term were 

neglected. Thus the term only includes the time-dependent portion. This was a compromise, 

because the spatial terms were dependent on the nodalization selected by the user to represent 

a particular system, and there was some user experience that suggested that this led to 

numerical errors. The main intent of this term is to capture sudden accelerating flows, and the

time-dependent term provides that capability. The coefficient, C, in this term is called the

coefficient of virtual mass. It is given as: 

C = 0.5 (1 + 2αg) / (1 - αg) for 0≤αg≤0.5 

and, 

C = 0.5 (3 - 2αg) / αg for 0.5≤αg≤1 

4.2.4.6 Release of Dissolved Gas 

The liquid agent in the storage vessel is stored under pressure of a noncondensable 

gas (nitrogen). It is assumed that the agent has been stored for a long enough period for the

liquid to be saturated with nitrogen (equilibrium condition). At equilibrium, the concentration 

of dissolved nitrogen in the liquid medium is proportional to the partial pressure of nitrogen 

in the gaseous space above the liquid. This is Henry’s law: 

xN

*

 = H PN

where xN

*

 is the equilibrium concentration (mole fraction of solute in solution), H is Henry’s

constant and PN is the partial pressure of nitrogen above the liquid. H is usually a function of

temperature and the liquid type. 

As the system is depressurized, nitrogen will come out of solution since the partial 

pressure decreases. According to Henry’s law, the corresponding concentration of nitrogen 

in the solution must also decrease. Hence the excess must be released. In a fast 

depressurization transient, the excess can still be in solution due to surface tension. 

Considering a micro-bubble of dissolved gas with radius rc, the difference between the 

pressure inside and outside the bubble is given by 

85
∆PNUC = 2σ/rc

where σ is the surface tension. As the pressure outside decreases, the pressure difference 

builds up, and until the point where surface tension can no longer hold the bubble together, 

and the excess nitrogen will then be released. This initiation mechanism is similar to

nucleate boiling, in which the liquid is superheated before boiling can occur. The radius of

the micro-bubble can be empirically determined by observing the point at which nitrogen 

release occurs. In this model, this radius is assumed known (user input). The release 

pressure is then given by 

Prel = Pi - ∆PNUC = Pi - 2σ/rc

where Pi is the initial system pressure. 

The above determines the initiation of release. The rate of release is determined as 

follows. Henry’s law determines the equilibrium concentration at given conditions. As the 

system is depressurized, the corresponding partial pressure decreases. The corresponding 

equilibrium concentration will be less. The solution is therefore at a state of non-equilibrium. 

The release is a means to achieve the new equilibrium again. Thus, the release rate is

postulated to be proportional to the difference between the actual concentration and the 

equilibrium concentration (dictated by Henry’s law) as follows: 

Γrel = Γ0 (xN – xN

*

) 

where Γ0 is some release coefficient, xN is the actual concentration, and xN

*

 is the equilibrium

concentration as given by Henry’s law. Thus 

Γrel = Γ0 (xN – H PN) 

The release will be zero when the actual concentration approaches the equilibrium 

concentration, as it should be. The release model is similar to other transfer phenomena, 

such as heat transfer and mass transfer. 

The model requires that user specify, as input, the parameters rc and Γ0, and the initial 

pressure Pi and temperature Ti. 

The numerical implementation of the release rate model poses a problem, since the partial 

pressure is not the local value, but the value above the liquid-gas interface. A compromise is 

made by using the local total pressure in place of the partial pressure in calculating the 

equilibrium concentration using Henry’s law. This over-predicts the equilibrium

concentration, and hence under-predicts the release rate. It may be compensated for by using 

a larger release coefficient. Thus, the release rate is numerically implemented as, for the kth

volume: 

Γrel, k = Γ0 (xN,k – H Pk) 

86
Note that HPk is the mole fraction of nitrogen in solution. To convert to mass of nitrogen per 

unit mass of solvent, one gets 

Mass fraction = HPk WN / (1-HPk)Wf

where WN is the molecular weight of nitrogen (28.0) and Wf is the molecular weight of the

fluid type. The actual concentration is related to the solute density as 

xN,k = (ρs/αfρf )k

The solute generation rate is therefore the negative of the release rate: 

Γs,k = - ΓN,k = -Γ0 { (ρs/αfρf )k– H Pk WN/ (1-H Pk)Wf} 

The above is only valid if the local concentration is greater than the equilibrium

concentration. Otherwise the release rate should be zero. Hence 

Γs,k = min (Γs,k, 0.0) 

Also, the amount of release in one time step should not be more than the total amount

dissolved in the volume. Thus the maximum release rate is given by 

Γrel,k, max = ρS/∆t 

where ∆t is the current time step. Hence the final expression for the solute generation term is 

Γs,k = max {min (Γs,k, 0.0), -ρS/∆t} 

The Henry’s constant is estimated from the PROFISSY Code (se Yang et. Al., 1997) as 

follows. Assuming Henry’s law to be valid, the equilibrium concentration is given by 

xN = HPN

In PROFISSY, the concentration of nitrogen in the agent liquid (xN, mole fraction) is 

calculated given the system pressure (P), temperature (T), total volume (V), and total agent 

mass (ma). The partial pressure of nitrogen can be determined from

PN = P – Pv 

where Pv is the saturation pressure of the agent fluid at T. Table 4.2.1 lists the output from 

PROFISSY, and the calculation of H as described above. Three fluid systems are

considered: HFC-227ea, HFC-125 and Halon 1301. As can be seen, the solubility of 

nitrogen in all three agents is very similar. Assuming only temperature dependence, the 

Henry’s constant can be correlated as a quadratic, as shown in Figure 4.2.3. The analytical 

form of the correlation is given by 

H = H0 + C1T + C2T2

87
where H is in mol/Pa, and T in K. 

and, 

H0 = 2.347767 x 10-7

C1 = -1.55063 x 10-9

C2 = 2.957799 x 10-12

Figure 4.4 shows the calculated bottle pressure for various release coefficients. The 

critical radius was specified at 7.5e-9 m. The initial concentration was estimated as 0.02916 

kg N2/ kg liquid phase. This value is comparable to the reported initial concentration of Test 

146. At about 0.1 seconds, as the pressure decreases to a value of about 4 MPa, there is a

noticeable change in the rate of pressure decrease. This is the pressure at which gas begins to 

be released, and experimental results generally indicate that this would happen at about the 

same pressure for a given fluid type and system configuration. As discussed earlier, this 

pressure is controlled by the value of the critical gas bubble radius, rc, chosen by the user. 

The lower the value of rc specified, the lower the pressure at which gas release would be 

initiated. 

The value of Γ0 controls the rate of gas release. At low values of Γ0, the rate of

release is small, and the pressure continues to decrease monotonically below the pressure at 

which gas release was initiated. As the value of Γ0 is increased, there is a noticeable decrease 

in the depressurization rate at the gas release pressure (at about 0.1 seconds). This is the 

expected trend. Test runs of the code indicate that this method of modeling the gas release 

works relatively smoothly, and the code appears to handle a wide range of release 

coefficients without numerical difficulty. 

One of these sensitivity case, corresponding to Γ0 of 10,000 kg/m3

-s and an initial

concentration of 0.02196 kg N2/kg Liquid, was used to further check out the model. Figure 

4.2.5 shows the calculated pressure profile as compared to data. Figure 4.2.6 shows the void 

fraction responses for the various cells used to model the supply bottle. Volume 203-01, 

shown in the figure, is at the top of the bottle, and Volume 203-10 is at the bottom. The 

figure shows a relatively smooth emptying of the bottle, with the void fractions in the lower 

cells increasing to 1.0 subsequent to the void fractions in the cells above them. The void 

fraction in the last cell, Volume 203-10, increases to 1.0 after all the other cells, at about 1.0 

second. This represents the time in the calculation when the bottle is emptied. 

This calculation was also used to check out other parameters. The solute

concentration profiles showed smooth changes, following expected trends. A mass balance 

was carried out to ensure that the code conserved the mass of nitrogen, between the initial 

amount and the additional amount released from the liquid. Essentially, these checks ensured

that the gas release model was implemented in the code as intended, and the mechanistic gas

release model produced the expected trends. There is of course some uncertainty in how to 

specify the correct values for the parameters rc and Γ0. However, as shown in other sections

of this report, comparisons to experiments with different fluids and system configurations 

indicate reasonably good agreement between code calculations and experimental data, with 

these parameters varying over a relatively narrow range. This has been useful in arriving at a 

range of recommended values for these parameters. Further, from a design perspective, while 

88
variation of these parameters causes a noticeable variation in the initial pressure history, there

is not much variation in the time required to empty the supply bottle. Hence it is believed that 

this gas release model, used with the range of recommended values for rc and Γ0, would 

adequately serve the needs of a design engineer evaluating the performance of a suppressant 

delivery system. 

4.2.5 SPECIAL PROCESS MODELS

This section describes certain models in the code that simulate special processes. The 

user can select many of these processes via input. A detailed description of these process 

models can be found in the base code documentation (RELAP5, 1995). The information 

below is a brief summary, to provide some basic information and where appropriate, some

user guidelines on the use of these models. 

4.2.5.1 Choked Flow

Choking is defined as the condition where the mass flow rate becomes independent of 

the downstream conditions. Choking occurs because acoustic signals can no longer propagate 

upstream as the fluid velocity equals or exceeds the propagation velocity. In typical FSP 

applications, choking can occur at the bottle exit, as soon as the discharge begins, and also at 

the discharge location from the system. The choked flow model in FSP is by Trapp and 

Ransom(1982) and is used to calculate the mass discharge from the system. The model can 

also be used to predict the existence of and calculate choked flow at interior points in the 

system. 

The user can select the choking model as a part of the junction control flag 

information in the input file. The user can select this option at all junctions. In fact, the

default option will select the choked flow model. It is recommended, however, that the user 

only select this model at the system exit location(s) and at the junction representing the exit 

of the supply vessel. The model should be turned off at all other internal junctions. This is 

because the sonic velocity (propagation velocity) in two-phase flow can be a strong function

of the void fraction, and the choked flow model imposes an external criterion on the code’s 

normal numerical solution method. This can lead to difficulties in code execution and a

drastic reduction in the computational time step. 

Also note that the input for junctions allows the use of discharge coefficients for

subcooled liquid, saturated fluid and single phase vapor. Current assessment of the code 

against experiments suggests that best results are obtained by using a subcooled discharge 

coefficient of 1.0, a two=phase discharge coefficient of 0.8 to 0.9, and a single phase vapor 

discharge coefficient equal to the two-phase discharge coefficient. 

4.2.5.2 Horizontal Stratification and Entrainment 

Flow at low velocity in a horizontal pipe can be stratified due to buoyancy forces. The 

horizontal flow regime map in the code contains a mechanistic model to predict when 

stratification would occur. When the flow is stratified, the area average pressures are affected 

by nonuniform transverse distribution of the phases. The normal code equations are based on 

assuming uniform transverse pressure distribution. To account for the nonuniform pressure 

89
distribution, adjustments are made to the momentum equations to consider area average 

pressures for the vapor and liquid phases, and the interfacial pressure between them. This

model is always active in the code, but becomes significant only when the code’s 

stratification criterion is met. The code then calculates a liquid level in the channel(circular

pipe geometry is assumed) and uses this information to estimate the hydrostatic pressure in 

the transverse direction. 

Also, when the flow is stratified, the void fraction flowing in a junction may be 

different than the upstream volume void fraction. A typical example would be when there is 

an area contraction. With stratified conditions in the upstream volume, if the level is below 

the pipe centerline, the void fraction in the junction would be significantly higher than the 

upstream volume. And if the level is above the pipe centerline, the junction void fraction 

would be lower than the upstream volume. Also, in addition to the formation of a level and

its effect, the vapor can entrain liquid through the junction, and similarly, the liquid can pull 

through some gas. The code considers these possibilities, and adjusts the junction void 

fraction to be different than the upstream volume void fraction when the upstream volume is

calculated to be stratified.

The user does not actively engage these models. It is part of the code’s internal

structure and will automatically be invoked as needed. 

4.2.5.3 Abrupt Area Change 

A general piping network system can contain sudden area changes and orifices. The

code contains built-in models for the user to access in simulating these geometries. The basic 

hydrodynamic model in the code is formulated for slowly varying flow area variations. 

Hence, special models are needed when the flow area undergoes an abrupt change, such as at 

an orifice. 

Basically, the user can select the abrupt area change option at any junction in the 

model, as a part of the junction control flag input information. The code will calculate the 

energy loss coefficients associated with the abrupt area change internally. In a case where the 

flow is going through a sudden contraction, the code will also internally compute the losses

associated with the formation of a vena contracta. Along with selecting the abrupt area 

change option, if the user also enters form loss coefficients, these will be added to the

internally calculated form losses. 

4.2.5.4 User-Specified Form Losses 

In the junction control flag information card, the user can specify a junction to be a 

smooth or abrupt area change junction. The user can also specify form loss coefficients at the 

junction. For standard geometries, and configurations such as elbows and tees, these form

loss coefficients can be obtained from standard handbooks. If the user has frictional pressure

drop data for a particular system at known flow rates, that can be used to back out form loss 

coefficients.

90
The user can also enter flow dependent form loss coefficients. Typically, the form 

losses would be a function of the Reynolds Number. The code input allows a standard form

and the user has to enter the desired constants for the Reynolds Number dependence. This is 

shown in the code input manual, in Appendix A. 

4.2.5.5 Cross Flow Junction 

The code’s numerical method is generally formulated using one-dimensional 

elements. However, there are applications where an approximate treatment of transverse flow 

can provide an improved physical simulation. The crossflow junction component provides 

this feature.

A typical place to use the crossflow junction would be in modeling a standard 90o

Tee. The crossflow junction would be used to model the connection to the branch line, as 

shown in Figure 4.2.7a. In this case, the momentum flux in the side branch is assumed to be 

perpendicular to the main flow stream. Thus the main stream momentum flux does not 

contribute to the crossflow momentum formulation. Essentially, in the crossflow junction 

formulation, there is no transport of axial direction momentum due to the flow in the 

transverse direction. 

The user selects the crossflow junction option by specifying the appropriate

connection code when describing the ‘From’ and ‘To’ volumes for a junction. This is 

described in the input manual in Appendix A. 

4.2.6 COMPONENT MODELS AND SYSTEMS 

The code contains a variety of generic components that are used to build system 

models. These include components such as Single Volume, Time-Dependent Volume, Single 

Junction, Pipe, Valve, etc. Most of these are self-explanatory. This section provides some

additional information on building a system model and the use of some of the components 

and features. 

4.2.6.1 Branch 

There are several ways to simulate branching flows. Several single junctions can be 

connected to the exit of a single volume to simulate flow splits. A BRANCH component is

also provided in the code. Essentially, this functions no different than a single volume to

which more than two junctions are connected. Primarily, it is a user-convenience, which 

allows all the junction information to be entered along with the volume information in one 

component instead of using several separate components. A crossflow junction can be used 

to connect a Branch component to another component. Two typical uses are in modeling 

branching flows in a one-dimensional branch, and a Tee branch. 

A one-dimensional branch is one where more than two normal junctions connect to a

volume, and where it is assumed that multidimensional effects are small compared to system

interaction effects. Typically, this would be used in branched flows that occur in headers or

91
plena. Fluid entering a header and then splitting into several parallel paths would be modeled 

as a one-dimensional branch. 

A Tee branch is formed by the arrangement shown in Figure 4.2.7a, where the 

crossflow junction was used. In this case, the Volume V2 may be specified as a Branch 

component, and the crossflow junction connecting to it could be defined by the appropriate 

connection code, described in the input manual in Appendix A. The Tee connection could 

also be modeled without a crossflow junction, by using a Branch component with normal

junctions, as shown in Figure 4.2.7b. This approach has the advantage that fewer volumes are

used. The disadvantage is that if the Tee is really a Standard 90o

 Tee, then it may not be the 

best representation of the physical process. In typical FSP applications, very little difference 

has been observed between the two approaches. Hence, either approach seems acceptable. 

Note that a volume with more than two junctions connected to it is treated by the code

as a Branch component, even if it has not been defined as such in the input. What this really 

means is that the volume velocity in this volume is determined as a momentum flux weighted 

average of the flows in the junctions connecting to it. Also, the volume cross-sectional area is 

apportioned between the various streams in proportion to the volumetric flow in these 

streams. The consequence is that, if one stream has a very high velocity, there would be 

excess error in conserving momentum (and mechanical energy) around this volume. One 

place where this has an impact is if one of the junctions connecting to the Branch is the exit 

junction for the system and it is choked at a high velocity. To avoid this situation, it is 

recommended that the user allow some straight pipe length connected to all the junctions of a

Branch. It is also recommended that when the BRANCH component is used, the junctions of

the component use the abrupt area change option. 

4.2.6.2 Valves 

Valves are quasi-steady models that are used to either specify an option in a system

model or to simulate control mechanisms. They are essentially treated as single junctions and 

the input required to describe them is similar. There are two categories of valves, those that 

open and close instantly or those that open and close gradually. The base code contains 

several types of valves. Only two types are envisioned to be useful in modeling typical fire 

suppressant delivery systems, and these are described here. 

A Trip valve is one that is fully open or fully closed based on a trip setting. The 

opening or closing will occur in one time step following the trip signal. The trip signal can be 

structured to latch the valve in an open or closed position, if desired, as described in the input 

manual in Appendix A. 

A Motor Valve has the ability to control the junction flow area as a function of time.

The operation of this valve is controlled by two trips, one to open and one to close. A userspecified rate parameter controls the rate at which the valve area changes. The abrupt area 

change option is normally used with this component to allow the code to automatically 

compute the form loss coefficient as a function of the valve area. 

92
Typical use of this component in FSP applications would be to describe the junction

that connects the suppressant supply vessel to the piping system. The valve could be opened 

as fast as desired by the user to simulate the initiation of the discharge transient. 

4.2.6.3 Trip Systems 

The Trip system consists of the evaluation of logical statements. Each trip statement 

has a true or false result, and an associated variable, TIMEOF. This variable is the time in the 

transient at which the trip became true and has a value of –1.0 when the trip is false. This 

variable can be used to set time delays based on events during the transient. Within the code

structure, the trip system only evaluates the logical statements(such as , is transient time

greater than 2.0 seconds?) The decision of what action is needed based on the trip status

resides within other models. In the above example, the trip system would assign a true or 

false value to a particular trip number depending on whether the transient time was less or 

greater than 2.0 seconds. If this trip number is invoked by a valve component, the valve 

models would decide whether to close or open the valve, at what rate, etc.

Using the trip system is fairly self-explanatory, and directions are provided in the 

input manual in Appendix A. 

4.2.6.4 Control Systems 

The control system provides the capability to evaluate simultaneous algebraic and 

differential equations. This allows simulation of control systems typically used in 

hydrodynamic systems. Another use is to define auxiliary output quantities, which can then 

be extracted from the output or restart files as tables or plots. 

In typical FSP applications, there are no complicated trip and control logic required. 

Hence, control systems are primarily used to define auxiliary output variables. For example,

enthalpy is not a normal code variable, and if it is needed for some reason, it can be 

computed using several control systems (such as the SUM and MULTIPLIER components) 

to calculate it from the internal energy, pressure and specific volume, which are all normal 

code variables. Another example would be if a user wanted to know the Froude Number at

the system exit, for evaluating mixing characteristics in the discharge area. Similarly, the 

INTEGRAL component could be used to keep track of suppressant mass delivered. A variety 

of control systems are provided in the input manual in Appendix A, with directions on how 

to enter the needed cards in the input file. 

Table 4.1: PROFISSY Results 

HFC-227ea/N2 Halon/N2 R125/N2

T(K) H (mol/Pa) H(mol/Pa) H(mol/Pa)

190 n/a 4.70357E-08 n/a

200 n/a 4.20328E-08 n/a

210 n/s 3.87112E-08 3.64965E-08

220 3.76614E-08 3.67283E-08 3.79995E-08

93
230 3.61976E-08 3.58968E-08 3.36686E-08

240 3.52156E-08 3.29835E-08 3.36637E-08

250 3.47006E-08 3.44772E-08 3.46937E-08

260 3.3033E-08 3.1454E-08 3.41243E-08

270 3.34344E-08 3.26477E-08 3.24578E-08

280 3.28759E-08 3.2787E-08 3.47492E-08

290 3.28942E-08 3.21827E-08 3.39495E-08

300 3.35448E-08 3.35049E-08 3.76173E-08

310 3.48921E-08 3.74043E-08 3.90425E-08

320 3.57496E-08 3.96115E-08 4.44921E-08

330 3.78611E-08 4.87664E-08 5.39104E-08

340 4.14481E-08 n/a n/a

350 4.59199E-08 n/a n/a

360 5.41519E-08 n/a n/a

n/a: outside PROFISSY data range for the fluid 

94
Vector Node, 

Mass and Energy

Cell 

j-1 j 

CL

Flow 

Scalar Node, 

Momentum

Cell 

K L 

Figure 4.2.1: Staggered mesh numerical method 

. 

95
Figure 4.2.2: Flow Chart of Generation Program 

EXIT Inconsistent

EXIT 

Output

Data 

Call Table5

given T, P 

Call Table4

given P 

Consistent

Call Table3

given T 

Test Input 

 Data 

Read in Pressure and 

Temperature Data

Provide Data 

Statements

BEGIN 

96
Figure 4.2.3: Curve Fit for Henry Constant 

97
Figure 4.2.4: Effect of Release Coefficient on Calculated Pressure Response 

98
Figure 4.2.5: Comparison of Calculated and Measured Bottle Pressure for Test 146 

99
Figure 4.2.6: Void Fraction Responses in Lower Portion of Bottle 

100
Figure 4.2.7a: 90-degree Tee model with crossflow junction

Volume

V2

Volume

V1

Junction J1 Junction J2

Volume

V2

(Branch) Junction J1

Junction J2

Volume

V1

Volume

V3

Volume

V3

Junction J3

Crossflow 

Volume

V4 

Figure 4.2.7b: Tee model with normal junctions 

101
4.3 CODE ARCHITECTURE 

This section provides an overview of the code architecture and briefly summarizes 

computer adaptability, code organization, input processing, and transient operation. 

4.3.1 COMPUTER ADAPTIBILITY 

The FSP code was developed from the base code RELAP5/MOD3.2, which was 

written in FORTRAN 77 for a variety of 64-bit and 32-bit computers. Here, a 64-bit 

computer is one in which floating point, integer, and logical quantities use 64-bit words; a 

32-bit machine uses 32-bit words for those same quantities but also allows 64-bit floating 

point operations. Examples of 64-bit computers are Cray, Cyber-NOS-VE, and DEC/Alpha

machines. Examples of 32-bit computers include IBM mainframes, such as a 3090, 

workstations, including those from DEC, HP, IBM, SGI, and Sun, and many personal 

computers. 

In developing FSP from the base code, RELAP5/MOD3.2, the source code was

obtained from the Idaho Engineering National Laboratory. This code was adaptable to other 

computing environments. The code was configured to work on a PC with WINDOWS 95 

operating system and tested to ensure that the code options would function properly. The 

hardware chosen for installation was a laptop PC with a Pentium 233 MHz processor, 64 MB

RAM and 3.2 GB hard disk drive. The compiler used was the DigitalTM Visual FORTRAN 

compiler (Standard Edition, Version 5.1) for WINDOWS 95 operating system. 

4.3.2 TOP LEVEL ORGANIZATION 

The coding is modular, using top-down structuring. The various models and 

procedures are isolated in separate subroutines. The top level structure is shown in Figure 

4.3.1 and consists of input (INPUT), transient/steady-state (TRNCTL), and stripping (STRIP) 

blocks. The input block (INPUT) processes input, checks input data, and prepares required 

data blocks for all program options and is discussed in more detail in Section 4.3.3. The 

transient/steady-state block (TRNCTL) handles both transient and the steady-state options. 

The steady-state option determines the steady-state conditions if a properly posed 

steady-state problem is presented. Steady state is obtained by running an accelerated transient 

until the time derivatives approach zero. Thus, the steady-state option is very similar to the 

transient option but contains convergence testing 

algorithms to determine satisfactory steady state, divergence from steady state, or cyclic 

operation. If the transient technique alone were used, approach to steady state from an initial 

condition would be identical to a physical transient from that initial condition. Pressures, 

densities, and flow distributions would adjust quickly, but thermal effects would occur more

slowly. To reduce the transient time required to reach steady-state, the steady-state option

artificially accelerates heat conduction by reducing the thermal 

capacity of the conductors. The transient/steady-state block is discussed in more detail in 

Section 3.4. Note that in running typical FSP problems, the user will not be using wall heat 

conductors and will typically choose the transient option, allowing a ‘null transient’ time for

steady-state initialization.

102
The strip block (STRIP) extracts simulation data from a restart plot file for

convenient passing of simulation results to other computer programs, for plotting or for

examining the results in some other fashion. 

4.3.3 INPUT PROCESSING 

The code provides detailed input checking for all system models using three input 

processing phases. The first phase reads all input data, checks for punctuation and typing 

errors (such as multiple decimal points and letters in numerical fields), and stores the data

keyed by card number such that the data are easily retrieved. A list of the input data is 

provided, and punctuation errors are noted. 

During the second phase, all input data are processed. Some processed input is stored 

in fixed common blocks, but the majority of the data are stored in dynamic data blocks that

are created only if needed by a problem and sized to the particular problem. In a NEW-type 

problem, which is the typical problem for users of FSP, dynamic blocks must be created.

Also during this phase, restart data from a previous simulation are read if the problem is a 

RESTART type. Most FSP problems will not use this feature. In any case, in RESTART 

problems, dynamic blocks may be created, deleted, added to, partially deleted, or modified as 

modeling features and components within models are added, deleted, or modified. Extensive 

input checking is done, but at this level, checking is limited to new data from the cards being 

processed. Relationships with other data cannot be checked because the latter may not yet be 

processed. As an illustration of this level of checking, junction data are checked to determine

if they are within the appropriate range (such as positive, nonzero, or between zero and one) 

and volume connection codes are 

checked for proper format. No attempt is made at this point to check whether or not 

referenced volumes exist in the problem until all input data are processed. 

The third phase of processing begins after all input data have been processed. Since 

all data have been placed in common or dynamic data blocks during the second phase, 

complete checking of interrelationships can proceed. Examples of cross-checking are

existence of hydrodynamic volumes referenced in junctions and heat structure boundary 

conditions; entry or existence of material property data specified in heat structures; and 

validity of variables selected for minor edits, plotting, or used in trips and control systems. 

As the cross-checking proceeds, cross-linking of the data blocks is done so that it need not be 

repeated at every time step. The initialization required to prepare the model for the start of

the transient advancement is done at this level. 

Input data editing and diagnostic messages can be generated during the second and/or

third phases. Input processing for most models generates output and diagnostic messages

during both phases. Thus, input editing for these models appears in two sections. As errors 

are detected, various recovery procedures are used so that input processing can be continued 

and a maximum amount of diagnostic information can be furnished. Recovery procedures 

include supplying default or benign data, marking the data as erroneous so that other models 

do not attempt use of the data, or deleting the bad data. The recovery procedures sometimes 

generate additional diagnostic messages. Often after attempted correction of input, different 

diagnostic messages appear. These can be due to continued incorrect preparation of data, but 

103
the diagnostics may result from the more extensive testing permitted as previous errors are 

eliminated. 

The results of the input processing appear in the output file. If errors are found that 

cannot be recovered from, the code stops, with a message stating errors were found during 

input processing. More specific information on the type of error found and its possible 

location is also provided in the output file. If no input errors are found, the code prints a 

message in the output file stating that the input processing was successful. This information 

is the first part of the output file. Hence, when using the code, even if the code execution

appears to have been completed, before proceeding to extract output variables for plotting, it 

is a good idea for users to check the output file to verify that the input processing was 

successful. 

4.3.4 TRANSIENT PROCESSING 

Figure 4.3.2 shows the functional modular structure for the transient calculations, 

while Figure 4.3.3 shows the second-level structures for the transient/steady-state blocks or 

subroutines. The subroutine TRNCTL shown in Figure 4.3.1 consists only of the logic to call 

the next lower level routines. Subroutine TRNSET performs final cross-linking of 

information between data blocks, sets up arrays to control the sparse matrix solution,

establishes scratch work space, and returns unneeded computer memory. Subroutine TRAN, 

the driver, controls the transient advancement of the solution. Nearly all the execution time is

spent in this block, and this block is the most demanding of memory. Nearly all the dynamic

data blocks must be in the central memory, and the memory required for instruction storage 

is high, since coding to advance all models resides in this block. When transient advances are 

terminated, the subroutine TRNFIN releases space for the dynamic data blocks that are no 

longer needed. 

A description is next presented of the functions of all of the modules (subroutines) 

driven by TRAN (see Figure 4.3.3). The time step control module (DTSTEP) determines the

time step size, controls output editing, and determines whether the transient advancements 

should be terminated. During program execution, this module displays such information as 

CPU time, problem time, time step size, and advancement number on the standard output, 

usually a terminal screen. 

The trip system module (TRIP) evaluates logical statements. Each trip statement is a 

simple logical statement that has a true or false result. The decision of what action is needed

resides within the components in other modules. For example, valve components are 

provided that open or close the valve based on trip values. 

The equation of state boundary volume model (TSTATE) calculates the

thermodynamic state of the fluid in each hydrodynamic boundary volume (time-dependent 

volume). This subroutine also computes velocities for the time-dependent junctions. 

The heat structure model (HTADV) advances heat conduction/transfer solutions. It 

calculates heat transferred across solid boundaries of hydrodynamic volumes. This module 

104
will not be used if the user does not enter information in the input file regarding the presence 

of wall heat structures. Most FSP users will not use heat structures. 

The hydrodynamics module (HYDRO) advances the hydrodynamic solution. 

The reactor kinetics module (RKIN) will not be activated by FSP users It is a part of

the base code meant for the calculation of nuclear reactor kinetics. 

The control system module (CONVAR) provides the capability of simulating control 

systems typically used in hydrodynamic systems. It consists of several types of control 

components. Each component defines a control variable as a specific function of time

advanced quantities. These control variables are described in the code input manual in 

Appendix A. 

The time-advanced quantities include quantities from hydrodynamic volumes, 

junctions, valves, heat structures, trip quantities, and the control variables themselves. This

permits control variables to be developed from components that perform simple, basic 

operations. 

105
INPUT TRNCTL STRIP 

FSP 

Figure 4.3.1 Top Level Structure 

106
TIME 

STEP 

CONTROL 

EQUATION

OF STATE 

BOUNDARY

VOLUMES 

HYDRODYNAMICS CONTROL 

SYSTEM 

TRIP 

SYSTEM 

HEAT

STRUCTURE 

REACTOR 

KINETICS 

(Not Used 

in FSP) 

DRIVER 

Figure 4.3.2 Modular Structure of Transient Calculations 

107
TRNSET TRAN TRNFIN 

DTSTEP TSTATE HYDRO CONVAR 

TRIP HTADV RKIN 

TRNCTL 

Figure 4.3.3 Transient (Steady-State) Block Structure 

108
4.4 CODE INSTALLATION AND EXECUTION

This section describes how to install the FSP code on a personal computer with a 

WINDOWS 95 or 98 operating system, and how to execute the code. The code is transmitted 

on a CD-ROM with several files and folders in it. This section also describes what is on the

CD-ROM. 

In order to run a problem, the user needs to build an input file and name it with an 

extension ‘.i’. Instructions on how to build this file are provided in the code input manual,

shown in Appendix A. A sample problem is shown in Section 5 and can be used as a guide to 

building the input file. This file would be filname.i. The user can use as many characters as 

needed to adequately name this file, but generally keeping this to 8 characters is convenient, 

making sure to not include any blank spaces in the filename. A sample input file is contained

on the CD-ROM and can be located by searching for the file with the ‘.i’ extension. The

procedure for editing this file or for running the FSP code with this input file is described in 

sections 4.4.1 and 4.4.2 below. 

When the code is executed with the input file, filename.i, the code will create two 

files. One file is the output file, called filename.o. The user can examine this file by choosing 

the Edit option described in Section 4.4.1. A sample output file is contained in the CD-ROM. 

The second file is the restart file, called filename.r. This file cannot be examined on 

the screen because it is a binary file. However, it is the file needed for the user to extract 

information for plotting any of the code output variables. This file can also be used to restart 

the problem and run it for a longer transient time, but most problems are expected to be fairly 

small. Hence, it is recommended that the user not attempt to restart a problem, but simply 

rerun the problem from the beginning, to a longer transient time, as needed. 

To extract output variables for plotting, the user needs to build a file with the

extension, ‘.inp’, called filename.inp. This file is simply a listing of the variables that the

user wants to plot. The first line of this file is a title card, followed by one line for each 

variable desired. For example, if the user wanted to extract the pressure in component

number 203 as a function of time, one line of this file will read ‘p 203010000’. A sample is 

contained in the CD-ROM, and it is recommended that the user modify this file as needed.

Note that when this is done, the user should save the file with the same filename as used for 

the input file, but with the extension ‘.inp’. The procedure for editing this file is the same as

for editing any of the other files, as shown in Section 4.4.1. The procedure for extracting (or 

stripping) the output variables is shown in Section 4.4.3. When this stripping procedure is 

executed, the code will create a folder with the extension ‘.ace’, called filename.ace. Inside

this folder, the code will place a set of files, each file containing a set of X-Y pairs. Each of 

these files will be named with the name of the output variable requested. For example, the

file p.203010000, will contain a set of time and pressure values for component number 203 

which the user defined in the input file. 

To create a plot of any of the output variables chosen, the user may use any 

convenient software that the user is familiar with. The ACE graphics package, provided on 

the CD-ROM, offers one convenient method. The procedure for using this graphics software 

109
package is shown in Section 4.4.3. The user manual for this graphics package is contained on 

the CD-ROM.

4.4.1 SYSTEM REQUIREMENTS 

• Display settings should be at least 256 colors, 800x600 pixels (needed for graphics 

program Ace (aka xmgr) 

TCP/IP networking protocol must be installed for the Starnet X-Server software required

for the graphics program Ace (aka xmgr). See Start Menu – Settings – Control Panel – 

Network. 

• 

• WinZip is required for installation. (WinZip 7.0 Evaluation Copy is available as

winzip70.exe.) WinZip is available from www.winzip.com. 

The installation files require approximately 15 Mb on Drive C. This estimate includes 

2.5 Mb for WinZip and 5.5 Mb for the Starnet X-Server. 

• 

• Pentium II or III and Windows 95, 98 or NT Version 4. 

4.4.2 INSTALLATION 

Steps 1 through 5 create a scratch directory, load the necessary software, and unpack the 

FSP Version 1.0 archive in C:\fspv1. 

1. Check that the CDROM contains the following files. In the figure below, Windows 

Explorer was used to check the contents of the CD-ROM, which was loaded on Drive 

E. 

2. Create the C:\tmp directory required for execution (fsp, rstrip.) If input files will 

reside on another drive, then a \tmp directory must also be created on that drive. The 

\tmp directory is used for scratch files. 

110
3. This step can be skipped if X-Server software such as Hummingbird Exceed is 

already resident. Install the evaluation copy the Starnet x-win32 X-Server software

(needed for graphics program Ace (aka xmgr). Execute x-win40.exe from the

CDROM. Choose a typical installation and the default location. Enter “demo”

(without the “ ”) as the serial number when requested. X-win32 is available from 

www.starnet.com. 

4. This step can be skipped if WINZIP or a WINZIP compatible file utility program is 

already resident. Install the evaluation copy the WinZip file utility software (needed 

to decompress and load the fspv1.zip archive file). Execute winzip70.exe from the

CDROM. Choose a typical installation and the default location. Choose to have 

WinZip start with the “Classic interface” when asked. Note that WinZip is only 

needed for the next step in the installation. Once the installation has been 

satisfactorily completed, WinZip can be removed from the system via the Control 

Panel Add/Remove Programs dialog if desired. 

5. Use WinZip to extract all the files from fspv1.zip to C:\ using the path names in the 

WinZip archive (be sure to select the radio button for “All files” and check the 

checkbox for “Use folder names”).

Steps 6 through 10 create actions for specific file types (file extensions) in Windows

Explorer. These actions are then available on the pop-up menu that is launched by right 

clicking on a file. These actions can then be used to execute FSP and related programs

against the selected file.

Notes: 

When browsing for Wordpad in the subsequent steps, the Wordpad.exe is usually located 

in C:\Program Files\Accessories. 

In the steps below, when adding the “%1” in the subsequent steps, there should be a space 

after the .bat or .exe and before the “%1”. Also, the “ “ around the %1 are to be included 

because they protect pathnames that may include blank characters. 

6. Use Windows Explorer View Menu - Folder Options – File Types tab to setup actions 

for the new file type *.i files. 

File Description Action Command

FSP Input Edit Browse and select Wordpad then add 

“%1” 

 Run FSP C:\fspv1\fsp.bat “%1” 

7. Use Windows Explorer View Menu - Folder Options – File Types tab to setup an 

action for the new file type *.o files.

File Description Action Command

111
FSP Output Edit Browse and select Wordpad then add 

“%1” 

8. Use Windows Explorer View Menu - Folder Options – File Types tab to setup actions 

for the new file type *.inp files. 

File Description Action Command

Rstrip Input Edit Browse and select Wordpad then add 

“%1” 

Run Rstrip C:\fspv1\fsp.bat “%1” 

9. Use Windows Explorer View Menu - Folder Options – File Types tab to setup an 

action for the new file type *.out files. 

File Description Action Command

Rstrip Output Edit Browse and select Wordpad then add 

“%1” 

10. Use Windows Explorer View Menu - Folder Options – File Types tab to add the Run 

Ace action for the existing file type Folder. 

File Description Action Command

Folder Run Ace C:\fspv1\ace.bat “%1” 

4.4.3 FSP VERSION 1.0 FILES 

Table 4.1 lists the files in the WinZip archive including the date and time last 

modified, file size and pathname for installation. Installation Step 5 (above) uses this 

pathname to load these files in the proper directory structure. 

The .bat files supplied with FSP depend upon the directory structure created under 

C:\fspv1 by Installation Step 5. 

The directories C:\fspv1\fsp\test files\ and C:\fspv1\fsp\test files\r5mbg146.ace 

contain example inputs and outputs. These files are not necessary for FSP Version 1.0 

execution and may be deleted if so desired. All other files are needed for execution of FSP 

Version 1.0, Rstrip and the graphics program Ace (aka xmgr). 

112
4.4.4 CODE EXECUTION 

4.4.4.1 Editing a File 

The following file types may be viewed or edited with Wordpad: 

File Type Extension

FSP Input .i 

FSP Output .o 

Rstrip Input .inp 

Rstrip Output .out 

For example, view the test directory C:\fspv1\fsp\test files with Windows Explorer. Then 

select an input file such as R5mbg146.i with one left click of the mouse. Then, right click 

the selected file with the mouse. A pop-up menu will appear as shown below. This menu 

will include the actions created for *.i files in Installation Step 6. 

Selecting the Edit action will launch Wordpad to view or edit the file. 

113
4.4.4.2 FSP Execution 

To execute FSP bring up the pop-up menu for a FSP Input file (*.i) as described 

above and select the Run FSP action. A MSDOS Window will open and messages will be 

displayed as the \tmp\runit.bat file executes the FSP application as shown below. For long 

runs this Window may be minimized. After the run is completed, the Window name will 

indicate “Finished”. Then, the Window may by closed. 

Output and restart-plot files will be created in the same directory. These files will have the

same base name as the input file. The output file will have a .o extension and the restart-plot 

file will have a .r extension as shown below. 

114
4.4.4.3 Rstrip Execution 

To execute Rstrip bring up the pop-up menu for an Rstrip Input file (*.inp) as

described above and select the Run Rstrip action. A MSDOS Window will open and 

messages will be displayed as the \tmp\runit.bat file executes the Rstrip application as shown 

below. 

This Window may be minimized; however, the \tmp\runit.bat script prompts user before 

deleting all the files in a directory. If the Window is minimized the prompt will not be seen

and the \tmp\runit.bat script will remain waiting for a response. 

Answer the “Are you sure (Y/N)?” prompt in the MSDOS Window with a y or Y and Enter

if the directory name indicated is of the form \tmp\ace##\.

After the run is completed, the Window name will indicate “Finished”. Then, the Window 

may by closed. 

An output with the same base name as the input file will be created in the same directory. 

The output file will have a .out extension. 

A sub-directory with the same base name as the input file will be created in the same

directory. This sub-directory will have a .ace extension. The sub-directory will contain a file 

for each minor edit variable stripped from the restart-plot file as shown below. The names of 

the individual files will be formed using the FSP minor edit variable name and the FSP 

notation for volume, junction or control variable number. For example file p.203010000 will 

contain the time and pressure values for volume 20301. 

115
4.4.4.3 Ace (aka xmgr) Execution 

To execute the Ace graphics package bring up the pop-up menu for an (*.ace) 

directory as described previously and select the Run Ace action. A MSDOS Window will

open and messages will be displayed as the \fspv1\ace.bat file brings up x-win32 X-Server as 

shown below. 

Dismiss the X-Win32 Demo message by clicking on the OK button. 

116
Ace will then display an empty plot area. Use the File menu and select Read Sets. A Read 

Sets file selection dialog will then pop-up. Left Click on one of the minor edit data files

shown in the Read Sets file selection dialog and then left click on the OK button. A plot of

the data will then be created as shown below. 

When finished looking at plots of the FSP results, use the File menu and choose Exit to close

the Ace graphics program. Also close the X-Win32 X-Server that runs minimized by right 

clicking on its icon in the taskbar and choosing Close. 

Table 4.4.1 - Contents of the WinZip Archive fspv1.zip 

117

File Name Modified Size Path

Ace.bat 1/19/0010:50 PM 326fspv1\

Ace.pif 1/20/0001:48 AM 967fspv1\

Fsp.bat 1/19/0011:11 PM 212fspv1\

Fsp.PIF 1/20/0001:46 AM 967fspv1\

Zzz.exe 1/19/0011:06 PM 340,011fspv1\

Tpfh2o 1/18/0011:09 PM 129,440fspv1\fsp\envrl\

Fsp.exe 1/20/0001:35 AM 2,105,344fspv1\fsp\fsp_mod\Release\ 
File Name Modified Size Path

Tpfco2 1/19/0010:02 PM 62,344fspv1\fsp\nf\

Tpfhalon 1/19/0010:02 PM 85,136fspv1\fsp\nf\

Tpfr125 1/19/0010:02 PM 117,544fspv1\fsp\nf\

TpfHFC-227ea 1/19/0010:03 PM 100,048fspv1\fsp\nf\

R5mbg146.i 1/19/0011:54 PM 8,194fspv1\fsp\test files\ 

R5mbg146.inp 1/19/0011:52 PM 164fspv1\fsp\test files\ 

R5mbg146.o 1/20/0001:46 AM 263,305fspv1\fsp\test files\ 

R5mbg146.out 1/20/0001:46 AM 937fspv1\fsp\test files\ 

R5mbg146.r 1/20/0001:46 AM 521,368fspv1\fsp\test files\ 

R5mbg146_hp.o 1/19/0010:46 AM 266,429fspv1\fsp\test files\ 

Chokej.204000000 1/20/0001:46 AM 6,529fspv1\fsp\test files\r5mbg146.ace\ 

Florgj.204000000 1/20/0001:46 AM 6,529fspv1\fsp\test files\r5mbg146.ace\ 

Mflowj.204000000 1/20/0001:46 AM 6,538fspv1\fsp\test files\r5mbg146.ace\ 

p.203010000 1/20/0001:46 AM 6,524fspv1\fsp\test files\r5mbg146.ace\ 

Velgj.204000000 1/20/0001:46 AM 6,535fspv1\fsp\test files\r5mbg146.ace\ 

Voidgj.204000000 1/20/0001:46 AM 6,529fspv1\fsp\test files\r5mbg146.ace\ 

Rstrip.exe 1/20/0012:27 AM 319,488fspv1\rstrip\

C 9/22/9711:21 AM 3,184fspv1\xmgr-exceed\bin\

Hcabout.dll 9/22/9711:21 AM 153,600fspv1\xmgr-exceed\bin\

Hcaboutn.nls 9/22/9711:21 AM 30,208fspv1\xmgr-exceed\bin\

Hclctl.dll 9/22/9711:21 AM 18,944fspv1\xmgr-exceed\bin\

Hclctl3d.dll 9/22/9711:21 AM 26,112fspv1\xmgr-exceed\bin\

Hcllocal.dll 9/22/9711:21 AM 12,288fspv1\xmgr-exceed\bin\

Hclnls.dll 9/22/9711:21 AM 8,192fspv1\xmgr-exceed\bin\

Hclt.dll 9/22/9711:21 AM 6,656fspv1\xmgr-exceed\bin\

Hclwprod.dll 9/22/9711:21 AM 18,944fspv1\xmgr-exceed\bin\

Hclwutil.dll 9/22/9711:21 AM 31,232fspv1\xmgr-exceed\bin\

Hclxpmon.dll 9/22/9711:21 AM 47,616fspv1\xmgr-exceed\bin\

Hclxport.dll 9/22/9711:21 AM 39,424fspv1\xmgr-exceed\bin\

Nls.dir 9/22/9711:21 AM 1,788fspv1\xmgr-exceed\bin\

Pout.dat 3/19/9809:27 PM 12,216fspv1\xmgr-exceed\bin\

Xerrordb 9/22/9711:21 AM 8,604fspv1\xmgr-exceed\bin\

XkeysymDB 12/21/9909:35 PM 4,696fspv1\xmgr-exceed\bin\

Xlib.dll 9/22/9711:21 AM 279,040fspv1\xmgr-exceed\bin\

Xmgr.exe 9/21/9710:14 PM 1,785,856fspv1\xmgr-exceed\bin\

xt.dll 9/22/9711:21 AM 185,856fspv1\xmgr-exceed\bin\

Number of files:43 

Total 

Size: 7,035,864Mb

118
4.5 INPUT PREPARATION AND SAMPLE PROBLEM 

This section describes input preparation and code execution for a typical sample

problem. The Lehigh test configuration is chosen for simulation. The information presented 

here, along with the input manual shown in Appendix A, can be used to understand how to 

develop an input model of a desired suppressant delivery system using the FSP code. The

intent of this section is to describe how the user can set up an input file, and examine relevant 

output variables to analyze the performance of a suppressant delivery system. 

4.5.1 System Configuration 

The Lehigh University experimental set up consists of a test section made of a 200-

cm long 2” Schedule 80 commercial steel pipe. The test section is considerably longer than 

typical supply vessels, but the specific geometry is chosen to enable an accurate 

measurement of the liquid level inside the test section, and subsequent determination of the 

discharge mass flow rate from the test section. The bottom of the test section is equipped 

with a Marotta quick-opening valve. The discharged fluid is piped through a 11-ft long 1⁄2” 

Schedule 80 steel pipe to a recovery tank. The recovery tank is about 10” in diameter and

has a volume of approximately 28 liters. The pressure in the storage tank is expected to

range from 1 to 3 atm at room temperature. Cooling tubes are installed inside the storage 

tank to achieve the desired downstream pressure and temperature between tests and also to 

condense vapor and return the liquid agent for subsequent tests. The top of the test section is

connected to a nitrogen bottle, for initial pressurization. Figure 4.5.1 shows the schematic of

the experimental set-up. Liquid nominally occupies 60% of the test section. 

4.5.2 Development of FSP Model 

In developing the input model, the physical system is typically discretized into a 

series of cells or volumes (called SNGLVOL components). These cells or volumes are

interconnected via junctions (called SNGLJUN components). The physical dimensions of 

the system are specified via the geometrical input data associated with these volumes and 

junctions. Initial fluid thermodynamic conditions are specified for each cell, and initial flow 

conditions are specified for each junction. Appendix A provides the details of the required

input data. To further facilitate the development of an input model, FSP has other built-in 

special components that can be used to represent typical piping systems. These components 

include PIPE, VALVE, and BRANCH components. A PIPE component is simply a

collection of a number of cells connected by junctions, and is most appropriate to model a

physical piping section. This is more convenient to use than specifying a number of 

SNGLVOL and SNGLJUN components to represent each cell and flow path inside the pipe

section. 

For the sample problem described here (i.e., the Lehigh apparatus), the vertical test

section is represented by a PIPE component 100 with 10 cells. The horizontal discharge pipe 

is modeled by another PIPE component 300 with 15 cells. A VALVE component (203) is

used to model the Marotta valve connecting the bottom of the test section to the discharge

119
pipe. The collection vessel is modeled by a SNGLVOL component 500. Finally, the 

connection between the end of the discharge pipe and the collection vessel is represented by

a SNGLJUN component 403. Note that the outlet of the system in this case is the storage 

tank, and it has been modeled using a SNGLVOL component. This is because, in this

experiment, the pressure in the storage tank would increase during the experiment, and at 

some point would have an effect on the system due to its back pressure on the system. This 

model tries to account for this effect and hence models the storage tank as an ‘active’ 

volume. In a design simulation using the FSP code, the outlet would discharge into a large 

room or the atmosphere, whose pressure would not be affected by the discharge of the agent.

To model that situation, the user should use a TMDPVOL component to represent the outlet 

volume, set at a constant pressure. 

The FSP nodalization diagram is shown in Figure 4.5.2. The following presents the 

details of input data development for the various components used to represent the Lehigh 

apparatus. The major reference for this section is the Crane Technical Paper No. 410, 24th

Printing, (1988). The major components of the experimental set-up are modeled as follows. 

Test Section: PIPE componet 100 

The test section is a 2-m long 2” Schedule 80 commercial steel pipe. In the particular test 

simulated (Test 3 from Configuration 1), 61.9% of the test section is occupied by liquid agent 

(i.e., initial liquid level is 1.238 m from test section bottom). From Crane (1988), the ID of 

the test section is 

ID = 1.939” = 0.04925 m, from Crane, page B-16 

Refer to Appendix A, the required input data for this component are as follows (all input are 

in SI units per Appendix A unless otherwise specified) 

Component No.= 100 

No. of Volumes = 10 

Volume Area = π(ID)2

/4 = π(0.04925)2/4 = 1.9050e-3 

Junction Flow Area = 1.905e-3 or default 

Volume Length 

In order to model the initial liquid height exactly, the section is discretized such that the top 4

cells have a total length of 0.762 m, and the bottom 6 cells 1.238 m. Hence, the volume

length input data are given by: 0.200 for cells 1 to 3; 0.162 for cell 4; 0.238 for cell 5; and 

0.200 for cells 6 to 10. Note that the test section is oriented downwards, with cell No. 1 at

the top and cell No. 10 at the bottom.

Other Geometrical Input 

Volume of Volume = 0.0, default 

Vertical Angle = -90o

Elevation Change = Volume Length x sin(Vertical Angle) = negative of the volume lengths 

Roughness = 1.5e-4 ft = 4.572e-5 for commercial steel, from Crane, page A-23 

Hydraulic Diameter = ID = 0.04925 or default 

Forward Loss Coefficient = 0.0, straight pipe 

120
Reverse Loss Coefficient = 0.0, straight pipe 

Volume Control Flag = 00000 

Junction Control Flag = 000000 

Volume Initial Condition Option = 04 for Volumes 1 – 4, and 13 for remaining volumes 

Initial Conditions: 

As mentioned before, the test section is initially 61.9% agent filled. The test section is

subdivided into 10 cells with the bottom 6 cells filled with liquid, and the top 4 cells filled

with a mixture of agent vapor and nitrogen. Thus, for Volumes 1-4, the initial condition

option No. 4 (see Appendix A) is used to accommodate the presence of noncondensable gas 

(i.e., nitrogen). The theoretical equilibrium gaseous (vapor + noncondensable) quality is 1.0 

(since no liquid is present). For reasons attributed to the numerical solution method of the 

code, however, it is recommended that a value close to, but less than 1.0 be used. Using a 

value of 1.0 has led to the code grinding down to small time steps and in some cases, stop

running altogether. Experience shows that a value of 0.98 mitigates the problem, without 

affecting the physics of the model, and is hence used in this sample problem. 

For Volumes 5-10, only liquid (with dissolved nitrogen) is present. Hence initial condition

option No. 13 is used which allows the user to specify the pressure and temperature, and

account for dissolved nitrogen. For this sample problem, the initial filled pressure is 4.180 

MPa, and the initial system temperature is 305 K.

Junction Initial Condition Option = 1, specify mass flow 

Junction Initial Conditions: zero flow 

The initial concentration (based on mass) of dissolved nitrogen is calculated as follows. 

Initial Concentration of Dissolved Nitrogen 

As discussed in Section 4.2.4.6, FSP uses the total pressure in Henry’s law to calculate the 

mole fraction of the dissolved nitrogen in the liquid solution. Theoretically, the partial 

pressure of nitrogen should have been used, but it poses difficulty for numerical 

implementation. To be consistent with this approximation, the initial equilibrium

concentration of dissolved nitrogen should also be calculated using the total pressure. Thus, 

the initial mole fraction of dissolved nitrogen is calculated as

xN,i = HPi

where H is the Henry’s constant given in Section 4.2.4.6 and Pi is the initial system pressure. 

The required mass concentration (mass of dissolved nitrogen per unit mass of liquid) for 

input is: 

X = xN,i MN/(1-xN,i)Ma = HPi MN/(1-HPi)Ma

where MN and Ma are the molecular weight of nitrogen and liquid agent. The liquid agent in 

this case is HFC-227ea, which has the formula, C3HF7. So, in the sample problem, this is

evaluated as: 

121
H = 2.347767 x 10-7 –1.55063 x 10-9 (305) + 2.957799 x 10-12

 (305)2 = 3.698 x 10-8

X = (3.698 x 10-8)(4180000) 28 / [1-(3.698 x 10-8)(4180000)](170) = 0.03012 

Since dissolved gas is only present in the liquid, the initial dissolved gas concentration is

entered as 0.0 for Volumes 1 through 4, and as 0.03012 for Volumes 5 through 10. 

Discharge Valve 

The Marotta valve is a quick opening valve (with a stroke time of about 0.05 s). Discharge 

data for this valve were obtained by experiments at Lehigh University. The data show that 

the equivalent loss coefficient associated with the valve is about 1.5. Refer to Appendix A, 

the required input data are as follows: 

Component No.= 203 

Component Type = valve 

From Connection Code = 100010000, outlet of Component 100 

To Connection Code = 300000000, inlet of Component 300 (discharge pipe to be discussed 

later) 

Area = default (i.e., minimum of adjacent volume areas) 

Forward Loss Coefficient = 1.5, based on test data supplied by Lehigh University. 

Reverse Loss Coefficient = 1.5. 

Junction Flag = 000000 

Discharge Coefficients = 1.0 for subcooled, 0.9 for both two-phase and single-phase vapor. 

These values have been determined to be optimum during the data assessment stage. 

Initial Condition Option = 0, specify velocity 

Initial Conditions = zero velocity 

Valve Type = trpvlv, for valve with negligible stroke time

Open Trip No. = 501, where Trip 501 is true when t>t1 where t1 is the steady state

initialization period. A short steady state initialization period is recommended so that the 

system can adjust itself for hydrostatic pressure balance prior to opening of the valve. In this 

sample problem, a steady-state initialization period of 0.5 seconds has been allowed. 

Discharge Pipe 

The discharge pipe is a 11 ft long 1⁄2” Schedule 80 commercial steel pipe. Hence, from Crane 

(1988), the ID is given by 

ID = 0.546” = 1.3868e-2 m. from Crane, Page B-16 

Refer to Appendix A, the required input data are as follows: 

Component No. = 300 

Component Type = pipe 

Number of Volumes = 15, to yield approximately 0.2-m grid size 

Volume Area = π(ID)2

/4 = π(0.013868)2/4 = 1.5105e-4 

Junction Flow Area = 1.5105e-4, or default 

Volume Length = 3.3528/15 = 0.2235 

122
Volume Volume = 0.0, default 

Vertical Angle = 0o

, horizontal, from Figure 5.1.1 

Elevation Change = Volume Length x sin(Vertical Angle) = 0.0 

Roughness = 1.5e-4 ft = 4.572e-5 for commercial steel, from Crane, page A-23 

Hydraulic Diameter = ID = 0.013868 or default 

Forward Loss Coefficient = 0.0, straight pipe 

Reverse Loss Coefficient = 0.0, straight pipe 

Volume Control Flag = 00000 

Junction Control Flag = 000000 

Initial Conditions: 

Normally, the discharge pipe is predominantly air-filled at the collection vessel pressure, 

prior to the discharge test. Hence, the initial condition option No. 4 should be used. The 

quality is theoretically 1.0 (no liquid present). As discussed earlier, 0.98 is used in place of 

1.0. 

Volume Initial Condition Option = 4 

Volume Initial Conditions: P = 7.200e5 Pa, T=305 K, Quality=0.98 

In case that the back pressure and initial temperature do not allow the specification of a

mixture of vapor and noncondensable, it is recommended that superheated vapor condition

be used. For instance, if the back pressure is greater than the saturation pressure at the initial 

temperature, it would not be feasible to specify a mixture of vapor and noncondensable using 

Option No. 4. In that case, the given conditions will support the specification of single-phase 

superheated vapor. The user should then use Option No. 3 and input the given back pressure 

and temperature. 

It should be mentioned that allowing air (or noncondensible) to be present downstream of the 

discharge valve causes computation difficulties occasionally. By replacing the air-filled

condition by superheated vapor at the same pressure and temperature (if possible), the code’s 

performance is decidedly improved. In situations where the specified back pressure and 

initial temperature do not allow the existence of single phase vapor, the back pressure may be 

adjusted accordingly to permit the existence of single phase vapor. Since the flow is 

expected to be choked at the discharge valve, the actual magnitude of the back pressure 

should have only minimal impact on the discharge transient. The user may need to 

experiment with both options 3 and 4 when computational difficulties are encountered. 

Junction Initial Condition Option = 1, specify mass flow rate 

Junction Initial Conditions = zero flow 

Inlet to Storage Tank 

The other end of the discharge pipe is connected to a storage tank, as shown in Figure 5.1.1. 

Refer to Appendix A, the required input data are as follows: 

Component No. = 403 

Component Type = sngljun, for single junction 

From Connection Code = 300010000, outlet of Component 300 

123
To Connection Code = 500000000, inlet of Component 500 (storage tank to be discussed 

later) 

Junction Area = 1.5105e-4, same as discharge pipe area 

Forward Loss Coefficient = 1.0, sudden expansion, from Crane, Page A-29 

Reverse Loss Coefficient = 0.5, assume sharp-edged, from Crane, Page A-29 

Junction Flag = 000000 

Discharge Coefficients = 1.0 for subcooled, 0.9 for two-phase and single-phase vapor. These 

values have been determined during the assessment stage. 

Initial Condition Option = 1, specify mass flow rates 

Initial Conditions = zero flow 

Storage Tank

The storage tank is approximately 10” in diameter and has a volume of about 28 liters (or 

0.028 m3

). During a test run, its pressure will increase due to the build up of discharge fluid.

Between tests, the tank and its content may be cooled to bring the system back down to a

desirable pressure and temperature level. Refer to Appendix A, the required input data are as 

follows: 

Component No. = 500 

Component Type = snglvol, for single volume 

Volume Flow Area = π/4 (ID)2

 = π/4 (0.254)2 = 0.05067 

Length of Volume = 0.0, default 

Volume of Volume = 28 liters = 0.028 

Inclination Angle = -90o

Elevation Change = - (Volume/Area) = -(0.028/0.05067) = -0.5525 

Roughness = 1.5e-4 ft = 4.572e-5, assume commercial steel, from Crane, page A-23 

Hydraulic Diameter = 0.0, default 

Control Flag = 00000 

Initial Conditions 

The discussion of the initial conditions in the discharge pipe applies here. 

Initial Condition Option = 4 

Initial Conditions: P = 7.200e5 Pa, T=305 K, Quality=0.98 

Note that for numerical stability, the user may need to approximate the storage tank 

conditions as single-phase vapor at the desired pressure and temperature, as discussed before. 

Nitrogen Release Model 

The nitrogen release model has been discussed in Section 4.2.4.6. It requires four user input 

data: the release coefficient (Γ0), critical radius (rc), initial pressure (Pi) and the initial system

temperature (Ti). Note that Pi and Ti are the conditions just prior to discharge. They may be 

different from the so-called fill pressure and fill temperature. To evaluate Pi and Ti given the 

fill values, a code such as PROFISSY (Yang et al) could be used. Both the release

coefficient and critical radius are semi-empirical. Section 2.4.6 has presented some

sensitivity results. For the sample problem here, the following input data are specified: 

124
Γ0 = 2000 

rc = 10-8

Pi = 4.180 x 106

, same as the initial system pressure 

Ti = 305, same as the initial system temperature 

Trip Logic 

The code allows the user to specify trip logic to perform certain logical operation, such as 

opening or closing a valve. In this problem, the Marotta valve is modeled as a trip valve with 

negligible stroke time. The valve is initially close. It opens at 0.5 s (0.5 s is a steady state 

initialization period). This behavior can be modeled in FSP by specifying one trip that opens 

at 0.5 s. For time less than 0.5 s, the trip is false, and the valve will remain at its initial 

specified position (which is close for this application). Refer to Appendix A, the following 

input data for the trips are used: 

Trip 501: t ≥ 0.5. Open Trip. This trip is true when time is greater than 0.5 s. 

Control Variables 

Control variables (CV) are used to monitor certain system responses, such as the collapsed 

liquid level in the test section and the integral discharge flow mass. This section documents

the development of these control variables. 

Collapsed Liquid Level: CV 010 

CV010 = Σ Vol. Lengthi (voidf)i, i = 100010000 through 100100000 

Integrated Discharge Mass from bottle: CV 020 

CV020 = ∫ mflowj @ 203000000 dt 

Integrated Discharge Mass from discharge pipe: CV 021 

CV021 = ∫ mflowj @ 403000000 dt 

4.5.3. Input File Listing

In this sample problem, we have named the input file as run1-3.i. A listing of run1-

3.i is given below. 

= Lehigh Experiment, Configuration 1, Run 3 

* 

***************************************************** 

* 

* Test Conditions 

125
* 

* Fluid = HFC-227ea 

* Initial Pressure = 4180 kPa 

* Initial Temperature = 305.0 Deg.K. 

* Initial Liquid Volume = 2358 ml 

* Initial Liquid Height = 1.238 m 

* (Based on x-section area of 1.9053e-03 m2) 

* Initial Downstream Pressure = 720.0 kPa 

***************************************************** 

* 

* Model Options 

* 

* CD (two-phase) = 0.9 

* Gamma0 = 2000.0 

* Crit. Radius = 1.00e-8 m. 

* 

***************************************************** 

* 

100 new transnt 

101 run 

102 si si 

* 

* Gamma0 CritRad Pfill Tfill 

* 

106 0.0 1.00e-08 4.180e06 305.0 

* 

110 nitrogen 

* 

* RefVol Elev Fluid System 

* 

120 100010000 0.0 HFC-227ea testsec 

* 

* Time Step Control 

* 

* End Min Req Flags Minor Major Restart 

* 

201 0.5 1.0e-6 0.001 00003 100 500 500 

202 3.5 1.0e-6 0.0001 00003 500 5000 5000 

* 

* Minor Edit Variables Requested 

* None requested, Cards 301-399 not necessary 

* 

* Trips 

* 

501 time 0 gt null 0 0.5 l *open discharge valve at 0.5 

secs. 

502 time 0 lt null 0 0.0 l *always false 

* 

****************************************************************** 

* Test Section, Component 100 

****************************************************************** 

* Comp.Name Type 

1000000 testsec pipe 

* nv 

1000001 10 

* Vol.Area nv 

1000101 1.9053e-3 10 

* Jun.Area nj 

1000201 0.0 9 

126
* Vol.Len nv 

1000301 0.200 03 

1000302 0.162 04 

1000303 0.238 05 

1000304 0.200 10 

* Vol.Vol nv 

1000401 0.0 10 

* Vert.Ang nv 

1000601 -90.0 10 

* Elv.Chn nv 

1000701 -0.200 03 

1000702 -0.162 04 

1000703 -0.238 05 

1000704 -0.200 10 

* Rough Dh nv 

1000801 4.57e-5 0.0 10 

* kf kr nj 

1000901 0.0 0.0 9 

* Vol.Flag nv 

1001001 00000 10 

* Jun.Flag nj 

1001101 001000 9 *no choking option 

* IC Pres Temp Qual Dummy nv 

1001201 004 4.180e6 305.0 0.98 0.0 0.0 4 

1001202 013 4.180e6 305.0 0.0 0.0 0.0 10 

* Vel(0)/Flow(1) 

1001300 1 

* Liq Vap Int nj 

1001301 0.0 0.0 0.0 9 

* 

* Solute Concentration, KgN2/KgLiq 

* 

1002001 0.0 04 

1002002 0.03012 10 

* 

******************************************************************** 

* Marotta discharge Valve, Component 203 

******************************************************************** 

* Comp.Name Type 

2030000 break valve 

* From To Area kf kr Flag Cds 

* Liq 2p 

Vap 

2030101 100010000 300000000 1.5105e-4 1.5 1.5 000100 1.00 0.9 0.9 

* IC Liq Vap Int 

2030201 0 0.0 0.0 0.0 

* Valve Type 

2030300 trpvlv 

* Trip No. 

2030301 501 

* 

******************************************************************** 

* Discharge Pipe, Component 300 

******************************************************************** 

* Comp.Name Type 

3000000 dispipe pipe 

* nv 

3000001 15 

* Vol.Area nv 

127
3000101 1.5105e-4 15 

* Jun.Area nj 

3000201 0.0 14 

* Vol.Len nv 

3000301 0.2235 15 

* Vol.Vol nv 

3000401 0.0 15 

* Vert.Ang nv 

3000601 0.0 15 

* Elv.Chn nv 

3000701 0.0 15 

* Rough Dh nv 

3000801 4.57e-5 0.0 15 

* kf kr nj 

3000901 0.0 0.0 14 

* Vol.Flag nv 

3001001 00000 15 

* Jun.Flag nj 

3001101 001000 14 *no choking option 

* IC Pres Temp Qual Dummy Dummy nv 

3001201 004 7.200e5 305.0 0.98 0.0 0.0 15 

* Vel(0)/Flow(1) 

3001300 1 

* Liq Vap Int nj 

3001301 0.0 0.0 0.0 14 

******************************************************************** 

* Inlet to Storage Tank, Component 403 

******************************************************************** 

* Comp.Name Type 

4030000 inlet sngljun 

* From To Area kf kr Flag Cds 

* Liq 2p Vap 

4030101 300010000 500000000 1.5105e-4 1.0 0.5 000000 1.00 0.9 0.9 

* IC Liq Vap Int 

4030201 1 0.0 0.0 0.0 

******************************************************************** 

* Storage Tank, Component 500 

******************************************************************** 

5000000 storage snglvol 

* Area Len Vol Azi Incl Elev Rough Dh 

Flag 

5000101 0.05067 0.0 0.028 0.0 -90.0 -0.5526 4.572e-5 0.0 00000 

* IC Pres Temp Qual 

5000200 4 7.200e5 305.0 0.98 

******************************************************************** 

* 

** Control Systems 

* 

* Collapsed Liquid Level in Test Section (m) 

* 

20501000 level sum 1.00 0.0 1 

20501001 0.0 0.2 voidf 100010000 

20501002 0.2 voidf 100020000 

20501003 0.2 voidf 100030000 

20501004 0.162 voidf 100040000 

20501005 0.238 voidf 100050000 

20501006 0.2 voidf 100060000 

20501007 0.2 voidf 100070000 

128

20501008 0.2 voidf 100080000 
20501009 0.2 voidf 100090000 

20501010 0.2 voidf 100100000 

* 

* Integral discharge mass from bottle (kg) 

* 

20502000 intbrk integral 1.00 0.0 1 

20502001 mflowj 203000000 

* 

* Integral discharge mass from pipe (kg) 

* 

20502100 intpipe integral 1.00 0.0 1 

20502101 mflowj 403000000 

* 

* 

. end of input listing 

129
4.5.4 Results & Data Comparison

A new case and keep track of the changes separately 

The input file was named run1-3.i. The output from FSP consists of an ASCII output 

file (called run1-3.o) and a BINARY restart/plot file (run1-3.r). The ASCII output file 

records the major thermal hydraulic variables at time interval specified by the user through 

the Major Edit Frequency. The BINARY restart/plot file records most of the calculated

variables at time interval specified by the user through the Restart Edit Frequency. Refer to 

Appendix A for details. The ASCII output file is useful in examining the conditions of the 

entire system at a given time. The BINARY restart/plot file is useful in examining the 

transient response of a particular variable. For general data comparisons, the restart/plot file 

is used. An extraction program has been provided to allow the extraction of any valid 

variable from the restart file and has been described in Section 4. A list of the valid variables 

in FSP is given in Appendix A. Once the variable is extracted, it can be plotted with or 

without the corresponding data for examination. 

As discussed in Section 4, the first step is to create an input file containing the list of 

the variables the user wants to extract from the restart file for plotting. This file should be 

named run1-3.inp. In this example, we have chosen to examine the test section pressure, the 

collection vessel pressure, and the discharge mass flow rate. The following is the listing of

run1-3.inp

rst01 c:\mydocu~1\fire-s~1\Lehigh~1\run1-3.r

p 100010000 

p 500010000 

mflowj 203000000 

Note that the first input line (the title card) specifies the full DOS path of the restart/plot file, 

and renames it to a file named “rst01”. The user may not modify the name “rst01”. 

Upon execution of the Rstrip program (see Section 4.4.4.3), a directory run1-3.ace will be

created under the current working directory. The following three files will be created by 

Rstrip in the run1-3.ace directory: p.100010000, p.500010000, and mflowj.203000000. 

Each of these files is simply the transient history of each variable (i.e, variable versus time). 

If data comparison is desired, the corresponding data as a function of time must be created by 

the user. The format of any data file should be the same as the FSP data file, i.e., each data 

line contains time-data pair. There is no restriction to the naming of these user-created data 

files. In this example, we have created the following test data files: pbot.dat (measured test

section pressure), and pcv.dat (measured collection vessel pressure). It is recommended that 

the data files be placed in the same directory as the FSP data files. In this example, pbot.dat

and pcv.dat are placed in the run1-3.ace directory. 

Next, the Ace program can be executed to plot any of the FSP data file created by Rstrip and

any experimental data file created by the user. The execution instruction has been given in 

Section 4.4.4.4. For instance, the compare the calculated and measured test section pressure,

follow Section 4.4.4.4 instruction to open the Ace window. Use the File menu and select 

130
Read Sets. A Read Sets file selection dialog window will pop-up. Left Click on 

p.100010000 and then left click OK. Left click on the AS button on the left menu bar to 

choose Automatic Scaling. A plot of the calculated pressure in cell 100-01 (i.e, the top of the 

test section) will appear. Return to the Read Sets file selection dialog window, Left Click on 

pbot.dat and then left click OK. The measured test section pressure will appear in the same

plot. At this point, the user can add title and legend, change the line style of each curve, etc. 

The Ace graphic program allows the user to perform many tasks to customize each plot. The 

user should consult the Ace manual for instructions, or experiment with the various action 

items under the File, Data, Graph, Page and View menus. 

When finished looking at the plots, use the File menu and choose Exit to close the Ace

graphic program. The user can route the plot to a printer that supports postscript printing. 

Alternatively, the user can write the plot file to a postscript file. To generate the postscript 

file, use the File menu and choose Printer setup. A Printer setup window will pop up. 

Choose Postscript landscape in the “Device” dialog, File in the “Print to” dialog, type in the 

name of the postscript file to be created in the “Print Control String” box (by convention, this 

file should have extension ‘.eps’) and depress the “generate EPS” button. Left click the Print 

button. A file with the specified name will be created. This postscript file can then be 

imported to a document, or printed out later. 

As an example, Figure 4.5.3 shows the comparison of the calculated and measured test 

section pressure for the sample problem, obtained using the above procedure. Looking at the 

results of the comparison, the user may want to conduct a sensitivity study by varying an 

input variable, such as rc or Γ0. In this case, the user would change this parameter in the input 

file and re-run the case. It would be prudent to rename the input file when running.

131
Safety Relief Valve 

Level 

Indicator 

200 cm

N2

Test Section

2” Sch.80

Storage Tank

Discharge Valve 

Assembly 10” dia 

Discharge Pipe (1/2” Sch.80) 

11 ft

Figure 4.5.1: Schematic of Lehigh Test Apparatus 

132
PIPE 300 

VALVE 203

PIPE 100 Source Vessel 

SNGLVOL 

500 

SNGLJUN 403

Discharge Pipe

Valve 

Assembly 

Collection Vessel 

Figure 4.5.2: FSP Nodalization Diagram for Lehigh Tests 

133
Figure 4.5.3: Comparison of Calculated and Measured Bottle Pressure for Sample Problem

134
4.6 ASSESSMENT AGAINST EXPERIMENTS 

This section discusses some of the code assessments carried out against experiments.

Early assessment of the FSP code was against Halon-1301 discharge experiments reported by 

Elliott et al (1984). These tests were carried out with three different system configurations. 

These included a bottle-nozzle arrangement, a bottle with straight pipe, and a bottle with a

straight pipe branching to two separate exit legs. Test 146 from Elliott et al (1984) was with a 

bottle-nozzle arrangement. This test was used in early assessment of FSP. The results have 

been shown in Section 4.2.4.6 of this report, in the discusssion of the gas release model. The 

other configurations tested by Elliott et al (1984) were also simulated during various 

developmental stages of the code. The results generally showed good agreement between 

code calculations and experimental data. A re-assessment of the final code version against 

data from these configurations was not carried out. This was primarily because the tests were

conducted with Halon-1301, and the intent of code development was for alternate fluids. The

expectation is that the final version of FSP would be able to simulate these tests quite well. 

The experiments at Lehigh University, carried out as a part of this code development

project, examined several piping configurations, with HFC-227ea and with HFC-125. The 

code has been assessed against these experiments. These experiments are particularly useful

for code assessment because, in addition to the usual measurement of pressure, the

experiments were able to estimate the mass flow rate of discharge, and also measured the 

void fraction near the pipe exit. The results of the assessment are being presented separately, 

and generally show good agreement against the test data. 

Other experimental data on suppressant discharge through typical delivery systems have been 

compiled and reported by Grosshandler et al, 1994, and by Gann, 1995. This section presents

the assessment carried out against some of the experiments reported in Gann, 1995. The

experiments were conducted at NIST, and included a series of discharge tests using several 

agents (including Halon 1301). 

4.6.1 FSP ASSESSMENT AGAINST NIST EXPERIMENTS 

4.6.1.1 NIST Experimental Apparatus 

The basic design of the NIST experimental apparatus consisted of a storage vessel (or 

pressure vessel), discharge piping, and a recovery tank. A nitrogen make-up tank was

connected to the pressure vessel to simulate constant-head discharge if needed. The pressure 

vessel was fitted with a quick release valve at the bottom. Figure 4.6.1 shows the test 

apparatus. The data included pressure vessel pressure and pressure at various locations along 

the discharge piping. Several discharge piping configuration were used. These included: (1) 

3.5-m straight discharge; (2) 3.5-m with 90-degree bend; (3) 4.5-m with 90-degree bend; (4) 

3.5-m straight with expansion from 9.5-mm ID to 15.9-mm ID; (5) 3.5-m straight with 

contraction from 15.9-mm ID to 9.5-mm ID; and (6) 4.5-m, 15.9-mm ID Tee configuration.

Two pipe sizes were used for all configuration, except otherwise stated, and these were 9.5-

mm ID and 15.9-mm ID. The various configurations are shown in Figure 4.6.2. 

135 
Four agents were tested in these experiments. These were Halon 1301, HFC-125,

HFC-227ea and CF3I. The assessment presented here is only for the HFC-227ea (or HFC227ea) experiments. The Table below lists the tests conducted with HFC-227ea. Note that the 

test identifier in the first column is actually the Figure Number in Appendix II of the NIST 

Report (Gann, et al, 1995). The various configurations shown in the last column are shown in

Figure 4.6.2. 

Test 

Identifier 

Pressure 

(Mpa) 

Fill Volume

(%) 

Configuration 

B-54 2.75 56.5 straight pipe, 9.5 mm ID, 3.5 m long 

B-56 4.12 60.2 straight pipe, 9.5 mm ID, 3.5 m long 

B-58 4.12 49.2 3.5 m, 90-deg bend, 9.5 mm ID 

B-59 4.12 49.2 4.5 m, 90-deg bend, 9.5 mm ID 

B-60 4.12 49.0 straight pipe with contraction 

B-61 4.12 49.0 straight pipe with expansion 

B-66 2.75 54.1 straight pipe, 15.9 mm ID, 3.5 m long 

B-68 4.12 49.4 straight pipe, 15.9 mm ID, 3.5 m long 

B-71 4.12 49.4 3.5m, 90-deg bend, 15.9 mm ID 

Of the above tests, four tests were simulated with FSP. These are Tests B-59, B-68, B-61 and 

B-60. Together they represent a variety of piping configurations. The results of the 

simulations are presented below. 

4.6.1.2 FSP Simulations of NIST Experiments 

calculate release rates of dissolved gas. Hence the model was implemented with two usersupplied parameters. These are the critical gas bubble radius, rc, and the gas 

A FSP input model for each discharge piping configuration was developed. The 

following geometrical information can be obtained from Gann et al, 1995: 

Storage Vessel Volume: 4.06 x 10-3 m3

Storage Vessel ID: 100 mm

Discharge Valve: 31.8 mm inlet to 44.5 mm exit 

Discharge Valve Stroke: order of 10 ms

Nitrogen Make-up Tank: 88 x 10-3 m3

Nitrogen Charge Tubing: 19.1 mm ID 

Nitrogen Solenoid Valve: 30 kPa pressure difference opening, stroke 50 ms

Discharge Piping Sizes: 9.5 mm and 15.9 mm ID 

Discharge Piping Lengths: 3.5 and 4.5 m 

Recovery Tank Volume: 15 x 10-3 m3

The above information was used to developed the input models for FSP for each of 

the tests. The models are not described in detail or presented here. They were very similar,

conceptually, to the sample problem described in detail in Section 5 of this report. The

supply vessel (called storage tank in these tests) and discharge piping were each modeled 

with PIPE components, with a VALVE component in between them, simulating the initiation 

136
of the discharge process. In modeling the discharge piping, junctions were located at the 

physical locations of area changes (contraction or expansion) and bends. Appropriate loss 

coefficients were entered to simulate the energy losses at these locations. Figure 4.6.3 shows 

the general FSP model used for these simulations. 

These cases were all run with the value of rc in the range of 7.5e-9 to 1.0e-8. The

discharge coefficients at the system exit were chosen as 1.0, 0.8 and 0.8 for the subcooled,

two-phase, and single phase vapor conditions, respectively. The value of Γ0 was the same for 

all the tests, at 10,000. Figures 4.6.4 through 4.6.7 respectively show the calculated and 

measured pressure responses in the supply vessel, for Test B-59, B-68, B-61 and B-60. Not

shown here are the pressure histories at various locations along the discharge piping. 

Generally, at these locations, the measured pressure would increase sharply as the discharge 

was initiated, and then decrease gradually as the fluid left the system. This trend was picked 

up by the code. Overall, the code picks up the qualitative trends in the pressure history. If the

time for complete discharge of the agent is taken as the time at which the pressure drops to its

lowest value, the code seems to calculate this generally within about 25%. A better 

agreement between calculation and data could probably be obtained by varying rc and Γ0. 

This has not been attempted here, primarily because there are inherent uncertainties in these 

parameters, and a code user would normally vary these parameters in design calculations to

get conservative results for design purposes. Generally, a higher value of Γ0 would release 

more dissolved gas earlier in the transient. This would generally delay the complete 

discharge of the bottle slightly. Users can utilize this information to develop conservative

designs in applying the FSP code. 

137
Figure 4.6.1: NIST Test Apparatus (adapted from Gann et al, 1995)

138
FLOW CONFIGURATION

Receiver (a) 

3.5 m, straight

(b) 

3.5 m, with 90o

 bend 

Receiver 

(c) 

4.5 m, with 90o

 bend 

Receiver 

Receiver (d) 

Expansion, 9.5x15.9 mm

Receiver (e) 

Contraction, 15.9x9.5 mm

Receiver (f) 

Tee 

Figure 4.6.2: Piping Configurations for NIST Experiments 

139
sngljun 301

snglvol

205 

Pipe 300

Valve 

Pipe 

203 

tmdpvol

500 

Valve 405 

Nitrogen Make

up Tank 

Test Section

Val Discharge Pipe ve 203 

Recovery Tank

Figure 4.6.3: FSP Input Model

140
Figure 4.6.4: Comparison of Calculated and Measured Bottle Pressure for B-59 

141
Figure 4.6.5: Comparison of Calculated and Measured Bottle Pressure for B-68 

142
Figure 4.6.6: Comparison of Calculated and Measured Bottle Pressure for B-61 

143
Figure 4.6.7: Comparison of Calculated and Measured Bottle Pressure for B-60 

 

144
5. Comparison of Code Predictions to Lehigh Experimental Data 

5.1 Comparison with Experimental Data Collected at Lehigh Test Facility 

In this section a detailed comparison of the computer code predictions and the 

experimental values for a rapid transient discharge of a suppressant flow will be examined. 

To fully test the code, all aspects of an actual suppressant discharge must be considered. For 

this reason we feel the Lehigh experimental data should be the benchmark for comparison 

since it includes the most extensive experimental data to date. This is due to the fact that we 

were able to obtain data for local temperatures, void fraction, and mass flow rate which were 

previously unavailable. In addition, we have extensive data for a number of different flow 

conditions and piping configurations. 

When using the computer code to make predictions, we first had to determine a few 

operating parameters to represent the operating conditions of a specific run. A detailed list of

these parameters are given in the Code Description Section of this report. Some of these 

parameters are straight-forward and pertinent to the piping configuration, such as diameters, 

lengths, etc., and operating conditions such as initial pressure, temperature, etc. However, 

some of these parameters are part of the modeling used in the code. One of these models is 

the gas dissolution model, which is discussed in the Code Description Section. Dependent on 

the fluid used during the discharge, the gas dissolution model requires the following 

constants: 

1 – The critical bubble radius when nitrogen is released from the fluid, 

2 – The release rate of nitrogen. 

A third parameter, that needs to be defined as input to the code, is the discharge 

coefficient for choked flow conditions. When using the computer code to make predictions, 

we first had to determine a few parameters. First off, there is the critical radius, which can be 

found for a specific run by estimating when the release of nitrogen begins. This is the point 

near the beginning of the run when the pressure stops dropping rapidly, represented as a 

change in slope in the pressure histogram. With this pressure and the surface tension of the 

fluid, a critical radius can be calculated as described in the Section 4.2.4.6. After examining 

the first set of experimental runs (#1-#7), it was determined that the value was in the range of

5-20 nanometers. As for the second parameter, the rate of degassing, it was found that the 

range was anywhere from 1000-10,000 depending on the system. Last, the discharge 

coefficient must be chosen. Since the default is 1.0, we knew that the value must be

relatively close to this value. For ease in use and for consistency, it is desirable to choose 

one value for each of these parameters for use when running the predictive computer code.

To determine the appropriate values for these parameters we had to examine the sensitivity of

changing their values and the effect they made on the code’s output and ultimately, compare

them to our experimental data. For the comparison, experimental run #4 was chosen as the

test sample. 

The first step in our sensitivity study was to examine the effect on the computer code 

output when changing the critical radius. To isolate the effect of changing this parameter, we 

kept the values for the gas release rate and the discharge coefficients constant at 1000 and 

145
0.8, respectively. Our findings showed that a lower value for the critical radius caused the 

slope change in the source vessel pressure trace to occur later. Therefore, this parameter is 

important for determining when the gas begins to release. By comparing the code outputs for 

values of 7.5 & 10 nanometers to the experimental results, as seen in Figure 5.1 we chose a 

final value of 10 nanometer for our critical radius. The resultant trace is the red line in the 

figure, with the blue line representing the 7.5 nanometer output and the black line showing 

the experimental data. 

The second step in this study was to specify a value for the rate of nitrogen gas 

release. Once again we kept the parameters that were not being studied constant, the 

discharge coefficient was 0.8 and the critical radius was 7.5 nanometers. This parameter 

influences the speed at which the nitrogen gas is released once the process begins, therefore 

changing the shape of the pressure histogram. This can be seen in Figure 5.2, where the 

black line is the experimental data, the blue line is a value of 1000, the red line is a value of 

2500, and the green line is a value of 5000. It is obvious that at lower values (blue) the 

pressure trend remains straight and uniform, but at higher values (green) the speed of release 

of the nitrogen causes a pressure recovery or “bump” in the trace. From this comparison 

alone, it is hard to determine the best value to use for this parameter, so we made a similar 

comparison using void fraction. This is shown in Figure 5.3, with an increasing average void

fraction as the value for the gas release rate parameter increases. What this figure shows us,

that the other did not, is that a value just below the trace for the value of 2500 would seem to 

work the best. With observations from the above two figures, we chose a gas release value

of 2000. 

The last step in this procedure was to pick a value for the discharge coefficient. This 

parameter influences when and where choking will occur in the system. We tried values

ranging from 0.75-1.15 and found that there was only a negligible difference in our output. 

Therefore we decided to go with a value of 0.9 for the discharge coefficient. Now that all of 

our operating parameters are chosen, a final comparison with the experimental data of Run 

#4 was conducted. The result shows good agreement and is shown in figures 5.4 & 5.5, with 

the red lines being the prediction and the black lines being the experimental data. While we 

have only compared the code output to one experimental run, we feel that it was necessary to 

chose only one value for each parameter so that in future comparisons we do not have to take 

into account any effect changing these values might have on the output. Also, this run was 

chosen over the others because it is one of the repeatability tests and is in the range of typical 

suppressant discharge applications. 

In this next section, the focus will be on comparing the computer code prediction with 

an entire experimental run. The experimental run used in this comparison is the same run as

used earlier in Section 3, run #5. The first quantity we will look at is the cumulative mass in 

the source vessel as time progresses, shown in figure 5.6. The prediction shows excellent 

agreement with the experiment. The only noticeable differences being the final mass value 

and the speed of the discharge represented by the different slope gradients. Even with these 

differences, the prediction gives an end of liquid flow at approximately 1800 milliseconds

compared to a measured time of 1600 milliseconds in the experiment, a very good estimate 

of discharge time. 

146
The next comparison examined in this study was on the void fraction experienced at 

the end of the piping system before the discharge reached the collection vessel. Again, the 

agreement between the prediction and experiment was good as figure 5.7 shows. While the 

experiment is influenced by the local flow regime giving the output a fluctuating appearance,

the prediction has a smooth trace due to the fact that it calculates the void fraction on mass 

balance terms giving an average void fraction value for how much liquid and vapor are 

present in the system. A similar trend can be seen in both, a sudden drop as the pipe is filled 

with liquid and then an increase as nitrogen comes out of solution and raises the average 

amount of gas in the flow until the run ends and there is no more liquid flowing. The only

glaring difference is in the final value reached at the end of the run. The prediction never

goes back to an all-vapor reading of 1.0, informing us that in the simulation there is still 

liquid present in the pipe. In a realistic situation, the liquid would eventually drain into the 

collection vessel since it is the lowest part of the system, but the prediction is finding this

value through a mathematical simulation so it may not account for this. Neglecting this 

problem, a good representation of the void fraction can be found using the prediction. 

The most important comparison for this study would have to be a look at the pressure 

traces, since this is what drives the discharge and delivers the suppressant to the fire. Figures 

5.8-5.12 show the prediction compared to the experiment for the pressures in the source 

vessel, at pipe position #1, at pipe position #2, at pipe position #3, and in the collection 

vessel. They all agree very well with the experimental data with the only discrepancy being 

the final pressure at the end of the run. While the shape of the curves varies slightly, they all 

show the same trend and a similar slope as the pressure decreases during the discharge. This 

observation is most exciting since most of the experimental data available is in the form of

pressure histograms, allowing for the code to be more widely used in comparisons and in 

simulating the output of various system configurations. 

 

While these comparisons were available with all of the runs conducted in this study, 

there were also special measurements made in some cases. For instance, we took

experimental measurements of the pressure drop across different pipe fittings. While the

prediction used a literature value for the friction loss factor associated with the fittings, we 

can still compare what the pressure drop would be over these fittings. This gives us an idea 

of whether the literature values are accurate for a two-phase flow through different pipe 

configurations. Figures 5.13-5.15 show the experimental data for run #9 compared to the 

predicted data for the pressure drop across a capped tee, 90° elbow and connecting union, 

respectively. This run used the more simple of the two alternative piping configurations 

shown previously in Figure 3.11 of section 3. They all show good agreement, with the 

prediction claiming that the pressure started to fade earlier than the experiment and gradually 

fell to zero where the experiment jumped down when the liquid ran out and the vapor/gas 

flow started. We do not see the sudden drop in pressure at the end of the liquid flow that is 

evident in the experiment. Although, the prediction shows an end of the run in the same

general vicinity as the drop in the experimental trace, proving that this is most likely the end 

of the liquid flow. Once again, the code predicts the physical results reasonably well even 

with complicated piping layouts. 

The last comparison we examined dealt with the temperature measurements.

Unfortunately, the predictive code does not have the fluid to pipe heat transfer turned on, 

simply meaning that it does not take into account the fact that there is heat transfer with the 

147
metal walls surrounding the flow. This would cause the flow to either heat up or cool down 

depending on the ambient conditions in the room. It would also be influenced by the heat 

capacity of the metal walls coming into contact with the flow (i.e. vessels and pipes). 

Without this included, the only real temperature change will be from the any phase change in

the system or any compression or expansion of the contents of the system. In specific, the 

degassing of the nitrogen and the adiabatic expansion of the nitrogen near the end of the run 

would be the main contributors. But since the code is not influenced by the ambient

temperature of the system, the code predicts a much lower temperature than was found in the 

experimental runs. The comparisons in this analysis use experimental run #16, which used 

HFC-125 as a test fluid in which the largest temperature change experience in the 

experiments. In Figure 5.16, a comparison is made between the experimental measurement of 

the source vessel temperature with the shielded thermocouple which was at the bottom of the 

vessel and immersed in the fluid and the fluid temperature predicted by the code. The trend is 

similar except for the magnitude and speed of the temperature drop, but the fact that the 

shapes are alike is promising that at least the phenomena behind the temperature drop is 

correct. Figure 5.17 shows a temperature comparison for the pipe position #1 in the piping 

system. Once again, the fluid temperatures are being shown and a similar result is obtained. 

Although, in this case, the differences are more pronounced than in the source vessel 

comparison. It is clear that the point when the temperatures began to drop is very close, but 

the magnitude of the drop is very different. As we expected, the heat transfer coefficient in

the pipe is greater than anywhere else in the system and therefore the short-comings of not 

including the fluid to pipe heat transfer in the code are most obvious in this region. The

comparison of the gas temperature in the collection vessel is shown in Figure 5.18. The 

collection vessel shows a slightly different result, a temperature drop at the beginning of the 

run was experienced due to the fact that the liquid was evaporating to raise the pressure

enough to reach its saturation equilibrium. This trend stopped when the pressure rose high

enough that there was no longer any need for the liquid to evaporate. Although the magnitude 

of the drop is different, both the prediction and the experiment show the same insight into 

whether the measurements were in the right magnitude and trend. 

After study of all of the comparisons between the predictions and the experimental 

runs the general agreement is very good. In most cases, the prediction gives a good 

estimation of the physical characteristics of a two-phase flow in a model system. The code 

can even handle two different test fluids with reasonable accuracy. It also predicted three

different piping configurations well, showing its versatility in modeling any potential system.

While these points support the code’s advantage for use with suppressant systems, it must 

also be noted that there were a few limitations of the code. The most notable being that it 

will not always work for any initial conditions inputted for a specific test fluid. While this is

most likely due to a limitation with the thermodynamics package, it still leaves a small

amount of uncertainty when running the code. The user must slightly deviate from the actual 

run condition in order to fall in a region of the thermodynamics package that will yield an 

output. For example, in experimental run #5 the initial conditions were entered into the code 

through the input deck, but an error occurred before the program could complete the 

simulation. But, it turns out by changing the ambient temperature from 31o Celsius to 30o 

Celsius the simulation works fine. This is only a minor change, but it still makes operation 

of the simulation slightly cumbersome. 

148
Another area that caused a problem was when we tried to simulate our low pressure

runs (2000-2500 kPa initial source vessel pressure). There were two problems, which arose

during these tests. The first was similar to the problem stated previously, something with the 

thermodynamics package. These experimental runs (#1 & #8), were conducted with low 

downstream pressures so that flashing could occur in the pipe directly after the start of the

run and the actuation of the valve. The code would not allow these conditions to run the 

simulation due to the fact that it believed that there was non-equilibrium in the downstream

areas. The downstream pressure had to be raised to a pressure closer to the saturation

conditions of the test fluid in order for the code to operate. The other problem we found with 

these low pressure test runs was in regard to the high fill level we used in the source vessel.

In experimental runs #1, #8, #11, #14, & #18 a fill level of 80%-90% was used to try and cut 

down on the amount of nitrogen present in the system. The hope was to keep the system

pressure low enough to facilitate flashing of the suppressant during discharge. This proved 

to be a problem with the code’s predictive output. While this initial condition would not 

cause the code to stop running with an error, it did cause an odd result to be predicted. For 

example, in run #1, a total mass of test fluid (HFC-227ea) filled into the source vessel at the 

beginning of the run was about 4.9 kilograms. However the code prediction never even got 

close to that value during its predicted discharge, it reached a final value of about 3.1 

kilograms. The only explanation we can add at this time is the possibility that the code was 

expected to only be predicting discharges that held a fill level of about 60%, instead of the

test fill level of about 85%. But, this problem would cause a significant difference in all 

aspects of the prediction to be encountered. 

The last problem that we encountered was with the high pressure runs (4750-

5000 kPa initial source vessel pressure). Once again, we believe this had something to do 

with the properties package. When the initial conditions were fed into the code through the 

input file an error would occur immediately, causing the program to stop. We found that the 

only way to get these runs to continue was by lowering the source vessel temperature 

considerably, while keeping the downstream temperature around the ambient value. Even 

then we could not always reach the desired operating pressure. For example, in run #6, the 

source vessel temperature was dropped from 31o Celsius to 2o Celsius and we could only get 

an initial source vessel pressure of 4860 kPa in place of 4910 kPa to operate. When these 

runs finally ran, they still only continued for half of the time that the experimental runs were 

conducted for. It seems that the computer code would finish early in the run (i.e. 1500 

milliseconds for run #6) most likely due to the temperature dropping below the range of the 

thermodynamics package. Even with these observed problems, the code did very good job 

when predicting most runs in the application area, an initial source vessel pressure between 

3000-4500 kPa with a fill level of about 65%. This is the first code that has been able to be 

compared to void fraction data and instantaneous mass inventory data and it does an 

excellent job in estimating these values. Plots showing the comparisons between the 

predicted and experimental data can be found in Appendix C for 18 of the test case runs at 

the Lehigh University facility. 

149
0

500

10001500200025003000350040004500

-500 0 500 1000 1500 2000 250

Time [msec]

Pressure [kPa]

r = 7.5 n

r = 10 nm 

experime

Gas Release Rate = 

Discharge Coeffic

Figure 5.1: Comparison of Experimental and Predictive Source Vessel Pressure Data in Sensitivity Study on 

for Conditions of Run #4. 

150 
0

500

10001500200025003000350040004500

-500 0 500 1000 1500 2000 2

Time [msec]

Pressure [kPa]

Gamma = 1000 Gamma = 2500 Gamma = 5000 Experimental

Critical Ra

Discharge C

Figure 5.2: Comparison of Experimental and Predictive Source Vessel Pressure Data in Sensitivity Study o

Value for Conditions of Run #4. 

151
0.0

0.20.40.60.81.01.2

-500 0 500 1000 1500 2000 2

Time [msec]

Void Fraction

Gamma = 1000 Gamma = 2500 Gamma = 5000 Experimental

Critical Ra

Discharge C

Figure 5.3: Comparison of Experimental and Predictive Void Fraction Data in Sensitivity Study on Gas Re

Conditions of Run #4. 

152
0

500

10001500200025003000350040004500

-500 0 500 1000 1500 2000 25

Time [msec]

Pressure [kPa]

Prediction

Experiment

Figure 5.4: Comparison of Experimental and Predictive Pressure Data using Final Values for Code Opera

Conditions of Run #4. 

153
0.0

0.20.40.60.81.01.2

-500 0 500 1000 1500 2000 25

Time [msec]

Void Fraction

Prediction

Experiment

Figure 5.5: Comparison of Experimental and Predictive Void Fraction Data using Final Values for Code Ope

Conditions of Run #4. 

154
0.0

0.51.01.52.02.53.03.5

-500 0 500 1000 1500 2000 25

Time [msec]

Cumulative Mass [kg]

Experiment

Prediction

Figure 5.6: Comparison of Cumulative Mass for Experimental Run #5 and Predictive Computer C

155
0.0

0.20.40.60.81.01.2

-500 0 500 1000 1500 2000 25

Time [msec]

Void Fraction

Prediction

Experiment

Figure 5.7: Comparison of Void Fraction for Experimental Run #5 and Predictive Computer Co

156
0

500

10001500200025003000350040004500

-500 0 500 1000 1500 2000 2

Time [msec]

Pressure [kPa]

Prediction

Experiment

Figure 5.8: Comparison of Source Vessel Pressure for Experimental Run #5 and Predictive Compute

157
0

500

10001500200025003000350040004500

-500 0 500 1000 1500 2000 25

Time [msec]

Pressure [kPa] Prediction

Experiment

Figure 5.9: Comparison of Pipe Position #1 Pressure for Experimental Run #5 and Predictive Compu

158
0

500

10001500200025003000350040004500

-500 0 500 1000 1500 2000 25

Time [msec]

Pressure [kPa]Experiment

Prediction

Figure 5.10: Comparison of Pipe Position #2 Pressure for Experimental Run #5 and Predictive Compu

159
0

500

10001500200025003000350040004500

-500 0 500 1000 1500 2000 25

Time [msec]

Pressure [kPa] Experiment

Prediction

Figure 5.11: Comparison of Pipe Position #3 Pressure for Experimental Run #5 and Predictive Compu

160
0

500

10001500200025003000350040004500

-500 0 500 1000 1500 2000 25

Time [msec]

Pressure [kPa]

Prediction

Experiment

Figure 5.12: Comparison of Collection Vessel Pressure for Experimental Run #5 and Predictive Comp

161
0

50

100150200250300350

-500 0 500 1000 1500 2000 2500 

Time [msec]

Pressure [kPa]

Experiment

Prediction

Figure 5.13: Comparison of Pressure Drop Across a Capped Tee for Experimental Run #9 and Predictive C

162
0

50

100150200250300350

-500 0 500 1000 1500 2000 2500 

Time [msec]

Pressure [kPa]

Experiment

Prediction

Figure 5.14: Comparison of Pressure Drop Across a 90° Elbow for Experimental Run #9 and Predictive Co

163
0

20406080

100120140160180200

-500 0 500 1000 1500 2000 2500 

Time [msec]

Pressure [kPa]

Experiment

Prediction

Figure 5.15: Comparison of Pressure Drop Across a Union for Experimental Run #9 and Predictive Com

164
-15

-10

-5

05

1015202530

-500 0 500 1000 1500 2000 25

Time [msec]

Temperature [C]

Experiment

Prediction

Figure 5.16: Comparison of Source Vessel Fluid Temperature for Experimental Run #16 and Predictive Co

165
-15

-10

-5

05

1015202530

-500 0 500 1000 1500 2000 2

Time [msec]

Temperature [C]

Experiment

Prediction

Figure 5.17: Comparison of Pipe Position #1 Fluid Temperature for Experimental Run #16 and Predic

166
-20

-15-10

-5

05

1015202530

-500 0 500 1000 1500 2000 25

Time [msec]

Temperature [C]

Prediction

Experiment

Figure 5.18: Comparison of Collection Vessel Gas Temperature for Experimental Run #16 and Predictive C

167
6. SUMMARY 

This code development effort was aimed at developing a flexible computer program

that would allow prediction of the performance of fire suppressant fluid delivery systems,

including a variety of fluids and delivery system piping configurations. The approach taken 

was to start from an existing transient two-phase flow analysis code, available in the public 

domain, and modifying it appropriately. In this case, the starting code was RELAP5/MOD3.2 

(called the base code in this report), used for analysis of accidents and transients in the 

nuclear power industry (RELAP5, 1995). It represents the current state-of-the-art for 

modeling two-phase flows in piping systems, and was written for water-steam systems. It 

allows flexible modeling of piping systems, with the user being able to control the size of

control volumes in the numerical solution procedure, and also control the computational time

step. 

The major modifications made to the base code consisted of the following: 

1) Thermodynamic and transport properties were incorporated for Halon-1301, CO2, HFC227ea and HFC-125. The intent was to be able to benchmark the code against current 

systems and methods that operate with Halon-1301, and be able to use the code with 

other fluids being considered as replacements for Halon-1301. The thermodynamic

properties were implemented by using the REFPROP property package, developed by 

NIST (Gallagher et al, 1993), and incorporating it into the base code. The method by 

which this was carried out is flexible, and would facilitate incorporation of properties for

additional fluids, such as CF3I, if needed. The modifications were implemented such that 

the user can select one of the available fluids, as an input option. 

2) Constitutive relationships were modified as necessary to incorporate fluid-independent 

models. A review of the base code showed that many of the constitutive relationships in 

the base code could be directly adapted for fluids other than water. The areas where 

modifications were made were in the models for the calculation of interphase drag, and 

wall heat transfer. 

3) A mechanistic model was implemented for simulating the release of dissolved gas. This 

required incorporating source and sink terms in appropriate mass conservation equations, 

and developing a constitutive model for calculating the rate of gas release. There are 

inherent uncertainties in being able to calculate release rates of dissolved gas. Hence the 

model was implemented with two user-supplied parameters. These are the critical gas

bubble radius, rc, and the gas release coefficient, Γ0. The model calculates the release of 

dissolved gas as a function of local conditions. This allows prediction of gas release in a 

continuous fashion as a function of spatial location, as the suppressant fluid is discharged 

from a supply vessel and flows through a piping system. 

The resulting code, FSP, Version 1.0, is a one-dimensional, two-fluid model of two-phase 

flow, allowing the calculation of nonhomogeneous and nonequilibrium conditions between 

the phases. It retains all the flexibility of the base code in being to model, via user input, a 

variety of piping systems. It allows the user to select, via input, one of several fluids being 

considered currently as viable alternatives to Halon-1301. It was developed from a public 

domain code, and so was developed as a non-proprietary code. These were the characteristics 

168 
required in this code development project. The code has been installed successfully on 

personal computers with the WINDOWS 95/98 operating system. Some simple procedures 

were developed for code installation and execution and are shown in Section 4.4 of this 

report. 

An experimental program was also carried out in order to improve the understanding 

of the discharge process, to determine some related parameters, and to assess the code against 

the experimental data. The discharge test loop dimensions were selected to be similar to the 

actual systems used in applications. Special hardware and instrumentation features were 

included in the loop to allow the study of different operating conditions and to be able to 

measure parameters that have never been measured before. In the present test loop it was

possible to experimentally measure the dynamic discharge flow rate of the suppressant,

dynamic void fraction and fluid temperature in the discharge piping. Additional 

measurements were carried out to determine pressure losses across pipe fittings during this

transient two-phase discharge. To the best of our knowledge these are the first successful 

measurements of above parameters in a fire suppressant discharge loop. Experiments used 

water, CO2, HFC-227ea, and HFC-125 as discharge fluids. All of the suppressants, except 

CO2, were driven using compressed nitrogen at pressures up to 50 bars. 

The total duration of the discharge varied between 1.5 to 6.0 seconds depending on 

initial pressure and the piping configuration. Experimental tests showed the following

characteristics of discharge process:

- During the discharge, first the liquid suppressant is discharged, followed by vapor 

suppressant mixed with nitrogen gas, 

- Nitrogen dissolved in the suppressant comes out of the liquid at small range of 

pressures, and it takes some time to have all nitrogen come out of the liquid 

suppressant, 

- Void fraction at the end of the piping was found to stay around 0.5, 

- The discharge flow rate was found to be highest right after the discharge, and would 

decrease gradually toward zero at the end of the discharge. 

- When using HFC-227ea, the fluid temperatures in the discharge loop were found to 

be cooling only a few degrees Celsius, 

- When using HFC-125, the fluid temperatures in the discharge loop were found to 

drop from room temperature down to almost zero degrees Celsius. 

The code has been assessed against the above-described experiments on suppressant 

agent discharge, conducted with several fluids, in various piping configurations. The results 

show that the code is capable of calculating the pressure history and discharge flow rates 

similar to the measured data. The code is also able to pick up the qualitative behavior 

associated with the release of dissolved gas as the delivery system undergoes rapid 

depressurization during the discharge of the suppressant agent. Because of these 

characteristics, the code appears to be a viable tool for analysis or design of suppressant 

delivery systems. 

Because the code is very flexible in being able to model different piping systems via user 

input, it offers the user many options in developing an appropriate input model to describe 

the physical system. The input model essentially provides the initial and boundary conditions

169
for the equations that the code solves. This can become complicated. A detailed input manual

is appended to this report that describes how to build an input model. However, for 

convenience, a sample problem has been described in Section 4.5 of this report. It is 

recommended that this sample problem be used, along with the input manual, to develop an 

understanding of how the user can efficiently develop an input model for a desired 

application and carry out the code calculation. 

Wall heat transfer models are present in the code, and the user can use them to model the 

metal walls of the supply vessel and piping as passive heat sources. However, typical 

suppressant discharge calculations examine transient times of the order of a few seconds, 

which are small compared to the time needed for stored heat in the walls to have an impact 

on the transient. Hence, the sample problem shown in Section 4.5 and all the code 

assessments carried out so far have not modeled wall heat transfer. And, in order to keep the

input manual reasonably simple to use, the instructions of how to model wall heat transfer 

have not been included in the input manual. If the user encounters an application where this 

may be considered important, the documentation of the base code (RELAP5, 1995) provides 

this information. The user is cautioned that while the expectation is that these models will 

function as designed, they have not been exercised or assessed during the code development 

effort. 

One code limitation that has been encountered is that the code sometimes appears to have 

difficulty beginning the computation. A common cause that was found in these cases is the 

initial thermodynamic state of the fluid specified by the user. The sample problem described 

in Section 4.5 contains a discussion on this topic and the user is urged to review this 

discussion. 

Once the user builds an input model to describe a delivery system, there are essentially 

three input parameters that the user can vary to examine analytical uncertainties: a) The 

discharge coefficient directly affects the magnitude of the discharge flow rate. Note that the 

code uses this input value as a multiplier on the junction flow area. There are three separate 

discharge coefficients, for subcooled liquid, two-phase conditions, and single-phase vapor. 

Based on the assessments carried out, it appears that these should have values of 1.0, 0.8-0.9 

and 0.8-0.9 respectively. However, it is recognized that the discharge coefficient depends on 

the geometry of the exit. A smooth converging-diverging nozzle would have a different 

discharge coefficient than a sharp-edged orifice plate. The user should be aware of this when

specifying the discharge coefficient in the input model. Essentially, the discharge coefficient 

allows the user to account for uncertainties in the critical flow model of the code and the 

geometry of the junction. b) The critical gas bubble radius, rc controls the pressure at which 

gas release is initiated. The assessments carried out indicate that this parameter should be in 

the range of 7.5e-9 to 1.0e-8 m. The user should examine the assessment cases presented to 

evaluate whether this is reasonable or whether a value outside this range should be used for

this parameter. c) The gas release coefficient, Γ0, is directly proportional to the gas release 

rate. The assessments carried out indicate that this should be in the range of 103

 to 104

Kg/m3

-s. However, it is recognized that there is significant uncertainty in this parameter. 

The user should plan to vary this parameter for a specific application and gain understanding 

of how this would affect the overall time of suppressant discharge. 

170
Finally, the description of the code theory and structure, in Sections 4.2 and 4.3 of this 

report, are somewhat brief. This was intentional. The emphasis in these sections is to provide

the user with some basic information, and some additional guidelines to successfully build

the input model and run the code. A very detailed description of the code theory, code 

structure, and numerical solution method can be found in the documentation of the base code 

(RELAP5, 1995). An attempt has been made in this report to identify all the changes made to 

the base code, and provide sufficient description of the changes. 

171
7. REFERENCES 

Chen, J.C., Tuzla, K., and Sundaram, R., Development of Computer Code for Analysis of 

Fire Suppressant Fluid System, Technical Proposal submitted by Lehigh University, BAA 

Number 1773-6B, 1997. 

Chisholm, D., “A Theoretical Basis for the Lockhart-Martinelli Correlation for Two-Phase 

Flow”, Journal of Heat and Mass Transfer, 10, 1967, pp 1767-1778. 

Cokmez-Tuzla, A.F., Tuzla, K., and Chen, J.C., "High Speed Thermometry for Detection of 

Transient Liquid Contact on Superheated Surfaces", Temperature Its Measurement and 

Control in Science and Industry, 6 (ed. J. Schooley), American Institute of Physics, pp. 1173-

1177, 1993. 

Crane Handbook, “Flow of Fluids through Valves, Fittings, and Pipe”, Technical Paper No. 

410, Published by Crane Company, PA, 24th printing, 1988. 

DiNenno, Philip J. and Budnick, Edward K. Jr., “Halon 1301 Discharge Testing: A Technical 

Analysis”, National Fire Protection Research Foundation, Quincy, MA, 1988. 

Edwards, A.R. and O'Brien, T.P., "Studies of Phenomena Connected with the 

Depressurization of Water Reactors", Journal of British Nuclear Energy Society, Vol. 9, pp. 

125-135, 1979. 

Elliott, D. G., Garrison, P. W., Klein, G. A., Moran, K. M., and Zydowicz, M. P., “Flow of 

Nitrogen-Pressurized Halon 1301 in Fire Extinguishing Systems”, Jet Propulsion Laboratory 

Report 84-62, November 1, 1984. 

EPRI Report NP-2370, “The Marviken Full-Scale Critical-Flow Tests, Volume 1: Summary 

Report”, Electric Power Research Institute, Palo Alto, CA, 1982. 

Gallagher, J., McLinden, M., Morrison, G., and Huber, M., “NIST Thermodynamic 

Properties of Refrigerants and Refrigerant Mixtures Database (REFPROP)”, Version 4.0, 

November 1993. 

Gann, Richard G., Editor, NIST Special Publication SP 890, Volume 1, Section 8.6, 

September, 1995. 

GOTHIC Containment Analysis Package Version 5.0, Numerical Applications, Inc., User 

Manual (NAI 8907-02 Rev. 6), and Technical Manual (NAI 8907-06 Rev. 5), 1995. 

Grosshandler, William L., Gann, Richard G. and Pitts, William M., Editors, “Evaluation of

Alternative In-Flight Fire Suppressants for Full-Scale Testing in Simulated Aircraft Engine 

Nacelles and Dry Bays”, NIST Special Publication SP 861, April 1994. 

Hirt, C. W., Oliphant, T. A., Rivard, W. C., Romero, N. C., and Torray, M. D., “SOLALOOP: A Non-Equilibrium, Drift Flux Code for Two-Phase Flow in Networks”, Report No. 

L-7659, Los Alamos National Laboratory, Los Alamos, NM, 1979. 

172
Ishii, M., “One-Dimensional Drift-Flux Model and Constitutive Equations for 

Relative Motion between Phases in Various Two-Phase Flow Regimes”, Argonne National 

Laboratory Report ANL-77-47, 1977. 

Ishii, M., and Mishima, K., “Study of Two-Phase Fluid Model and Interfacial Area”, 

NUREG/CR-1873, ANL-80-111, 1980. 

Marviken Power Station, "The Marviken Full Scale Critical Flow Tests: Volume 1: 

Summary Report,", NUREG/CR-2671, MXC-301, United States Nuclear Regulatory 

Commission, Washington, D.C., also EPRI NP-2370, 1982 

Meyer, C. A., McClintock, R. B., Silvestri, G. J., and Spencer, R. C. Jr., ASME Steam

Tables, 1967. 

Mishima, K., and Ishii, M., “Flow Regime Transition Criteria Consistent with Two-Fluid 

Model for Vertical Two-Phase Flow”, NUREG/CR-3338, April 1984. 

Pitts, W. M., Yang, J. C., Gmurczyk, G. W., Cooper, L. Y., Grosshandler, W. L., and 

Cleveland, W. G., “Fluid Dynamics of Agent Discharge”, Chapter 3 of NIST Special 

Publication 861, Editors: William L. Groshandler, Richard G. Gann and William M. Pitts, 

Gaithersburg, MD, April 1994. 

RETRAN-03: A Program for Transient Thermal-Hydraulic Analysis of Complex Fluid Flow 

Systems, EPRI NP-7450, Volume 1, May 1992. 

TRAC-PF1/MOD2 Code Manual, Los Alamos National Laboratory, NUREG/CR-5673 or 

LA-1203-M, Volume 1-3, 1992. 

TRAC-BF1/MOD1 Code Manual: An Advanced Best-Estimate Computer Program for BWR 

Accident Analysis, Idaho National Engineering Laboratory, NUREG/CR-4356 or EGG2626, Volumes 1-2, 1992. 

TRAC-BF1/MOD1 Models and Correlations, INEL, NUREG/CR-4391 or EGG-2680, 1992. 

Trapp, J. A., and Ransom, V. H., “A Choked Flow Calculation Criterion for 

Nonhomogeneous, Nonequilibrium Two-Phase Flows”, International Journal of Multiphase 

Flow, 8, 6, 1982, pp 669-681. 

RELAP5 Development Team, “RELAP5/MOD3.2 Code Manuals, Volume 1: Code 

Structure, Models and Numerical Schemes”, U. S. Nuclear Regulatory Commission Report 

NUREG/CR-5535, 1995. 

RELAP5 Development Team, “RELAP5/MOD3.2 Code Manuals, Volume 2: User Manual”, 

U. S. Nuclear Regulatory Commission Report NUREG/CR-5535, 1995. 

173
Taitel, Y.A., Bornea, D., and Dukler, A.E., “Modeling Flow Pattern Transitions for Steady 

Upward Gas-Liquid Flow in Vertical Tubes”, AIChE Journal, Vol. 26, No. 3, May 1980, pp 

345-354. 

Tuzla, K., Palmer, T., Chen, J. C., Sundaram. R. K., and Yeung, W. S., “Development of 

Computer Program for Fire-Suppressant Fluid Flow”, Annual Technical Report to NIST, 

Lehigh University Institute of Thermo-Fluid Engineering and Science Report ITF 98-2, 

Lehigh University, Bethlehem, PA, December 1998. 

Van Wylen, G.J., and Sonntag, R.E., Classical Thermodynamics, 3rd Edition, John Wiley 

and Sons, 1985. 

Wallis, G.B., One-Dimensional Two-Phase Flow, McGraw Hill, New York, 1972. 

Welty, J. R., Wicks, C. E., and Wilson, R. E., Fundamentals of Momentum, Heat, and Mass 

Transfer, 3rd Edition, John Wiley & Sons, New York, pp. 210-211, 1984. 

Yang, J. C., Cleary, T. G., Vazquez, I., Boyer, C. I., King, M. D., Breul, B. D., Womeldorf, 

C., Grosshandler, W. L., Huber, M. L., Weber, L., and Gmurczyk, G. W., “Optimization of 

System Discharge”, NIST Special Publication 890, Gaithesburg, MD, pp 407-645, November 

1995. 

Yang, J. C., Vasquez, I., Boyer, C. I., Huber, M. L. and Weber, L., “Measured and Predicted 

Thermodynamic Properties of Selected Halon Alternative/Nitrogen Mixtures”, Int. J. 

Refrigeration., Vol. 20, No. 2, pp 96-105, 1997. 

Zuber, N., Tribus, M., and Westwater, J. W., “The Hydrodynamic Crisis in Pool Boiling of 

Saturated and Subcooled Liquids”, Second International Heat Transfer Conference, Paper 27, 

Denver, 1961.

174