National Aeronautics and NASA/TM—2017–218239

Space Administration
IS02
George C. Marshall Space Flight Center
Huntsville, Alabama 35812

Tank System Integrated Model: A Cryogenic

Tank Performance Prediction Program
L.G. Bolshinskiy
Jacobs ESSSA Group, Huntsville, Alabama

A. Hedayat, L.J. Hastings (retired), S.G. Sutherlin, and A.R. Schnell

Marshall Space Flight Center, Huntsville, Alabama

J.P. Moder
Glenn Research Center, Cleveland, Ohio

April 2017
 The NASA STI Program…in Profile

Since its founding, NASA has been dedicated to the • CONFERENCE PUBLICATION. Collected
advancement of aeronautics and space science. The papers from scientific and technical conferences,
NASA Scientific and Technical Information (STI) symposia, seminars, or other meetings sponsored
Program Office plays a key part in helping NASA or cosponsored by NASA.
maintain this important role.
• SPECIAL PUBLICATION. Scientific, technical,
The NASA STI Program Office is operated by or historical information from NASA programs,
Langley Research Center, the lead center for projects, and mission, often concerned with
NASA’s scientific and technical information. The subjects having substantial public interest.
NASA STI Program Office provides access to
the NASA STI Database, the largest collection of • TECHNICAL TRANSLATION.
aeronautical and space science STI in the world. English-language translations of foreign
The Program Office is also NASA’s institutional scientific and technical material pertinent to
mechanism for disseminating the results of its NASA’s mission.
research and development activities. These results
are published by NASA in the NASA STI Report Specialized services that complement the STI
Series, which includes the following report types: Program Office’s diverse offerings include creating
custom thesauri, building customized databases,
• TECHNICAL PUBLICATION. Reports of organizing and publishing research results…even
completed research or a major significant providing videos.
phase of research that present the results of
NASA programs and include extensive data For more information about the NASA STI Program
or theoretical analysis. Includes compilations Office, see the following:
of significant scientific and technical data
and information deemed to be of continuing • Access the NASA STI program home page at
reference value. NASA’s counterpart of peer- <http://www.sti.nasa.gov>
reviewed formal professional papers but has less
stringent limitations on manuscript length and • E-mail your question via the Internet to
extent of graphic presentations. <help@sti.nasa.gov>

• TECHNICAL MEMORANDUM. Scientific • Phone the NASA STI Help Desk at

and technical findings that are preliminary or of 757 –864–9658
specialized interest, e.g., quick release reports,
working papers, and bibliographies that contain • Write to:
minimal annotation. Does not contain extensive NASA STI Information Desk
analysis. Mail Stop 148
NASA Langley Research Center
• CONTRACTOR REPORT. Scientific and Hampton, VA 23681–2199, USA
technical findings by NASA-sponsored
contractors and grantees.
NASA/TM—2017–218239

Tank System Integrated Model: A Cryogenic

Tank Performance Prediction Program
L.G. Bolshinskiy
Jacobs ESSSA Group, Huntsville, Alabama

A. Hedayat, L.J. Hastings (retired), S.G. Sutherlin, and A.R. Schnell

Marshall Space Flight Center, Huntsville, Alabama

J.P. Moder
Glenn Research Center, Cleveland Ohio

National Aeronautics and

Space Administration

Marshall Space Flight Center • Huntsville, Alabama 35812

April 2017

i
 Acknowledgments

.The authors would like to extend their appreciation to Main Propulsion Branch (ER22) Management, Kathy
Henkel and Jay Russell, for providing opportunity and support to complete this work. Appreciation also goes to the
Jacob Technology ESSSA Group Team Leads, David Sharp, Mike Reynolds, and Sherri Spotswood for their support
during all stages of the Tank System Integrated Model (TankSIM) development and completion of this Technical
Memorandum. Special thanks go to Jonathan Stephens of ER24 for his excellent and very helpful technical discus-
sions, suggestions, and reviews during different stages of TankSIM work.

Available from:

NASA STI Information Desk

Mail Stop 148
NASA Langley Research Center
Hampton, VA 23681–2199, USA
757–864–9658

This report is also available in electronic form at

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

ii
 TABLE OF CONTENTS

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

2. ANALYTICAL MODEL DESCRIPTION ......................................................................... 2

2.1 Basic Model Approaches for a Settled Liquid ................................................................ 2

2.2 Model Approaches for an Unsettled Liquid ................................................................... 5

3. MASS AND ENERGY CONSERVATION EQUATIONS ................................................. 6

3.1 Ullage ............................................................................................................................ 6

3.2 Ullage Wall (Tank Wall Interfaced to the Ullage) .......................................................... 8
3.3 Bulk Liquid ................................................................................................................... 9
3.4 Liquid Wall (Tank Wall Interfaced to the Liquid) .......................................................... 9
3.5 Wall Liquid (Liquid Film on the Ullage Wall) ............................................................... 10

4. HEAT AND MASS TRANSFER IN THE TANK ............................................................. 11

4.1 External Heat Loads ...................................................................................................... 11

4.2 Ullage-Bulk Liquid Interface Heat and Mass Transfer .................................................. 11
4.3 Ullage Wall Heat and Mass Transfer ............................................................................. 13
4.4 Condensation on the Ullage Wall .................................................................................. 14
4.5 Bulk Vapor Condensation .............................................................................................. 17
4.6 Bulk Liquid Convective Heat and Mass Transfer ........................................................... 17
4.7 Bulk Liquid Boiling ....................................................................................................... 19
4.8 Conductive Heat Transfer Through the Wall Between the Ullage and Liquid
Regions .......................................................................................................................... 21
4.9 Fluid Thermodynamic Peoperties .................................................................................. 21

5. TANK PRESSURIZATION AND ULLAGE VENTING PRESSURE CONTROL ......... 22

5.1 Ullage Pressurization ..................................................................................................... 22

5.2 Ullage Venting Methods ................................................................................................ 24

6. TANK THERMODYNAMIC VENTING PRESSURE CONTROL SYSTEMS ............... 28

6.1 Thermodynamic Venting System Pressure Control Logic .............................................. 28

6.2 Spray Bar Thermodynamic Venting System ................................................................... 29
6.3 Ullage Droplets Heat and Mass Transfer ....................................................................... 39

iii
 TABLE OF CONTENTS (Continued)

7. TANK AXIAL JET THERMODYNAMIC VENTING PRESSURE CONTROL

SYSTEMS .......................................................................................................................... 47

7.1 Axial Jet Thermodynamic Venting System ................................................................... 47

7.2 Heat and Mass Transfer at the Ullage-Liquid Interface ............................................... 47
7.3 Helicoidally Coiled Tube-in-Shell Heat Exchanger Transfer ........................................ 50
7.4 Heat Transfer at Low Acceleration Conditions ............................................................ 52

8. COMPUTER MODEL DESCRIPTION .......................................................................... 55

8.1 TankSIM Program Base ............................................................................................... 55

8.2 Program Possibilities .................................................................................................... 55
8.3 Program Error Handling .............................................................................................. 56
8.4 TankSIM Functionality ............................................................................................... 57

9. PROGRAM VALIDATION AND MULTIPHASE MISSION EXAMPLE ..................... 63

9.1 Multipurpose Hydrogen Test Bed Data Validation ...................................................... 63

9.2 Multiphase Mission Example With Ullage Venting ..................................................... 73

10. SUMMARY AND FUTURE PROGRAM IMPROVEMENT ........................................ 75

10.1 Summary.................................................................................................................... 75
10.2 Future TankSIM Improvements................................................................................. 75

APPENDIX A—CONSERVATION EQUATIONS FOR MULTINODE

ANALYSIS OF CRYOGENIC STORAGE TANKS .................................... 76

A.1 Leibnitz’s Theorem ..................................................................................................... 76

A.2 Conservation of Mass ................................................................................................ 76
A.3 Conservation of Thermal Energy ............................................................................... 78

APPENDIX B—HEAT AND MASS TRANSFER INSIDE THE TANK ............................. 82

B.1 Ullage-Bulk Liquid Interface Temperature ................................................................. 82

B.2 Requirements for Using Flat Plate Correlations to the Cylindrical Segment
of the Tank ................................................................................................................. 83
B.3 Free Convection Heat Transfer From the Tank Wall Inside the Tank ......................... 85
B.4 Ullage Wall Condensation Model ............................................................................... 88

iv
 TABLE OF CONTENTS (Continued)

APPENDIX C—TANK THERMODYNAMIC VENTING PRESSURE CONTROL

SYSTEMS ..................................................................................................... 97

C.1 Liquid Film Heat Transfer ............................................................................................ 97

APPENDIX D—TANK AXIAL JET THERMODYNAMIC VENTING PRESSURE

CONTROL SYSTEMS ................................................................................. 99

D.1 Heat and Mass Transfer on the Ullage-Liquid Interface During Axial
Jet Mixing ..................................................................................................................... 99
D.2 Helicoidally Coiled Tube-in-Shell Heat Exchanger Transfer .......................................... 101
D.3 Heat Transfer at Low Acceleration Conditions ............................................................. 103

APPENDIX E—COMPUTER MODEL DESCRIPTION ..................................................... 105

E.1 Input Files Description ................................................................................................. 105

E.2 Variables Description for Input and Output Files ......................................................... 112

REFERENCES ........................................................................................................................ 119

v
 LIST OF FIGURES

1. TankSIM nodes .......................................................................................................... 2

2. Diagram of control volumes ...................................................................................... 3

3. Tank geometry transformation for unsettled liquid (no mixing) ................................. 5

4. Tank fluids transformation for unsettled liquid during mixing ................................... 5

5. Diagram of mass flows ............................................................................................... 6

6. Diagram of heat flows, not associated with mass transfer .......................................... 7

7. Horizontal heated downward-facing plate model ....................................................... 12

8. Liquid-ullage interface in finite volume tank .............................................................. 12

9. Calculation schematic for condensation on the upper dome ....................................... 15

10. Calculation schematic for condensation on the lower dome ....................................... 16

11. Definition sketch for lower dome convective model .................................................... 18

12. Calculation schematic for convection in the liquid on the upper dome wall ................ 18

13. Pressurization devices schematic: (a) Critical Venturi and (b) critical orifice ............. 22

14. Ullage venting pressure control logic .......................................................................... 25

15. Thermodynamic venting system pressure control logic ............................................... 29

16. Spray bar TVS schematic ........................................................................................... 30

17. Spray bar heat exchanger pressure drop calcualtion schematic ................................... 31

18. Spray bar heat exchanger transfer calculation schematic ............................................ 34

19. Spray injection tube schematic ................................................................................... 36

20. Resistance coefficient versus ratio of the injection tube to orifice velocities ................ 38

vi
 LIST OF FIGURES (Continued)

21. Droplet motion in the ullage ...................................................................................... 41

22. Droplet impingement model schematic ...................................................................... 43

23. Droplet behavior after impingement: (a) Spreading and (b) rebounding ................... 45

24. Axial jet TVS schematic ............................................................................................. 47

25. Axial jet TVS submergence cases schematic: (a) High submergence
and (b) low submergence ............................................................................................ 48

26. Helicoidally coiled tube-in-shell heat exchanger: (a) Exchanger measurements

and (b) calculation schematic ..................................................................................... 51

27. Liquid-ullage interface in MHTB tank with 90% fill level .......................................... 52

28. Interface and dry wall areas for hydrogen in MHTB tank: (a) 25% fill level
and (b) 90% fill level ................................................................................................... 53

29. Error message example ............................................................................................... 56

30. Custom error handling code segment ......................................................................... 56

31. Standard error handling code segment ....................................................................... 57

32. Program flowchart for regular mission phase ............................................................. 58

33. Initialization subroutine flowchart ............................................................................. 59

34. Heat exchanger transfer subroutine flowchart ............................................................ 60

35. TankSIM file structure ............................................................................................... 61

36. Ullage pressure, 25% fill level...................................................................................... 63

37. Liquid saturation pressure, 25% fill level..................................................................... 64

38. Ullage temperature, 25% fill level................................................................................ 64

39. Liquid temperature, 25% fill level................................................................................ 64

40. Ullage wall temperature, 25% fill level......................................................................... 65

vii
 LIST OF FIGURES (Continued)

41. Ullage pressure, 50% fill level...................................................................................... 66

42. Liquid saturation pressure, 50% fill level..................................................................... 66

43. Ullage temperature, 50% fill level................................................................................ 66

44. Liquid temperature, 50% fill level................................................................................ 67

45. Ullage wall temperature, 50% fill level......................................................................... 67

46. Ullage pressure, 90% fill level...................................................................................... 68

47. Liquid saturation pressure, 90% fill level..................................................................... 69

48. Ullage temperature, 90% fill level................................................................................ 69

49. Liquid temperature, 90% fill level................................................................................ 70

50. Ullage pressure, adjusted heat loads, 25% fill level ...................................................... 71

51. Liquid saturation pressure, adjusted heat loads, 25% fill level...................................... 71

52. Ullage temperature, adjusted heat loads, 25% fill level................................................. 72

53. Liquid temperature, adjusted heat loads, 25% fill level................................................ 72

54. Ullage wall temperature, adjusted heat loads, 25% fill level......................................... 72

55. Ullage and liquid saturation pressures (cycling ullage venting) ................................... 73

56. Ullage and liquid temperatures (cycling ullage venting) .............................................. 74

57. Ullage and liquid temperatures (cycling ullage venting, large scale) ............................ 74

58. Cylindrical to flat tank wall transformation ............................................................... 84

59. Relative difference between heat flow rates through circular and flat tank walls ......... 85

60. Minimal radius-to-height ratio value allows using the flat wall convective heat
flow model for cylindrical walls .................................................................................. 87

61. Relative Nusselt numbers difference versus radius-to-height ratio

for different Rayleigh numbers ................................................................................... 88

viii
 LIST OF FIGURES (Continued)

62. Condensation on the internal side of cylindrical wall segment ................................... 89

63. Condensation on the upper dome .............................................................................. 92

64. Angle a calculation schematic .................................................................................... 94

65. Condensation on the lower dome ............................................................................... 96

66. Heat and mass transfer in the liquid film .................................................................... 97

67. Geometrical parameters of axial jet ........................................................................... 100

68. Geometrical parmeters of helicoidally coiled tube-in-shell heat exchanger ................. 101

69. TankSIM address input file ........................................................................................ 105

70. TankSIM input file ..................................................................................................... 106

71. TankSIM mission profile input file ............................................................................. 106

72. TankSIM tank input file ............................................................................................. 108

73. TankSIM spray bar input file ..................................................................................... 109

74. TankSIM axial jet input file ........................................................................................ 110

75. TankSIM external venting line input file .................................................................... 111

76. TankSIM main output file .......................................................................................... 111

ix
 LIST OF TABLES

1. Recommended values for constant C in Lockhart-Martinelli parameter .................... 32

2. Heat leaks in 25% fill level test, adjusted heat loads .................................................... 70

3. Coefficient K from equation (165) .............................................................................. 83

4. Variables description for input and output files........................................................... 112

x
 LIST OF ACRONYMS

CFD computational fluid dynamics

J-T Joule-Thomson

MHTB multipurpose hydrogen test bed

NIST National Institute of Standards and Technology

PCS pressure control system

RefProp NIST reference fluid thermodynamic and transport properties database

TankSIM tank system integrated model

TM Technical Memorandum

TVS thermodynamic venting system

xi
 NOMENCLATURE

Areas, m2
Acw tank wall cross area on the ullage-liquid interface level
Alw wall-liquid interface
Au ullage surface area for low gravity cases
Auw wall-ullage interface
Awl wall liquid-ullage interface
Bo Bond number, dimensionless
Cp constant pressure specific heat, J/(kg·K)
Cv constant volume specific heat, J/(kg·K)
Dc average diameter of coil, m
Dcr critical bubbles diameter, m
Dd bubble departure diameter, m
Ddr ullage droplets diameter, m
Dhex heat exchanger hydraulic diameter, m
Djnz axial jet nozzle internal diameter, m
Dmi spray manifold internal diameter, m
Dsh internal diameter of coiled heat exchanger shell, m
e average internal energy, J/kg
G mass flux through cross-section of the heat exchanger, kg/(m2·s)
gr system acceleration ratio, dimensionless
h average enthalpy, J/kg

Heat Transfer Coefficients, W/(m2·K)

hboil boiling heat transfer coefficient inside heat exchanger
hcond ullage wall condensation heat transfer coefficient
hconv two-phase flow forced convection heat transfer coefficient inside heat
exchanger
hdrimp heat transfer coefficient from a tank wall to impingent droplets
hexm two-phase flow heat transfer coefficient inside the heat exchanger
hmex forced convection heat transfer coefficient inside spray manifold
hmext total heat transfer coefficient from spray manifold to heat exchanger

xii
 NOMENCLATURE (Continued)

Enthalpies, J/kg
hdrv vapor by liquid droplets temperature and pressure
hintv vapor by interface temperature and pressure
hlv vapor by liquid temperature and pressure
hsthe pressurization non-condensable gas from pressurization system
hstv pressurization vapor from pressurization system
huhe non-condensable gas in ullage by ullage pressure and partial pressure
huv vapor in ullage by ullage temperature and vapor partial pressure
hwlw vapor by wall liquid temperature and pressure
hfg latent heat on the liquid-ullage interface, J/kg

Others
Ja Jacobs number, dimensionless
k thermal conductivity, W/(m·K)
kmex spray manifold wall thermal conductivity heat transfer coefficient, W/(m2·K)
l* linear length scale for wall conduction heat transfer, m
M, m mass, kg

Mass Flow Rates, kg/s

m! mass flow rate
m! chill liquid taking from the tank for pipe system chilldown
m! cnj condensation initiated by axial jet
m! drboil droplets boil-off to ullage

m! drl unevaporated droplets falling to liquid

m! drwl unevaporated droplets move to wall liquid

m! evap evaporation from liquid interface to ullage
m! fir liquid taking from the tank for engine firing
m! hepress noncondensable pressurant to the ullage
m! hevent noncondensable gas vented from the ullage
m! lboil bulk liquid boil-off

xiii
 NOMENCLATURE (Continued)

m! lvent vented liquid

m! mixl total from injection pipe to bulk liquid
m! mixu total from injection pipe to ullage (droplets)
m! orif through spray bar orifice
m! pump liquid pumped from tank to thermodynamic venting system
m! ucndr bulk vapor condensation on the droplet

m! ucnl vapor condensation on the ullage bulk liquid interface
m! ucnwl vapor condensation on the wall liquid interface
m! uvent total ullage venting
m! vpress autogenous pressurant mass flow rate
m! vvent vapor vented from the ullage
m! wboil bulk liquid boiling on the wall boil-off
m! wlboil wall liquid boil-off
m! wll from wall liquid to bulk liquid
mfhe noncondensable gas mass fraction in ullage, dimensionless
mfv liquid vapor mass fraction in ullage, dimensionless
P pressure, Pa
Psat saturation pressure, Pa
Pr Prandtl number, dimensionless
Prl liquid Prandtl number, dimensionless
Prtp two-phase flow Prandtl number, dimensionless

Heat Flow Rates, W

qaddl additional heat flow from feed, pressurization, and venting lines to the liquid
qaddu additional heat flow from feed, pressurization, and venting lines to the ullage
qelw from environment to the liquid wall
qenl from environment to the liquid
qenu from environment to the ullage
qeuw from environment to the ullage wall
qintl from ullage-liquid interface to bulk liquid

xiv
 NOMENCLATURE (Continued)

qj axial jet volumetric flow rate, m3/s

qjcnl condensation mass flow through gas-liquid interface, initiated by axial jet
qlhex from liquid to the heat exchanger
qlit from liquid to the injection pipe
qltot total heat flow to the liquid from all sources
qlwct constant heat flow from construction elements to the liquid wall
qlwl from liquid wall to bulk liquid
qmg from the wall to the ullage in microgravity
qmhex from spray manifold to heat exchanger
qudr from ullage to droplets
quhex from ullage to the heat exchanger
quint from the ullage to the ullage-liquid interface
quit from ullage to the injection pipe
qutot total heat flow to the ullage from all sources
quwct constant heat flow from construction elements to the ullage wall
quwl ullage to wall liquid
quwu from ullage wall to ullage
quwlw conduction from the ullage wall to the liquid wall
quwwl from the ullage wall to the wall liquid
′
qunif heat flux from the environment, uniformly distributed on the tank
surface, W/m2
Retp two-phase flow Reynolds number, dimensionless
S* conduction shape factor, dimensionless
Sdr speed of droplets created by spray bar jets in ullage, m/s
s axial jet nozzle submergence, m
T temperature, K
Tenv environmental temperature, K
Tsat saturation temperature, K
U bubbles terminal velocity, m/s
V volume, m3
Vf volume available for the flow of fluid in the annulus, m3
Vtot total internal tank volume, m3
Vu ullage volume, m
X vapors quality inside heat exchanger, dimensionless

xv
 NOMENCLATURE (Continued)

Xlm two-phase flow Martinelli-Nelson parameter, dimensionless

γ coiled tube pitch, dimensionless
∆To axial jet nozzle outflow liquid subcooling, K
m viscosity, kg/(m·s)
mw liquid viscosity at wall temperature, kg/(m·s)
r density, kg/m3
rave average density for homogeneous fluid in macrogravity, kg/m3
σ Stefan-Boltzmann constant, W/(m2·K4)
σl surface tension, kg/s2
clw liquid wall heat experimental absorption coefficient for multilayer MLI,
dimensionless
cuw ullage wall heat absorption experimental coefficient for multilayer MLI,
dimensionless

Subscripts
ave average value
b bulk liquid
dr droplet in the ullage
drimp impinged droplet
he noncondensable gas
hest noncondensable gas storage
int bulk liquid-ullage interface
it injection tube
j axial jet
l liquid
lw liquid wall
orif spray bar orifice
tp two-phase
u ullage
uhe ullage noncondensable gas
uv ullage vapor
uw ullage wall
v vapor
w tank wall
wl wall liquid

xvi
 TECHNICAL MEMORANDUM

TANK SYSTEM INTEGRATED MODEL: A CRYOGENIC TANK

PERFORMANCE PREDICTION PROGRAM

1. INTRODUCTION

Accurate prediction of the thermodynamic state of the cryogenic propellants in launch

vehicle tanks is necessary for mission planning and successful execution. Cryogenic propellant stor-
age and transfer in space environments require that tank pressure be controlled. The pressure rise
rate is determined by the complex interaction of external heat leak, fluid temperature stratification,
and interfacial heat and mass transfer. If the required storage duration of a space mission is longer
than the period in which the tank pressure reaches its allowable maximum, an appropriate pressure
control method must be applied. Therefore, predictions of the pressurization rate and performance
of pressure control techniques in cryogenic tanks are required for development of cryogenic fluid
long-duration storage technology and planning of future space exploration missions.

This Technical Memorandum (TM) describes an analytical tool, Tank System Integrated
Model (TankSIM), which can be used for modeling pressure control and predicting the behavior of
cryogenic propellant for long-term storage for future space missions. TankSIM is written in Fortran
90 language and can be compiled with any Visual Fortran compiler. A thermodynamic vent system
(TVS) or ullage venting is used to achieve tank pressure control. Utilizing TankSIM, the following
processes can be modeled: tank self-pressurization, boiloff, ullage venting, and mixing.

The intended application of TankSIM is for performance prediction for cryogenic propel-
lant storage systems during in-space missions such as Mars and Lunar exploration and rendezvous
with asteroids and near-Earth objects, as well as ground-based thermal vacuum testing, under
different levels of system accelerations, external heat loads, and various thermal initial conditions.
Each mission profile can consist of several different phases including coast periods, propellant
transfer system chilldown, propellant tank pressurization, propellant settling, and main engine
operation. Each of these phases can have different accelerations, external heat loads, and durations,
and they can be repeated as many times as necessary in any sequence to create a full mission
timeline.

The TankSIM analytical model is described in section 2. Section 3 covers the basic mass
and energy conservation equations used in this program. The mass and heat transfer in a locked-up
tank during nonmixing regimes is presented in section 4. In section 5, the tank pressurization
and ullage venting pressure control are described. Section 6 introduces the spray bar TVS. An
axial jet TVS pressure control system is presented in section 7. In section 8, the computer model
for TankSIM is described. Section 9 covers comparison between TankSIM predictions and test
data as well as examples of full mission calculations. Discussions of results and plans for future
improvements of TankSIM are provided in section 10.
1
 2. ANALYTICAL MODEL DESCRIPTION

TankSIM is a lumped, multinode transient model based on mass and energy conservation
equations and correlation equations for mass and heat transfer calculations.

2.1 Basic Model Approaches for a Settled Liquid

2.1.1 Model Nodes

The TankSIM network with heat and mass transfer between the nodes is depicted in figure 1
with the following 13 nodes:

(1) Ullage.
(2) Ullage wall—upper head part (upper dome) of the tank wall interfaced to the ullage.
(3) Ullage wall—cylindrical part of the tank wall interfaced to the ullage; appears at low
and intermediate liquid level.
(4) Ullage wall—segment of the lower dome of the tank wall interfaced to the ullage
at level, lower than dome maximum height.
(5) Bulk liquid.

8 13 1

12
3
11

10

6
5

Figure 1. TankSIM nodes.

F1_1627

2
 (6) Liquid wall—bottom head part (lower dome) of the tank wall. interfaced to the liquid.
(7) Liquid wall—cylindrical part of the tank wall interfaced to the liquid.
(8) Liquid wall—segment of upper dome of the tank wall interfaced to the liquid when
filling is above the cylindrical part of the tank.
(9) Environment.
(10) Bubble in bulk liquid during boiling of the bulk liquid on the wall and in volume.
(11) Ullage—liquid interface.
(12) Wall liquid—liquid film on the ullage tank wall.
(13) Droplet in ullage during spray bar operation.

2.1.2 Control Volumes and Main Parameters

There are seven control volumes in the TankSIM model as shown in figure 2:

(1) Ullage—mixture of propellant vapors and noncondensable pressurization gas. Special

assumptions for the ullage were made:
– Temperature and volume of vapor and noncondensable gas are the same and equal
to the overall tank temperature and volume.
– Total ullage pressure is a sum of vapor and noncondensable gas pressures (Dalton’s law).
– Total ullage mass is a sum of vapor and noncondensable gas masses.
(2) Liquid—liquid propellant in the tank.
(3) Ullage-liquid interface—infinitely thin layer with its own temperature, containing no
mass and used for the heat and mass transfer between liquid and ullage (evaporation-condensation
processes) at the interface.

Droplets
Wall Liquid

Ullage

Ulage Wall

Interface

Liquid Wall

Liquid

Figure 2. Diagram of control volumes. F2_1627

3
 (4) Ullage wall—part of the tank wall interfaced to the ullage.
(5) Liquid wall—part of the tank wall interfaced to the liquid.
(6) Wall liquid—liquid film on the ullage wall.
(7) Droplets in the ullage during operation of the spray bar TVS.

2.1.3 Mass Conservation Equations

As shown in appendix A, the mass conservation equations can be written as

dm
= ∑ m! k , (1)
dt k

where m is mass of control volume or node, kg, and m! k is mass inflow and outflow rates, positive
for inflow and neagative for outflow, kg/s.

2.1.4 Energy Conservation Equations

As shown in appendix A, the internal energy and enthalpy equations for the control volume
without dissipation can be written as:

d(e ) dV
  m
dt
( )
= ∑ qk + ∑ m! k hk − e − P
dt
(2)
k k
and
d(h ) dP
m
dt
( )
= ∑ qk + ∑ m! k hk − h +V
dt
, (3)
k k
where

m – mass of control volume, kg

e – average internal energy of the control volume, J/kg
qk – heat flow rate, not associated with in and out mass flows, W
m! k – in and out mass flow rates, positive for inflow and negative for outflow, kg/s
h – average enthalpy of the control volume, J/kg.

Detailed equations for each control volume and node will be presented in section 3.

4
 2.2 Model Approaches for an Unsettled Liquid

2.2.1 Nonmixing Cases

For calculating heat-mass transfer in unsettled liquids during nonmixing phases, a concen-
tric spherical enclosure model is used. In this model, the outside and inside spheres represents the
tank shell and ullage, respectively, and the space between these spheres is a liquid layer. The vol-
umes of the outer and inner spheres are equal to tank and ullage volumes, and the volume of the
liquid layer is equal to the bulk liquid volume, as depicted in figure 3.

Wall

Ullage

Liquid

Figure 3. Tank geometry transformation for unsettled liquid (no mixing).

F3_1627

2.2.2 Mixing Cases

For an unsettled liquid in a mixing case, a homogeneous fluid filling the entire tank volume
was assumed (fig. 4). Thermodynamic properties such as density, thermal conductivity coefficient,
specific heat capacity, and others where calculated as mass weighted values by the following
equation:

Cphe mhe + Cpv mv + Cpl ml

Cpmix = . (4)
mhe + mv + ml

Ullage
Homogeneous
Fluid
Liquid

Figure 4. Tank fluids transformation for unsettled liquid during mixing.

F4_1627

5
 3. MASS AND ENERGY CONSERVATION EQUATIONS

3.1 Ullage

To determine the ullage mass, temperature, and pressure, conservation of mass and energy
equations (5) and (6) are utilized.

3.1.1 Ullage Mass Conservation

The ullage mass conservation equation with respect to figure 5 can be written as:

dmu
= m! evap + m! lboil + m! wboil + m! wlboil − m! ucnwl
dt

( )
+ m! drboil − m! ucndr − m! cnj − m! vvent + m! hevent + m! vpress + m! hepress . (5) ( )
m⋅ hepress
m⋅ vpress

m⋅ vvent
m⋅ drboil
m⋅ hevent m⋅ drwl

m⋅ ucndr m⋅ wlboil
m⋅ mixu
m⋅cnj
m⋅ drl m⋅ ucnwl
m⋅ evap
m⋅ ucnl

m⋅ lboil m⋅ wll
m⋅ wboil

m⋅ mixl

m⋅ pump
m⋅ fir
m⋅ chil m⋅ lvent

Figure 5. Diagram of mass flows.

F5_1627

6
 dmu dmv dmhe
Because mu = mv + mhe , then = + , or m! u = m! v + m! he .
dt dt dt

So, equation (5) can be rewritten as:

m! u = m! evap + m! lboil + m! wboil + m! wlboil − m! ucnl − m! ucnwl

+ m! drboil − m! ucndr − m! cnj − m! vvent − m! hevent + m! vpress + m! hepress . (6)

3.1.2 Ullage Energy Conservation

As depicted in figures 5 and 6, the energy conservation equation for the ullage follows as:

mu
( )=q
d eu
+ quwu − qu int − quwl − quit − qudr − quhex
enu
dt

dV
+ ∑ m! k hk − eu ∑ m! k − Pu u . (7)
k k
dt

qenu
dVwl
Pwl
dt
dVu
Pu dVdr
dt Pdr
dt
qudr qeuw

quwu quwl
quit
quwwl
quint quhex
quwlw

qintl qenl
qelw
qlhex dVl
Pl
qlwl qlit dt

Figure 6. Diagram of heat flows, not associated with mass transfer.

F6_1627

7
In equation (7), ∑ m! k = m! u is total ullage mass change rate and
k

∑ m! k hk = m! evaphint v + m! wboil hlv + m! lboil hlv + m! wlboil hwlv + m! drboil hdrv

k

+ m! hepress hsthe + m! vpress hstv − m! ucnwl huv − m! ucndr huv − m! cnj huv

( )
− m! vvent huv + m! hevent huhe . (8)

All terms described by equation (8) are not shown in figure 6.

Assuming vented gas has the same composition as the ullage (uniformly distributed species
within the ullage), defining mass fractions for the vapor, mfv, and noncondensable gas, mfhe, can be
written as

mfhe = mhe / mu = m! hevent / m! vent , then mfv = mv / mu = m! vvent / m! vent . (9)

Rewriting the energy balance for the venting gas gives the following:

m! vvent huv + m! hevent huhe

m! vvent huv + m! hevent huhe = m! uvent
m! uvent
( )
= m! uvent mfv huv + mfhe huhe . (10)

Note, that mfv + mfhe = 1.

Similarly, average ullage internal energy, eu , can be written as:

mv euv + muhe euhe

eu =
mu
(
= mfv euv + mfhe euhe . (11) )

3.2 Ullage Wall (Tank Wall Interfaced to the Ullage)

3.2.1 Ullage Wall Energy Conservation

For the tank wall interfacing to the ullage, the following can be written from the energy
conservation equation (see fig. 6):

( )=
d euw
muw
dt
∑ qn (12)
n
or

8
 muw
( )=q
d euw
− quwwl − quwu − quwlw . (13)
euw
dt

3.3 Bulk Liquid

3.3.1 Bulk Liquid Mass Conservation Equation

The mass conservation equation for a bulk liquid is as follows (see fig. 5):

dml
= − m! evap − m! lboil − m! wboil + m! wll + m! drl + m! mixl − m! pump − m! fir − m! chill − m! lvent . (14)
dt

3.3.2 Bulk Liquid Energy Conservation Equation

For a bulk liquid, the energy equation can be written as:

ml
( )=q
d el dV
+ qlwl + qintl − qlit − qlhex + ∑ m! i hi − el ∑ m! i − Pl l , (15)
enl
dt i i
dt

where ∑ m! i = m! l .
i

From figure 5,

∑ m! i hi = − m! evaphintl − m! lboil hl − m! wboil hl + m! wll hwl + m! drl hdrl

i

+ m! mixl hmixl + m! cnj hcnj − m! fir hl − m! chil hl − m! vent hl − m! pumphl . (16)

3.4 Liquid Wall (Tank Wall Interfaced to the Liquid)

3.4.1 Liquid Wall Energy Conservation

From figure 6, the energy conservation equation for the tank wall interfacing to the liquid
can be written as:

( )=
d elw
mlw
dt
∑ qn (17)
n
or

mlw
( )=q
d elw
+ quwlw − qlwl . (18)
elw
dt

9
 3.5 Wall Liquid (Liquid Film on the Ullage Wall)

3.5.1 Wall Liquid Mass Conservation

The wall liquid mass conservation equation follows:

dmwl
= m drwl + m ucnwl − m wlboil − m wll . (19)
dt

3.5.2 Wall Liquid Energy Conservation

In accordance with figures 5 and 6, the energy conservation equation can be written
as follows:

mwl
( )=q
d ewl dV
+ quwl + m! ucnwl hul + m! drwl hdr − m! wlboil hwl − ewl ∑ m! i − Pwl wl , (20)
uwwl
dt i
dt

where ∑ m! i is equal to the right side of equation (19).

i

10
 4. HEAT AND MASS TRANSFER INSIDE THE TANK

4.1 External Heat Loads

External heat loads are heat flows from the environment to the tank. The following five
external heat load mechanisms are incorporated into the TankSIM:

(1) Radiation to tank walls, calculated by the Stefan-Boltzmann law and corrected
by the average heat absorption coefficient, used for the tank thermal insulation.
(2) Heat flux uniformly distributed on the external surface of the tank wall.
(3) Constant heat flows to the ullage and to the liquid from supporting elements.
(4) Constant heat flows to the ullage wall and liquid wall from supporting elements.
(5) Heat flows to the tank walls from venting, filling, and pressurization lines. Heat flow
from the environment to the ullage wall (wall interfaced to the ullage) calculates by:

qeuw = ⎡σχ uw Tenv

⎣
4
( 4
− Tuw )
′ ⎤ Auw + quwct + qaddu . (21)
+ qunif
⎦

Similarly, heat flow from the environment to the liquid wall (tank wall interfaced to the
liquid) is as follows:

⎣ (
qelw = ⎡σχ lw Tenv
4 4
− Tlw )
′ ⎤ Alw + qlwct + qaddl . (22)
+ qunif
⎦

4.2 Ullage-Bulk Liquid Interface Heat and Mass Transfer

4.2.1 Interface Temperature

Control volume representing ullage—liquid interface with no mass has only temperature,
which is used for calculating heat and mass transfer across the interface. Temperature of the inter-
face is calculated by the modified Alabovskii’s1 equation (detailed description is provided in app. B):

KTl + Tu Tl + Tu K
Tint = = , (23)
K +1 1+ 1 K
where

0.5
k ⎛µ ρ ⎞
K = l ⎜ u l Prl ⎟ (1+ m)
ku ⎝ µl ρu ⎠

qevap
m= is the ratio of evaporation heat to total heat flows through interface.
qintl

11
4.2.2 Heat and Mass Transfer at the Liquid-Ullage Interface

Mass of evaporated/condensed liquid/vapor from the ullage—liquid interface is calculated

using the energy jumping boundary condition suggested by Delhaye2 and Meserole et al.3:

quint − qintl
m! evap = . (24)
( ) (
hv Tint ,Pu − hl Tint ,Pu )
In this equation, positive/negative mass flow rate corresponds to the evaporation/condensation
from/to the interface. Both qu int and qint l flow rates are calculated by the heat conduction model
with assumptions that temperatures of the ullage and bulk liquid are linearly depend on the height:

(
k T −T ) k T −T
quint = 2 Aint u u int , qintl = 2 Aint l int l . (25)
( )
hu hl

The conduction model was chosen because there are no correlations for liquid-gas interface
heat flow rates in a finite volume. Usually, the Horizontal Heated Downward-Facing Plate model is
used (see ref. 4). It is based on the finite flat plate, heated from the bottom, in infinite space (fig. 7),
which is far from real conditions (fig. 8).

Figure 7. Horizontal heated downward-facing

F7_1627 plate model.

Figure 8. Liquid-ullage interface in finite volume tank.

F8_1627

12
 4.3 Ullage Wall Heat and Mass Transfer

The ullage wall is divided into three segments—the cylindrical segment and the upper
and lower domes.

4.3.1 Ullage Wall Convection for the Cylindrical Segment

The heat transfer coefficient for the cylindrical part is calculated using the Churchill and
Usagi5 blending equation, applied to both laminar and turbulent convection regimes. Assuming the
tank wall is a flat plate model (see app. B.3),

6 ⎤ 1/6
⎡
(
6
) (
Nu = Nul + Nut
⎢⎣ ⎥⎦ ) , (26)

where

2
Nul = (Nusselt number for laminar convection)
(
ln ⎡1+ 2 Cl Ra 1/4 ⎤
⎣ ⎦ )
0.671
Cl = 4/9
(Correction coefficient for laminar Nusselt number)
⎡1+ ( 0.492 Pr )9/16 ⎤
⎣ ⎦

Nut =
(
f Tuw Tu Ctv Ra1/3) (Nusselt number for turbulent convection)
(
1+ 1.4 × 109 Pr Ra )
0.13Pr 0.22
Ctv = (Correction coefficient for turbulent Nut)
( )
0.42
1+ 0.61Pr 0.81

⎛T ⎞ ⎛T ⎞
f ⎜ uw ⎟ = 1+ 0.078 ⎜ uw − 1⎟ (Correction coefficient for high wall-ullage temperature
⎝ Tu ⎠ ⎝ Tu ⎠ difference).

4.3.2 Upper and Lower Domes (Endcups)—Ullage Heat Transfer

Due to the lack of appropriate correlations, a heat conduction mode was used to calculate
heat flow between the upper and lower domes and the ullage:

2 Adu ku
qdu =
hu
( )
Tuw − Tu . (27)

13
 4.4 Condensation on the Ullage Wall

The condensation calculation for the ullage wall is based on the Gerstmann and Griffith 6,7
model, and two different modifications of the Nusselt vertical flat plate condensation model. Details
of the calculations are provided in appendix B.

4.4.1 Cylindrical (Barrel) Segment of the Tank

The condensation heat flow rate on the vertical cylindrical ullage wall calculates using
the Nusselt model with a Rohsenow modification 4, reintegrated for nonzero initial liquid film
thickness. With this modification, the average Nusselt number, Nu, is as follows:

1/4
(
⎡ g ρ ρ − ρ h′ L3 ⎤ )
Nu = 0.9428 ⎢
l l
(
v fg
⎢ kl µl Tsat − Tw ⎥
⎥
)
⎡
⎢⎣ ( b +1 − 2 b )( b + b +1 )1/2 + 2b3/4 ⎤ . (28)
⎥⎦
⎣ ⎦

The liquid film thickness, δ, can be calculated as:

1/4

⎛ 4µ k T − T L ⎞
δ = ⎜ l l sat w ⎟
( ) ⎛
b + b +
x⎞
1/2
, (29)
(
⎜⎝ ρl ρl − ρ v gh′fg ⎟⎠ ) ⎜⎝ L ⎟⎠

where x is a coordinate along the vertical surface as shown in figure 62 in appendix B, and L
is a height of the cylindrical segment

In equations (28) and (29), the dimensionless parameter, b, as a function of initial liquid
film thickness, δ0, is as follows:

b=
(
δ 04 ρl ρl − ρv gh′fg
. (30)
)
16 µl kl Tsat − Tw L ( )
4.4.2 Upper Dome (Endcup) Segment Interfaced to the Ullage

For the mass and heat transfer on the upper dome, the Gerstmann and Griffith6,7 condensa-
tion model on the bottom side of the inclined flat plate was used. Applying their model, the upper
dome is divided into two regions with respect to angle a: 0° ≤ a ≤ 20° and 20° < a ≤ 90°.

For 0° ≤ a ≤ 20°, the Nusselt number is:

0.90
Nu = , (31)
( )

Ra1/6 1+ 1.1Ra1/6

14
 1/2
ρlσ h fg ⎡ σ ⎤
where Ra = ⎢ ⎥ is a modified Rayleigh number.
( ) (
kl µl Tsat − Tw ⎢⎣ g ρl − ρv cos(α ) ⎥⎦ )

For 20° < a ≤ 90°, a modified Nusselt model was recommended6. In this model, the free-
falling acceleration, g, should be replaced by gsin(a) as shown in the following equation:

1/4
(
⎡ g sin(α ) ρ ρ − ρ h′ L3 ⎤ )
Nul = 0.9428 ⎢
⎢
l l
(
v fg
kl µl Tsat − Tw
⎥
⎥ )
⎡
⎣ ( b +1 − 2 b )( )
b + b + 1 1/2 + 2b 3/4⎤ , (32)
⎦
⎣ ⎦
1/4

⎛ 4µ k T − T L ⎞
δ =⎜ l l sat w ( ⎟
) ⎛
b + b +
x⎞
1/2
, (33)
( )
⎜⎝ ρl ρl − ρv g sin(α )h′fg ⎟⎠ ⎜⎝ L ⎟⎠
and

b=
(
δ 04 ρl ρl − ρv g sin(α )h′fg )
. (34)
16 µl kl Tsat − Tw L ( )
For the heat and mass transfer rates calculation, the upper dome is divided into layers
as depicted in figure 9. Each layer (for example, ABCD) has its own angle a, side area, Rayleigh
number, and liquid film thickness. The final liquid film thickness from the previous layer is used
as an initial thickness for the next one.

h1
hd B C
A D
α
x
g

Figure 9. Calculation schematic for condensation on the upper dome.

F9_1627

Angle a can be found as follows:

h
tg(α ) = d
(
h1 2 − h1 ) , (35)
r (1− h1)2

15
where
h
h1 = 1 .
hd

4.4.3 Lower Dome (Endcup)—Ullage Interface

As shown by Oosthuizen8, the condensation on the upper side of the inclined plate can be
described by modified Nusselt equations, where the free-falling acceleration is replaced by gsin(b).
As depicted in figure 10, the angle b is an angle between a tangent line and free-fall acceleration.
The condensation on the lower dome can be calculated as follows:

1/4
(
⎡ g sin( β ) ρ ρ − ρ h′ L3 ⎤ )
Nul = 0.9428 ⎢
⎢
l l
(
v fg
kl µl Tsat − Tw
⎥
⎥ )
⎡
⎢⎣ ( b +1 − 2 b )( b + b +1 )1/2 + 2b3/4 ⎤ (36)
⎥⎦
⎣ ⎦
and
1/4

⎛ 4µ k T − T L ⎞
δ =⎜ l l sat ( w
⎟
) ⎛
b + b +
x⎞
1/2
. (37)
( )
⎜⎝ ρl ρl − ρv g sin( β )h′fg ⎟⎠ ⎜⎝ L ⎟⎠

x
h1
hd

β
g

F10_1627
Figure 10. Calculation schematic for condensation on the lower dome.

In equations (36) and (37), parameter b is:

δ 4 ρl ( ρl − ρv ) g sin( β )hʹfg Cp (Tsat −Tw )

b= 0 , hʹfg = h fg + 0.68Ja, and Ja = l ,
16 µl kl (Tsat −Tw ) L h fg
(38)

where L is a height of the liquid inside the lower dome.

The calculation for the lower dome is similar to that for the upper dome, as described earlier.

16
 Angle b is described as follows:

tg( β ) =
r (1− h1)2 ,
hd h1( 2 − h1) (39)

h
where h1 = 1 .
hd
4.5 Bulk Vapor Condensation

Bulk vapor condensation occurs when the ullage temperature is lower than the vapor satu-
ration temperature that corresponds to the vapor partial pressure, Tu < Tsat(Pv). For this case, the
condensation mass flow rate is:

qutot
m! cond = . (40)
( ) (
hv Tu , Pv − hl Tu , Pv )

4.6 Bulk Liquid Convective Heat and Mass Transfer

For wall-liquid heat transfer modeling, the wall is divided into three segments—the cylindri-
cal segment and upper and lower domes.

4.6.1 Convection on the Cylindrical Segment of the Liquid Wall

The heat transfer coefficient for the cylindrical part of the liquid wall is calculated using the
same Churchill and Usagi5 blending equation as was applied for ullage-ullage wall (eq. (26)). In
this case, the Rayleigh and Prandtl numbers and correction coefficients are calculated using liquid
properties and liquid and wall temperatures.

4.6.2 Heat Transfer on the Lower Dome (Endcup) Interfaced to the Liquid

For heat transfer between the lower dome and bulk liquid, the horizontal heated upward
-facing plate with uniform wall and liquid temperatures model, described by Rohsenow4, is used:
10 ⎤ 1/10
(
Nu = ⎡ Nul
⎣⎢
)10
(
+ Nu t ) ⎦⎥
, (41)

where
1.4
Nul = – laminar Nu ,

⎣ (
ln ⎡1+ 1.4 0.835Cl Ra1/4 ⎤
⎦)

17
 ⎛ 1+ 0.0107 Pr ⎞
Nut = 0.14 ⎜ ⎟ Ra1/3 – turbulent Nu ,
⎝ 1+ 0.01 Pr ⎠
and
0.671
Cl = 4/9
⎡1+ (0.492 Pr)9/16 ⎤ – correction coefficient .
⎣ ⎦

As suggested by Goldstein et al.9, the length scale accepted in that model is the ratio
Awet /p, where Awet is a lower dome wetted area and p is a perimeter of a wetted segment (see
fig. 11).

Liquid Wall
Temperature Temperature Perimeter, p

Heated Plate

Wetted Area, Awet

Figure 11. Definition sketch for lower dome convective model.

F11_1627

4.6.3 Upper Dome (Endcup)—Liquid Heat Transfer

For the upper dome wetted segment, the modified Churchill and Usagi model (eq. (26))
can be used. In this model, as suggested by Rich10, free-fall acceleration g is replaced by g cos(ϕ),
where angle ϕ is the angle between the direction of acceleration and the upper dome tangent line as
depicted in figure 12. After replacing the free-fall acceleration, the Rayleigh number is:

Ra =
(
g cos(ϕ )βl ρl2Cpl Tl w − Tl L3 )
, (42)
kl µl
where L is a height at liquid in the upper dome.

h1 y

hd ϕ

x
g
r

Figure 12. Calculation schematic for convection in the liquid on the upper dome wall.
F12_1627

18
Angle ϕ can be found as follows:

tg(ϕ ) =
r (1− h1)2 and cos(ϕ ) = 1
, (43)
h1( 2 − h1)

hd 1+ tg2 (ϕ )

h
where h1 = 1 .
hd

4.7 Bulk Liquid Boiling

The liquid starts to boil when its saturation pressure is higher than the liquid pressure,
Psat > Pl.

4.7.1 Bulk Liquid Boiling on the Tank Wall (Heterogeneous Boiling)

The liquid is boiling on the wall when the temperature of the wall is higher than the liquid
saturation temperature, Tlw > Tsat(Pl ). To describe this process, a correlation containing the heat
flux from the wall to the liquid, including the dependency of the wall material properties, was
accepted. As reported by Grigoryev et al.11, there are two correlations that depend on thermal
properties of the liquid and the wall. The first correlation is used for nitrogen, oxygen, methane,
and other cryogenic liquids (excluding hydrogen, neon, and helium). In this case,

Qboil A13 ⎛ h fg ρv ΔT0 ⎞ ⎛ kl ρl ΔT ⎞

= C + C ρ
2 fg v ⎟ ,
h (44)
Alw A1 + 2m ⎜⎝ σ Tsat ⎟⎠ ⎜⎝ 1 µl ⎠
where
Qboil Alw = heat flux from the wall to the liquid

A1 = −m + m2 + γ M

C 1 = 8.010 −5

C2 = 10 −4.

The second correlation is suggested for low point boiling cryogens, such as hydrogen,
parahydrogen, neon, and helium, where klCpl ρl ≈ kwCpw ρw :

Qboil ⎛ h fg ρv ΔT0 ⎞ ⎛ kl ρl ΔT ⎞
= A22 ⎜ ⎟ ⎜ C3 + C4h fg ρv ⎟ (45)
and Alw ⎝ σ Tsat ⎠ ⎝ µl ⎠

19
 2
A2 = −0.25 ( m − 2ϕ D ) + ⎡⎣ 0.25 ( m − 2ϕ D ) ⎤⎦ + 0.5 (γ M + mϕ D ) , (46)

where C3 = 3 × 10 −2 and C4 = 3 × 10 −3 . All other parameters are the same for both of the heat
fluxes11:

ΔT kwCw ρw
D= ,
h fg ρv

kl
m = 178.884 ,
kwCw ρw

kl ΔT
M =2 ,
α ′h fg ρv

ΔT0 = 0.01 K,

α ′ = 0.01,

γ = 0.45447,

ϕ = 0.44597.

4.7.2 Bulk Liquid Pool Boiling (Homogeneous Boiling)

Liquid boils in a large volume when its saturation pressure is higher than the internal
pressure, Psat > Pl. In this case, the boiling mass flow rate can be calculated as:

ql tot
m! boil = . (47)
( ) (
hv Tl ,Pu − hl Tl ,Pu )

20
 4.8 Conductive Heat Transfer Through the Wall Between the Ullage and Liquid Regions

Conduction heat transfers between the wall interfaced to the ullage and the wall interfaced
to the liquid are calculated assuming temperature distribution inside both ullage and liquid wall
segments are linear:
A k
( )
qcond = cw w Tuw − Tlw , (48)
l*
where l*= 0.5htank is a length scale and Acw is a cross-sectional area of the tank wall at the liquid
level.

4.9 Fluid Thermodynamic Properties

To calculate thermodynamic properties as functions of pressure and temperature, the

NIST RefProp FORTRAN subroutines are used12. Using these procedures significantly decreases
calculation speed (several times) compared to the interpolation from property tables without
significantly increasing calculation precision. They are used inside the TankSIM code to ensure
greater accuracy in conservation of mass and energy only. In all other cases, the property tables
are accepted. This combination of procedures and property tables allows saving calculation time
and having high precision of energy and mass conservation.

For creating property tables, two FORTRAN programs were used—one for building propel-
lant tables and the other for pressurization (noncondensable) gas tables. Pressurization gas tables
includes superheated gas properties only, but propellant tables include three parts—saturation line,
superheated vapor, and subcooled liquid. During a table creation, the user can choose the pressure
and temperature ranges and their changing steps using precisions calculated for each part of the
table. It should be noted that, for superheated vapor and subcooled liquid, instead of temperature
itself, the normalized temperatures tnmv = (tu – tvsat)/dtv for vapor and tnml = (tlsat – tl)/dtl
for  liquid are used.

In these equations, tnmv and tnml are normalized temperatures, tu and tl are ullage and
bulk liquid temperatures, respectively, tvsat is the vapor saturation temperature by vapor partial
pressure, tlsat is the liquid saturation temperature by total ullage pressure, dtv is the vapor temper-
ature range which automatically begins from the vapor saturated temperature, tvsat, and ends with
maximal temperature decided by the user, and dtl is the liquid temperature range, which automati-
cally begins from the triple point and ends at the liquid saturation temperature. As one can see,
both ranges change from 0 to 1.

The TankSIM package includes folder ‘NIST_For’ that contains all needed programs
for  building property tables. The NIST_For program also gives a possibility to use an input file or
manually put the required information into the program.

21
 5. TANK PRESSURIZATION AND ULLAGE VENTING PRESSURE CONTROL

5.1 Ullage Pressurization

There are two types of pressurization procedures: (1) Using noncondensable gas and
(2) autogenous (using propellant hot vapor). The use of each procedure depends on the propellant
and specific requirements for that particular system.

5.1.1 Pressurization With Noncondensable Gas

The model of the noncondensable gas pressurization is based on critical (choked) gas flow
through orifice or Venturi between the storage of the gas and pressurized tank. In this case, the
mass flow rate can be calculated as described in the Miller’s Handbook13, as shown in figure 13.

k+1 k+1
⎛ 2 ⎞ k−1 π d 2 ⎛ 2 ⎞ k−1
m! press = Cd A k ρ1P1 ⎜ ⎟ = Cd k ρ1P1 ⎜ , (49)
⎝ k + 1⎠ 4 ⎝ k + 1⎟⎠

where

k = Cp/Cv – heat capacity ratio on the gas storage side, dimensionless

Cd – discharge coefficient, dimensionless
A – orifice hole or Venturi throat cross-sectional area; other parameters
can be seen in figure 13, m2
d – orifice hole or Venturi throat diameter, m
l – thickness of the orifice or length of Venturi cylindrical part, m
P1 – gas pressure on the gas storage side, Pa
r1 – gas density on the gas storage side, kg/m3
V1 – volume of the noncondensable gas storage, m3.

A, d, l

Pl , Vl Pl , Vl
A, d, l ρ
ρ l
l

(a) (b)

Figure 13. Pressurization devices schematic: (a) Critical Venturi and (b) critical orifice.
F13_1627

22
 Assuming that the process is to be isentropic, and because noncondensable gas is far from
the saturation state, the downstream pressure, density, and temperature can be found by the follow-
ing equations for perfect gas13.

k 1
⎛ 2 ⎞ k−1 ⎛ 2 ⎞ k−1 2
P2 = P1 ⎜ , ρ2 = ρ1 ⎜ , and T2 = T1 . (50)
⎝ k + 1⎟⎠ ⎝ k + 1⎟⎠ k +1

With the same assumptions, the final (after pressurization) parameters of the noncondens-
able gas in storage in accordance to Spalding and Cole14 can be found:

k k−1
⎛ ρ f in ⎞ ⎛ ρf in ⎞
Pfin = Pbeg ⎜ ⎟ , Tf in = Tbeg ⎜ ⎟ . (51)
⎝ ρbeg ⎠ ⎝ ρbeg ⎠

The final density of the noncondensable gas in storage can be found using the mass flow
rate from the storage to the tank, calculated by equation (49):

⎛ m! hepress Δt⎞
ρ fin = ρbeg − ⎜ ⎟ , (52)
⎝ V1 ⎠
where

∆t = pressurization time
Pbeg, Tbeg, rbeg = initial values
Pfin, Tfin, rfin = final values
m! he press = noncondensable pressurization gas mass flow rate, kg/s.

The discharge coefficient calculates differently for critical orifice and critical Venturi. For
orifice with 1 < l/d < 6, the discharge coefficient is Cd = 0.83932 according to Ward-Smith15. For
critical Venturi, the discharge coefficient, as suggested by Simmons16, depends on the Reynolds
number which is calculated using Venturi throat diameter.
7.9139
Cd = 1− , (53)
Re
4 m!
where Re = .
πµ d

5.1.2 Autogenous Pressurization

Autogenous pressurization usually is used during the engine firing. In this case, changing of
the ullage volume is determined by the propellant flow rate needed for the engine. Pressurization
mass flow rate is calculated using partial pressures of the vapor and noncondensable gas as follows:

mhe RheTu
Phe = , (54)
Vu fin

23
where

Vubeg is an initial ullage volume, m3.

m! fir
Vu fin =Vu beg + Δt is a final ullage volume, m3.
ρl

Partial pressure of the vapor calculates as Pv = Pu – Phe; then, vapor density under new
pressure, rvfin, can be found. Finally, the mass flow rate for pressurization vapors is calculated as:

ρv finVufin − mv beg
m! vpress = . (55)
Δt

5.2 Ullage Venting Methods

Ullage venting is used for pressure control operation during long-term propellant storage.
The venting system operates between a prescribed pressure interval, namely Pmin and Pmax. Pmin
and Pmax represent minimum and maximum ullage pressures prescribed during the venting
operation for ullage and/or saturated liquid. Ullage venting is modeled for the following cases:

(1) Cycling venting.

(2) Continuous venting through a constant orifice.
(3) Continuous venting through an adjustable orifice while keeping ullage pressure at
a prescribed constant value.

5.2.1 Cycling Venting

Ullage cycling venting is modeled for the following cases:

• Ullage venting if ullage and/or liquid saturation pressures reach/reaches prescribed Pmax.
• Ullage venting if ullage pressure reaches prescribed Pmax.
• Ullage venting if liquid saturation pressure reaches prescribed Pmax.

Details of cycling venting operation logic for the three cases is shown in figure 14. The input flag
allows the user to choose one of the control logic needed to be used in each particular case.

24
 Enter

1 3
Regime

Pu > Yes Pu > Pl >

Yes
Pumax Pumax Plmax

Vent ‘ON’ Vent ‘ON’ Vent ‘ON’

No No No
No PI > Plmax

No Pu < Pumin Pu < Pumin No Pl < Plmin No

PI > Plmin Yes

No
Yes Yes Yes Yes

Vent ‘OFF’ Vent ‘ON’ Vent ‘OFF’ Vent ‘OFF’ Vent ‘OFF’

Return

Figure 14. Ullage venting pressure control logic.

5.2.2 Liquid Boiling and Liquid Level Swell During Ullage Venting

During the ullage venting process, ullage pressure drops below liquid saturation pressure
and may lead to the following conditions:

• Homogenous boiling—when ullage pressure becomes lower than the liquid saturation pressure.
• Heterogeneous boiling—when the liquid saturation temperature drops lower than the tank wall
temperature, due to homogenous boiling, which leads to boiling on the wall.

In a settled liquid, vapor bubbles increase the liquid volume and liquid level (liquid level
swell) until they reach the liquid surface and join the ullage. The time that a bubble leaves the liquid
and becomes a part of the ullage region depends on several factors, most importantly, bubble size
and system acceleration.

25
 Basic assumptions used for the liquid level swell model follow:

• Heat transfer coefficient between the wall and bulk liquid is high, so the temperature difference
between the wall and bulk liquid is small.
• The nucleate boiling takes place on the wall.
• The terminal distance for the bubbles is assumed to be equal to the liquid level.

To determine the diameter of departure bubbles, an additional assumption for the departure
bubble diameter was made. If bubbles with wall departure diameter cover an area less than the wall
area, the terminal time for bubbles is calculated using the wall departure diameter. If the covered
area is higher or equal to the wetted wall area, bubbles merge to each other and the terminal time
is  calculated using the maximum bubble diameter allowed by the Bond number. This critical diam-
eter is calculated by equation (64). To calculate the bubble-covered area, the bubbles on the surface
are assumed to have hemispherical shapes.

The bubbles departure diameter calculates using the correlation, suggested by Cole and
Rohsenow17:

Bo1/2 = 4.65 × 10 −4 ( Ja *) 5/4 (56)

and
1/2
⎡ σ ⎤
Dd = 4.65 × 10−4 ⎢ ⎥ ( Ja *) 5/4 , (57)
(
⎢⎣ g ρl − ρv ) ⎥⎦

TsatCpl ρl
where Ja* = is a modified Jacobs number
ρv h dg

Dd is a bubble departure diameter, m.

According to Levich 18, the terminal velocity of rising bubbles is dependent on the bubble
diameter and divided into three regions:

(1) Small bubbles. These bubbles have nondeformed spherical shapes and rise by straight
lines. Their diameters correspond to Reynolds numbers, Re < 1, which leads to the following:

1/3
⎛ 3µl2 ⎞
D < 2⎜ ⎟ , (58)
( )
⎝ 2g ρl − ρv ρl µ ⎠

µl + µv
where µ = .
2 µl + 3µv

26
Velocity is calculated by the Rybczynski-Hadamard equation:

U=
( ρl − ρv ) gD2
µ . (59)
6 µl
For µv ≪ µl , inequality equation (58) and equation (59) become:

⎛
D < 2⎜
3µl2 ⎞
⎟
1/3
and U =
( ρl − ρ g ) gD2
. (60)

(
⎝ g ρl − ρv ρl ⎠) 12 µl

(2) Intermediate size bubbles, for which 1 < Re < 700. These bubbles have only a slightly
deformed spherical shape, and at Reynolds numbers close to 700, they begin to vibrate and rise
by a spiral motion instead of a straight line. Diameters of these bubbles are in the range:

1/3 1/5
⎛ 3µl2 ⎞ ⎡ 324 µ 2σ ⎤
2⎜ ⎟ ≤ D ≤ 2 ⎢ 2 l3 ⎥ . (61)
(
⎝ g ρl − ρv ρl ⎠) ⎢⎣ g ρl ⎥⎦

Velocity is calculated as:

U=
( ρl − ρv ) gD2
. (62)
36 µl

(3) Large bubbles—fully deformed bubbles with Re > 700. Their terminal velocity does not
depend on size:
1/5
⎡ 4σ 2 g ρ ⎤
U=⎢ v . (63)
2 ⎥
⎢⎣ 30 µl ρl ⎥⎦

Increasing their Re number by increasing size only can lead to bubble instability and bubble
breakup. The critical diameter of the bubbles, a diameter at which breakup of large bubbles begins,
can be calculated by the Levich equation (see ref. 18):

σ σ
Dcr ≈ 3 96 ≈ 4.4589 , (64)
( ) ( )
1/3 1/3
U 2 ρv ρl2 U 2
ρv ρl2

where velocity, U, is calculated by equation (63).

27
 6. TANK THERMODYNAMIC VENTING PRESSURE CONTROL SYSTEMS

The TVS typically includes a Joule-Thomson (J-T) expansion device, a two-phase heat
exchanger, and a mixing pump to destratify and extract thermal energy from the tank without sig-
nificant liquid losses. When tank pressure control cannot be achieved by mixing alone (bulk liquid
becomes saturated at the ullage pressure), a small amount of liquid extracted from the recirculation
flow is passed through a J-T valve where it is expanded to a lower pressure and temperature. The
cold two-phase mixture is then passed through the cold side of the heat exchanger, which extracts
thermal energy from the recirculation flow, and subsequently is vented out.

6.1 Thermodynamic Vent System Pressure Control Logic

The pressure control system (PCS) has two operational modes—mixing with and without
venting. Depending on mission requirements, modes can be combined in three different pressure
control logics:

(1) Mixing mode is released by switching recirculation pump ‘ON’ when ullage pressure
reaches the allowed maximum, Pu ≥ Pu , and switching ‘OFF’ when ullage pressure becomes
max
lower than the allowed minimum, Pu < Pu . After liquid saturation pressure reaches the allowed
min
maximum Pl , it holds liquid saturation pressure in the range of Pl to Pl independently
max min max
by switching the pump and venting ‘ON’ and ‘OFF’ at the same time.

(2) The PCS holds ullage pressure in the range of Pu to Pu with mechanical mixing
min max
(no venting) by switching recirculation pump ‘ON’ and ‘OFF’. After liquid saturation pressure
reaches Pu , hold ullage pressure in the required range by switching the pump and TVS ‘ON’ and
min
‘OFF’ at the same time.

(3) The PCS holds ullage pressure in the range Pu – Pu with mechanical mixing (no
min max
venting) by switching recirculation pump ‘ON’ and ‘OFF’. If ullage pressure exceeds Pu and
max
liquid saturation pressure is higher than Pu , switch pump and the TVS ‘ON’ – ‘OFF’ at the same
min
time. Switching between these cases is made by an input flag in the code.

28
 A flowchart of the subroutine, which releases these logics, is presented in figure 15.

Begin

1 Regime 3

Mixing Yes pu > pumax pu > pumax Yes Mixing plsat > plmin Yes

No No No Mixing &
Venting

No plsat < plmax or No pu < pumin pu < pumin No dPu/dt > 0 plsat < plmin No
dPu/dt>0

plsat < plmin

Yes Yes Yes Yes
Yes
Yes
Mixing & Mixing &
Mixing No Mixing No Mixing Venting No Mixing
No Venting

End

Figure 15. Thermodynamic venting system pressure control logic.

F15_1627

6.2 Spray Bar Thermodynamic Venting System

Spray bar TVS is used for pressure control in long-term storage cryogenic tanks under dif-
ferent acceleration conditions. As depicted in figure 16, it includes a recirculation pump, a TVS
venting valve, a liquid cooling device, a J-T valve, a spray manifold, a heat exchanger, an injection
tube manifold, injection tubes, and an external venting line with a backpressure device (orifice or
critical Venturi).

29
 Injection Pipe
Manifold
External
Venting Line
Injection Pipe

m⋅ Ivent
quit m⋅ mixu
quhex
Heat
Exchanger
Back Pressure
Device
qlit
Spray Bar qmex m⋅ mixl
Manifold
qlhex
m⋅ pump

m⋅ man

qpump

Joule-Thomson TVS Vent Recirculation

Valve Valve Pump

Figure 16. Spray bar TVS schematic.

F16_1627

6.2.1 Heat Exchanger Pressure Loss

The heat exchanger pressure loss model is based on the Lockhart-Martinelli19 two-phase
flow correlations.

To calculate the total pressure change inside the heat exchanger, the momentum conserva-
tion equation is applied. In accordance with figure 17, with constant cross-sectional areas, this
equation is:
dP ⎛ dP ⎞ ⎛ dP ⎞ ⎛ dP ⎞
= −⎜ −⎜ −⎜ ,

dz ⎝ dz ⎠ fr ⎝ dz ⎠ g ⎝ dz ⎟⎠ mom (65)
⎟ ⎟

where

⎛ dP ⎞
⎜⎝ ⎟ = ⎡α ρ + (1− α ) ρl ⎤⎦ g cos ω is the gravitational hydrostatic pressure loss
dz ⎠ g ⎣ v

30
 ω

p
p+d

exi
D
p

dz
e
mn
D
m⋅ man m⋅ Ivent

Figure 17. Spray bar heat exchanger pressure dropF17_1627

calcualtion schematic.

⎛ dP ⎞ d ⎧⎪ 2 ⎡ X 2 (1− X )2 ⎤ ⎫⎪
⎜⎝ ⎟ = ⎨G ⎢ + ⎥ ⎬ is the momentum pressure loss,
dz ⎠ mom dz ⎪⎩ ⎢⎣ ρvα ρv (1− α ) ⎥⎦ ⎪⎭

and

⎛ dP ⎞ is the frictional pressure loss.

⎜⎝ ⎟
dz ⎠ fr

In these equations,

G = two-phase mass flux through cross section of the heat exchanger, kg/(m2s)
X = vapor quality inside heat exchanger
a = Av  /A – void fraction – ratio of the gas flow cross-sectional area to the total
cross-sectional area.

According to reference 19, void fraction is calculated by:

0.36 0.07 −1
⎡ ⎛ ρv ⎞ ⎛ µl ⎞ 0.64 ⎤
⎢ ⎛ 1− X ⎞ ⎥
α = 1+ 0.28 ⎜ ⎟ ⎜µ ⎟ ⎜⎝ ⎟⎠ . (66)
⎢ ρ
⎝ l⎠ ⎝ v⎠ X ⎥
⎣ ⎦

31
 Using Lockhart-Martinelli19 and Martinelli-Nelson20 two-phase multiplier techniques, the
friction pressure loss gradient can be written as:

⎛ dP ⎞ 2 ⎛ dP ⎞ (67)
⎜⎝ ⎟⎠ = Φlo ⎜⎝ ⎟ .
dz fr dz ⎠ lo

In these equations ( dP dz )lo is a pressure gradient that would result if the liquid flowed
alone through the channel with the total mass flux, G.

Multiplier Φlo can be found from another multiplier Φl by equation (68):

Φ2lo = Φ2l (1− X )1.75 . (68)

As defined by Lockhart and Martinelli19, multiplier Φl 2 is:

 C 1 
Φ2l =  1 + + 2  , (69)
 X lm Xlm 

where Xlm is a Lockhart-Martinelli parameter, modified for boiling flow inside the pipe20:

0.571 0.143
⎛ρ ⎞ ⎛ µl ⎞ ⎛ 1− X ⎞
Xl m = ⎜ v ⎟ ⎜µ ⎟ ⎜⎝ ⎟ . (70)
⎝ ρl ⎠ ⎝ v⎠ X ⎠

The value of constant C in equation (69) depends on the liquid and gas flow regimes, which can be
found using Reynolds numbers for two-phase flow:

G(1− X )Dhex GXDhex

Rel = , Rev = . (71)
µl µv

The flow is laminar if liquid or gas Re < 2,000 and turbulent if Re ≥ 2,000. Recommended values
of the constant C are presented in table 1.

Table 1. Recommended values for constant C

in Lockhart-Martinelli parameter.
Liquid Gas C
Turbulent Turbulent 20
Laminar Turbulent 12
Turbulent Laminar 10
Laminar Laminar 5

32
 As suggested by Changhong et al.21, the constant C for narrow annuli depends on the
annuli hydraulic diameter, and can be calculated as:

(
C* = 21 1− e
−0.319Dhex
) . (72)
In equation (72), Dhex is an annuli hydraulic diameter and Dhex = Dexi – Dmne, where Dexi
is the internal diameter of the heat exchanger external pipe and Dmne is the external diameter
of the heat exchanger internal pipe or a spray bar manifold tube (see fig. 17).

With assumption that the two-phase flow in a heat exchanger is homogeneous, the pressure
gradient (dP/dz)lo can be written as reported by Carey22:

⎛ dP ⎞ 2 floG 2
⎜⎝ ⎟ = , (73)
dz ⎠ lo ρl Dhex
−0.25
⎛ GDhex ⎞
where flo = 0.079 ⎜ ⎟ .
⎝ µl ⎠

For the one-phase flow (liquid or gas) inside the heat exchanger, the pressure loss is calcu-
lated by modification of equation (65) without the momentum term, which can be neglected for this
kind of flow:

dP G2
=−f − ρ g cos ω , (74)
dz 2 ρ Dhex

where r is a density of liquid or gas, depending of what kind of fluid is flowing in the calculated
segment, and f is a friction factor, which calculates by the Ghabari-Farshad-Rieke correlation23:

−2.169
⎧⎪ ⎡⎛ ε D ⎞ 1.042 ⎛ 2.731⎞ 0.9152 ⎤ ⎫⎪
f = ⎨1.52 log10 ⎢⎜ hex +⎜
⎥⎬ , (75)
⎢⎣ ⎝ 7.21 ⎟⎠ ⎝ Re ⎟⎠ ⎥⎦ ⎪⎭
⎪⎩

GDhex
where e is the roughness of the internal surface of the spray manifold, m, and Re = is
a Reynolds number, calculated using liquid or gas viscosity. µ

This equation is applicable for a wide range of Reynolds numbers and relative roughness.
In addition, equations (74) and (75) are used for calculating pressure drop inside the spray bar
manifold, but in that case, hydraulic diameter is replaced by the manifold diameter.

6.2.2 Heat Exchanger Heat Transfer

The heat exchanger is a critical element since its function is to ensure optimal heat transfer
between the vented and recirculating fluids. The exchanger design must be capable of rejecting the

33
maximum environmental heat leak rate anticipated as well as simultaneously reducing the liquid
bulk temperature within a reasonable timeframe.

The model of heat exchanger is a multinode finite difference model that simulates two-phase
flow in a quasi-steady-state mode (see fig. 18). It is based on two-phase, fluid-forced convection and
nucleate boiling heat transfer inside narrow annuli, surrounded by two concentric pipe as shown
in figure 18.

Spray
Heat Bar
Exchange

beg
fin

i+1

Tbeg(i+1) = Tfin(i )

Pbeg(i+1) = Pfin(i ) beg

fin
Ebeg(i+1) = Efin(i ) hext hexm
kmex
ktex i
hmex
htex
beg
fin

m⋅ Ivent

m⋅ man

Figure 18. Spray bar heat exchanger transfer calculation schematic.

F18_1627

The total heat transfer coefficient from the spray manifold to the heat exchanger for each
node is calculated by a standard procedure4:

1
hmext = , (76)
1 1 δ
+ + man
hmex hexm kmex

where hmex is the forced convection heat transfer coefficient inside the spray manifold, W/m2 K;
it  calculates by the Dittus-Boelter correlation4:

kl m! man µ Cp
hmex = 0.023 Rel0.8Prl0.4, Rel = , Prl = l l , (77)
Dman Dman µl kl

34
where

kmex – thermal conductivity coefficient of the spray bar manifold pipe metal wall, W/m·K
dman – manifold wall thickness, m
Dman – spray bar manifold pipe internal diameter, m
m! man – spray bar manifold mass flow rate, kg/s
hexm – two-phase flow heat transfer coefficient inside the heat exchanger, W/(m2·K),
calculated by the Chen correlation24, hexm = hconv + hboil.

In this equation, hconv is a heat transfer coefficient for a forced convection calculated
by the Dittus-Boelter equation and modified for two-phase flow inside a narrow annuli:

kl 0.8 0.4
hconv = 0.023 Retp Prl , (78)
Dhex
where

( )
1.25
Retp – Reynolds number modified for two-phase flow, Retp = Rel ⎡⎣ F Xl m ⎤⎦
Xlm – two-phase flow Martinelli-Nelson parameter, calculated by equation (70)
Rel – liquid Reynolds number for a heat exchanger two-phase flow, calculated
by equation (71)
hboil – boiling contribution, calculated by the Forster-Zuber correlation corrected with
a suppression factor S(Retp)22:

⎡ k 0.79Cp0.45 ρ 0.49 ⎤
( ) l
hboil = 0.00122S Retp × ⎢ 0.5 l
⎢⎣ σ µl hlv ρv ⎥⎦
l
0.29 0.24 0.24 ⎥ ⎣ w ( )
⎡T − Tsat Pl ⎤
⎦
0.24
⎣ ( )
⎡ Psat Tw − Pl ⎤
⎦
0.75
. (79)

The suppression factor S(Retp) corrects the fully developed nucleate boiling prediction to
account for the fact that, as the convection effect increases in strength, nucleation is more strongly
suppressed. Both F(Xlm) and S(Retp) functions are calculated by Collier approximations25 as:

⎧1, −1
X lm ≤ 0.1
⎪⎪
( )
F X lm = ⎨ ⎛ 1 ⎞ −1
(80)
⎪2.35 ⎜ 0.213 + , X > 0.1
⎪⎩ ⎝ X lm ⎟⎠ lm

and

( )
S Retp =
1
1+ 2.56 × 10 −6 Retp
1.17
. (81)

35
6.2.3 Spray Injection Tube Model

The spray injection tube is used for delivering the sprayed liquid to the ullage and bulk
liquid. Its model is multinode, which assigns a node to each orifice (fig. 19).

quit
m⋅ i –1
dz
i –1 m⋅iti –1
Pu
m⋅ i
Liquid

Level ⋅
Piti i miti
Pl
⋅ ⋅
mi – miti
qlit
i +1
Dorif

Figure 19. Spray injection tube schematic.

Pressure loss for each node calculates by the following equation:F19_1627

2
dP ⎛ m! ⎞ 1
= − f ⎜ it ⎟ + ρl g cos ω , (82)
dz ⎝ Ait ⎠ 2 ρl Dit
where

f – friction factor, calculated by Ghabari-Farshad-Rieke correlation (eq. (75))

m! it – initial mass flow rate for each node, kg/s
Ait – cross-sectional area of an injection tube, m2
Dit – injection tube internal diameter, m.

36
 m! D
In equation (82), the Reynolds number is Re = it it , with w being the same angle as shown
in figure 17. Ait µl

The mass flow rate through orifices on each level is calculated by the following equation:

( )
1/2
⎡ 2ρ P − P ⎤
m! iti = Aorif ⎢ li t iti t ⎥ , (83)
⎢⎣ Kr ⎥⎦

where

Aorif – cross-sectional area of the orifice, m2

ρlit – density of liquid in the injection tube, kg/m3
Piti – pressure inside injection tube, Pa
Pt – pressure outside the injection line (tank pressure), Pa
Kr – resistance coefficient, which can be found in figure 20. This dependency was
approximated using Idelchik’s Handbook26.

The above parameters are specific for different orifices.

The maximum value of the velocities ratio can be found as:

Vit m! it ρlit Aorif m! Aorif

= = it . (84)
Vorif ρlit Ait m! orif m! orif Ait

At the the first orifice, where the mass flow rate in the injection tube is the highest, it can
be noted that Ait = norif Aorif . In other words, if there is no pressure loss in the injection pipe,
m! it = norif m! orif . Equation (84) then becomes equal to unity.

For other than the first orifice level, the flow rate in the injection tube becomes lower
because of liquid outflowing from the tube. In addition, if there is the pressure loss, the cross-
sectional area of the tube should be higher to increase pressure inside the injection tube, so the
velocities ratio will be less than the unity. In the multipurpose hydrogen test bed (MHTB), for
example, the spray bar injection tube maximum velocities ratio is equal to 0.77. From figure 20, the
resistance coefficient in this case is  equal to 2.76.

37
 2.85

2.8

2.75
Kr

2.7

2.65
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Vit /Vorif

F20_1627
Figure 20. Resistance coefficient versus ratio of the injection tube to orifice velocities.

It is very important to note that, in TankSIM, all pressure loss calculations (heat exchanger,
spray bar manifold, injection tube manifolds, injection tube, and so on) are made for straight wall
segments of lines. Curved and other restriction segments have to be modeled independently in the
form:
2
⎛ m! ⎞ 1
K r = ⎜ li ne ⎟ , (85)
⎝ Ali ne ⎠ 2 ρ ΔPli ne

and then Kr is used in TankSIM for calculations as input.

6.2.4 Recirculation Pump

As depicted in figure 16, the recirculation pump is operated to circulate liquid in the TVS.
Usually, several types of pumps are used for the TVS. They are different by their parameters and
dependencies between used energy, pressure raised, and mass flow rate. To avoid modeling the par-
ticular pumps and to generalize calculations, in TankSIM, only flow rate—pressure rise and flow
rate—energy functions are used in the form of fifth order equations:

ΔPpump = a0 + a1V! + a2V! 2 + a3V! 3 + a4V! 4 + a5V! 5 (86)

and
Wpump = b0 + b1V! + b2V! 2 + b3V! 3 + b4V! 4 + b5V! 5 ,
(87)

where

∆Ppump = pressure rise inside pump, kPa

Wpump = energy used by pump, W
V! = volumetric flow rate through the pump, m3/s.

38
Coefficients a0 – a5 and b0 – b5 are input numbers, which should be found before using TankSIM.

6.2.5 Joule -Thomson Cooling Device

A J-T device is used to cool a small quantity of bulk liquid, which is then used for condition-
ing propellant tanks during the long-duration missions. The modified Lee equation is used for a J-T
modeling. This equation is applied especially for the Visco Jet—multiple orifice devices, for which
the standard equation for a single orifice cannot be used, and is then modified by Pappel27 for two-
phase cryogenic fluids. For the metric system, the equation is:

4.6736 × 10−4
m! = K f
Lohm
( )
ρΔP 1− X out , (88)

where

!
m – mass flow rate through the J-T device, kg/s
Kf – coefficient, depended of the fluid: for hydrogen, is equal to 0.9,
for methane and nitrogen, 1 (refs. 27–29)
Lohm – hydraulic resistance of the J-T devices, given by a vendor30
r = Xin rv + (1 – Xin) rl – inlet two-phase fluid density, kg/m3
Xin – an inlet two-phase fluid quality
∆P – pressure drop across the J-T device, Pa
Xout = (htot – hl)/(hv – hl) – outlet two-phase fluid quality; this value is calculated assuming
isenthalpic fluid expansion through the J-T device
htot – total enthalpy at J-T inlet temperature and pressure, J/kg
hl and hv – enthalpies of the liquid and vapor, respectively; they are
calculated at the outlet temperature and pressure, J/kg.

6.3 Ullage Droplets Heat and Mass Transfer

Droplets of a subcooled liquid appear in ullage during the spray bar mixing or TVS work-
ing. To calculate the heat and mass transfer to ullage droplets from the ullage and ullage wall, the
follow assumptions were made:

• All droplets created by the equivalent orifices have the same size and velocity.
• Horizontal velocity of droplet does not change during its motion in the ullage.

The process of ullage droplets heat and mass transfer includes several stages:

• Jet atomization creates droplets that move towards tank walls.

• Moving droplets warm up to saturation conditions.
• Droplets boiling in the ullage.
• Droplets not boiled off in the ullage impinge ullage wall.
• Droplets warm up to saturation conditions on the wall (if not warmed up in ullage).
• Droplets boiling on the wall and adding mass to the ullage.

39
• Droplets not boiled off on the wall rebounding back to the ullage or adding to the wall liquid
film.
• All droplets remained in the ullage, rebounded from the wall to the ullage, fall down to the liquid,
or continue to interact with the ullage.

6.3.1 Generation of Droplets and Motion Parameters

The velocity of the droplet generated by an orifice can by calculated by the following
equation:

m! orif
vdr = , (89)
ρl corif Aorif

where corif is the discharge coefficient of an orifice and m! orif is the mass flow rate through
the orifice, calculated by equation (83).

The droplet diameter can be found using the Lyshevskiy31 correlation for atomization
of a liquid jet moving through gas:

− 0.266
⎛ρ ⎞
Ddr = 3.01dorif ⎜ u We ⎟ M 0.0733 , (90)
⎝ ρl ⎠

2
ρl dorif vdr µl2
where We =  and M = .
σl ρl σ l dorif

The number of droplets created by each orifice per second can be calculated as:

m! orif 6
n!orif = 3
. (91)
ρl π Ddr

In equations (89) through (91), subscript ‘l ’  refers to the liquid at the injection tube temperature
and pressure.

During traveling in ullage, droplets can reach the wall or fall down to the liquid depending
on acceleration conditions, initial droplet velocity, and distance from the injection tube to the tank
wall as shown in figure 21.

40
 vdr

vg
horif
hmin

Liquid
Level
L

Figure 21. Droplet motion in the ullage.

F21_1627

The shortest time, at which droplets begin to reach the wall, is tmin = L/vdr . The lowest height
of the orifice corresponds to this time and can be found as:

2
g ⎛ L ⎞ (92)
h min = ⎜ ,
2 ⎝ vdr ⎟⎠

where L is a distance between injection tube and tank wall.

For calculating the droplets to ullage heat transfer, average velocities of droplets reaching
and/or not reaching the wall need to be calculated as follows:

1/2
g
vdr = ⎡⎢vdr
2
+ horif ⎤
⎥⎦ , for droplets reaching the wall (93a)
⎣ 2
and

2 1/2
⎡
2 ⎛g L ⎞ ⎤
vdr = ⎢vdr + ⎜ ⎟ ⎥ , for droplets not reaching the wall. (93b)
⎢ ⎝ 2 vdr ⎠ ⎥
⎣ ⎦

41
6.3.2 Ullage Droplets Heat Transfer

The mass conservation equation for the ullage droplets according to figure 5 can be written
as:
dMdr
= ∑ m! i = m! mixu + m! ucndr − m! drboil − m! drl − m! drwl . (94)
dt i
Per figures 5 and 6, the energy conservation equation for ullage droplets is:

Mdr
( )=q
d edr
+ m! ucndr hul + m! mixu hmix − m! drboil hdrl − m! drl hdrl − m! d rwl hdrl
udr
dt
dVdr
−edr ∑ m! i − Pdr . (95)
i
dt

The heat transfer coefficient from ullage to a moving droplet can be calculated by Ranz
and Marshall correlation32 for falling droplets:

kl
hdr =
Ddr ⎢⎣ (
⎡2 + 0.6 Re
dru )
1/2 1/3 ⎤
Prdr , (96)
⎥⎦
where

ρ D v
Redru = u dr dr – Reynolds number, calculated for droplet moving in ullage
µu

µ Cp
Prdr = l l – Prandtl number for droplet liquid.
kl

During the convective transfer stage, which occurs exactly after droplet creation and before
droplet begins boiling, the heat flow rate can be calculated as:

( )
qconv = Adr hdr Tu − Tdr . (97)

On the boiling stage, the heat flow rate is:

( )
qboil = Adr hdr Tu − Tlsat , (98)

where Tlsat is the saturation temperature of the droplet liquid, calculated using ullage pressure.

42
6.3.3 Impinging Droplets—Wall Heat Transfer

For the heat transfer of droplets impinging on the wall, the Pasandideh-Fard et al.33 model
is used (see fig. 22). In this model, the diameter of impinging droplet and heat transfer coefficient
can be found.

Ddr

vdr
Dimp

vg

Figure 22. Droplet impingement model schematic.

F22_1627

The heat transfer coefficient of an impinging droplet is calculated by the following equation:

kl 1/2
himp = Redr Prdr0.4 , (99)
2Ddr

ρv D
where Redr = l dr dr is the droplet Reynolds number with a velocity, normal to the wall.
µl
The maximum diameter of an impinging droplet can be calculated as:

1/2
⎡ ⎤
Wedr + 12
Dimp = Ddr ⎢ ⎥ , (100)
⎣ (
⎢ 3(1− cosΘ ) + 4 We Re1/2
dr dr ) ⎥
⎦

where Θ is the wetting angle, which is close to zero for cryogenic liquids. In that case, equation
(100) will transform to:
1/2
D ⎛ 12 ⎞
Dimp = dr Re1/4
dr ⎜ 1+ ⎟ , (101)
2 ⎝ Wedr ⎠

2
ρl vdr Ddr
where Wedr = is a droplet Weber number with a velocity, normal to the wall.
σl

43
 Having a heat transfer coefficient and diameter of an impinging droplet, the heat flow rate
for one droplet can be calculated as:

2
π Dimp
qimp =
4
( )
himp Tuw − Tdr . (102)

The total heat transfer from a hot surface to an impinging droplet during the time droplet
spreading to its maximum diameter (Dimp in fig. 22) can be found by:

Qtot = q impts , (103)

where ts = 8Ddr / 3vdr.

After impinging, there are two scenarios for the droplets' behavior—depositing on the surface
or rebounding from it. According to Castanet34, to distinguish between liquid rebounding from
the solid surface and depositing on it, the so-called reference temperature, T  *, can be used:

T − Tsat
T* = w , (104)
Tlf − Tsat

where Tlf is a Leidenfrost temperature for impinging liquid.

If 0 <T  *<1 (which implies wall temperature is lower than Leidenfrost), liquid deposits on
the wall. In this case, the contact temperature is used for heat transfer calculations:

ε T +ε T
Twf = w w l l , (105)
εw + εl

where ε l = kl ρlCpl and ε w = kw ρwCp w .

When the reference temperature is higher than unity, after interacting with the wall
and warming up (or boiling), the liquid droplet reflects back to the ullage.

As it was found by Baumeister and Simon35, the Leidenfrost temperature for cryogenic
liquids, with good accuracy can by approximated as:

Tlf = 0.84375Tcrit , (106)

where Tcrit is a critical temperature of liquid.

44
6.3.4 Droplets Heat Transfer After Impinging in the Wall

As it was denoted earlier, after impacting into the wall, liquid remains on the wall and creates
liquid film, or rebound to the ullage as depicted in figure 23.

Vdr1

Vg1
Vdr
Vdr2
Vg2 Vg

Vdr3

Vg3

Vgf
(a) (b)

Figure 23. Droplet behavior after impingement: (a) Spreading and (b) rebounding.
F23_1627

6.3.4.1 Heat Transfer in the Liquid Film. After touching the wall and spreading on it,
the  sprayed droplets form a liquid film covering the wall, and during mixing, continually adds
mass to  it. For a description of the heat transfer in the liquid film, the model of Chun and Seban35
was  applied. By this model, heat transfer coefficients for the nonboiling liquid film are as follows:

−1/3
⎛ ν2 ⎞
hlam = 0.821⎜ l 3 ⎟ Re −0.22
f for a laminar flow (107)
⎝ gkl ⎠
and

−1/3 0.65
⎛ ν2 ⎞ ⎛ νl ⎞
hturb = 0.00381⎜ l 3 ⎟ Re 0.4
f ⎜ ⎟ for a turbulent flow, (108)
⎝ gkl ⎠ ⎝ al ⎠

where νl is a thermal diffusivity of the liquid in a film.

45
To specify a flow regime (laminar or turbulent), Chun and Seban suggested using film Weber
number:

2ρl δ ul2 1 µl2

We = = Re 2f . (109)
σ 16 ρlσδ

If We ≤ 1, there is a laminar film flow; otherwise, if We > 1, there is a turbulent flow. In equations
(107) through (109),

ρ δ u 4M!
Re f = l l = – film Reynolds number
µl µl

ρ 2 gδ 3
M! = – film mass flow rate per unit length of the wall, kg/m∙s
3µ

d – film thickness, m.

46
 7. TANK AXIAL JET THERMODYNAMIC VENTING PRESSURE CONTROL SYSTEMS

7.1 Axial Jet Thermodynamic Venting System

Axial jet TVS is used for ullage pressure control and liquid destratification during long-term
storage by mixing bulk liquid and condensation of vapor at the liquid-ullage interface. The axial jet
system is compact, lightweight, and capable of being integrated with a compact heat exchanger. As
shown in figure 24, it includes a recirculation pump, a TVS venting valve, a liquid cooling device,
a J-T valve, a heat exchanger, a jet nozzle, and an external venting line with a back-pressure device
(orifice or critical Venturi).

m⋅ cond qcond

m⋅ pump
m⋅ jet

Jet Nozzle

Back-Pressure
Device
m⋅ Recirculation
lvent
Pump
Venting Line
Joule-Thomson TVS Vent Valve
Heat Exchanger Valve

Figure 24. Axial jet TVS schematic. F24_1627

7.2 Heat and Mass Transfer at the Ullage-Liquid Interface

There are two different cases of axial jet heat transfer models—high and low nozzle submer-
gences, as depicted in figure 25. It can be shown analytically (see app. D.1) that these two cases can
be divided as follows:

47
 s ⎛ d ⎞
≥ 2.458 ⎜ 1− ⎟ (high nozzle submergence case) (110)
Dint ⎝ Dint ⎠
and

s ⎛ d ⎞
< 2.458 ⎜ 1− ⎟ (low nozzle submergence case) (111)
Dint ⎝ Dint ⎠

where d is a jet nozzle internal diameter.

Dint r
Ts U s (r )
Radial
Turbulent
U m (z) Jet

s s
Liquid
To Axial
α z Turbulent
Jet

d Dt d Dt

(a) (b)

Figure 25. Axial jet TVS submergence cases schematic: (a) High submergence
F25_1627
and (b) low submergence.

For high nozzle submergence, the jet fills the cylinder cross section before it reaches the
interface region and the turbulence near the interface is uniform over the cross section (fig. 25(a)).

According to Carey22 and Sonin et al.37, the average condensation mass flux from the ullage
to the liquid at the interface can be calculated as:

⎛ ⎞
! i 0.418exp −1.2 s , (112)
m! avh = M j ⎜
⎝ Dint ⎟⎠

where m! avh is the average mass flux at the interface at high nozzle submergence, kg/s∙m2.

In the case of low nozzle submergence, the turbulent velocity field is quite different from
the high submergence case, as illustrated in figure 25(b). Turbulent jet spreads from the nozzle, but
remains confined to the axial region. The axial jet impinges on the free surface, and the bulk flow

48
associated with it turns radially outward and forms a radial jet just below the interface. The radial
jet thickens and loses mean velocity (as well as turbulence intensity) as it moves outward, and is
eventually deflected down into the bulk liquid at the wall, after which it mixes with the bulk liquid
below. Mean flow in the bulk liquid outside the jet regions is low; however, some circulation and
mixing does result from entrainment into the two jet regions and the return flow toward the outlet
port. As shown by Thomas38, the condensation mass flux for that case can be written as:

⎡ s ⎤
j 2 (
! × ⎢ 0.26β − 0.18β − 0.077 β
m! avl = M 2 1 )
D
⎥ , (113)
⎣ int ⎦

where

m! avl – average mass flux at the interface at low nozzle submergence, kg/s∙m2

b1,b2 – empirical coefficients which are b1 = 0.34 and b2 = 0.24 by the Brown39 experiments
and b1 = 0.33 and b2 = 0.23 by Lin and Hasan40 numerical solution of the conservation
equations for the k-ε turbulence model.

In equations (112) and (113),

ρ c m! ΔT (1− Ja / 2 )
! = b pb j 0
M , (114)
j
ρ j h fg dDint Prb1/3

where

m! j – mass flow rate of the liquid flowing through the nozzle, kg/c
rj – density of the liquid flowing through the nozzle, kg/m3

Ja ≡ cpb ΔT / h fg – a Jacob number, calculated by bulk liquid parameters
∆T = Tsat – Tint – ullage-liquid interface subcooling, K
∆T0 = Tsat – T0 – a liquid subcooling at the nozzle outflow point, K.

Index “b” means that properties are calculated by bulk liquid temperature and pressure.

Average condensation heat transfer coefficient at the ullage-liquid interface can be

calculated by:
m! av h fg
h= . (115)
ΔT0

From equations (112), (114), and (115) for high nozzle submergence, it can be written:

⎛ s ⎞
havh = H × 0.418 exp ⎜ −1.2 . (116)
⎝ Dint ⎟⎠

49
For low nozzle submergence equations (113), (114), and (115),

⎡ s ⎤
(
havl = H i ⎢ 0.26β2 − 0.18β2 − 0.077 β1 )⎥ . (117)
Dint ⎦
⎣

In both equations (116) and (117), the coefficient H is:

ρb cpb m! j (1− Ja / 2 )
H= . (118)
ρj dDint Prb1/3

7.3 Helicoidally Coiled Tube-in-Shell Heat Exchanger Transfer

The model of heat and mass transfer in helicoidally coiled tube-in-shell heat exchanger
based on two-phase forced convection and nucleate boiling inside the coiled tube, which in this
case, is a venting line, and forced convection inside the annulus—created by internal spacer and
external shell with the coiled tube inside it, as depicted in figure 26(b).

The total heat transfer coefficient from the shell side to the venting liquid inside the coiled
tube calculates by a standard procedure4:

1
hm ext = , (119)
1 1 δ vp
+ +
hvpf hsvp kvp

where

hvpf – heat transfer coefficient from the venting pipe wall to the venting fluid, W/m2·K
hsvp – heat transfer coefficient from the shell side liquid to the venting pipe, W/m2·K
dvp – venting pipe wall thickness, m
kvp – thermal conductivity coefficient of the venting pipe metal wall, W/m·K.

The heat transfer coefficient for inside the tube can be calculated by the Chen correlation24
modified for a coiled tube:

⎛ Dhvp ⎞
hvpft = ⎜ 1+ 3.455
⎝
h (
Dc ⎟⎠ conv
)
+ hboil , (120)

where

2Dd
Dhvp = is an internal hydraulic diameter of coiled pipe, which usually
(D + d ) ⎡⎣3 − 4 − h ⎤⎦

50
has an elliptical shape with the major axis in the vertical direction (see app. D.3), and Dc is a coil
diameter as shown in figure 26(a).

2
⎛ D − d ⎞ where D and d are the internal major and minor axes of the ellipse, respectively.
h=⎜ ,
⎝ D + d ⎟⎠

Dec
Dc
r
d

p
D

H
Internal
Spacer

Coiled Tube
Dis
External
Dic
Shell
Des

(a) (b)

Figure 26. Helicoidally coiled tube-in-shell heat exchanger: (a) Exchanger measurements
F26_1627
and (b) calculation schematic.

The convective heat transfer coefficient inside the coiled pipe is calculated by the Dittus-
Boelter (see ref. 4) equation:

k 4 m! vent µCp
hconv = 0.023 Re 0.8Pr 0.4 , Re = , Pr = . (121)
Dhvp π Dhvp µ k

For a single-phase flow, the Reynolds and Prandtl numbers are calculated using thermody-
namic properties for a particular gas or liquid phase.

For the two-phase flow, the Reynolds number in equation (121) is replaced by:

4 m! vent
( )
1.25
Retp = Rel ⎡⎣ F Xl m ⎤⎦ , Rel = (1− X ) , (122)
π Dhvp µl

51
where X is vapor quality inside the coiled tube and F(Xlm)is the Reynolds number modification
coefficient calculated by equation (81).

In this case, the Prandtl number is calculated using liquid thermodynamic properties. The
boiling heat transfer coefficient is calculated by the Forster-Zuber equation (79) with the two-phase
Reynolds number from equation (122) and the suppression factor from equation (80).The heat
transfer coefficient from the liquid in the shell side to the venting pipe is calculated by the Dittus-
Boelter correlation:

k G D µ Cp
hsvp = 0.023 l Re l0.8Prl0.4 , Rel = shell hs , Prl = l l , (123)
Dhs µl kl
where

4 m! shell
Gshell = is the mass flux in the shell side, kg/(m2s)
π ⎡ Des
⎣
2
(
− Dis2 − 4Dc d + 2δ vp ⎤
⎦)
and

Dis2 − Des
2 ⎡
− Dd + 2δ vp ( D + d ) ⎤ 1+ γ −2
Dhs = 2 ⎣ ⎦ is the hydraulic diameter of the shell
( ) ( )(
2 Dis + Des + D + d + 4δ vp 3 − 4 − hs ) 1+ γ −2
side.

Details of the Gshell and Dhs equations above are provided in appendix D.3.

7.4 Heat Transfer at Low Acceleration Conditions

At low accelerations, the heat transfer for the settled liquid can be modeled by the correla-
tions described earlier with the special calculation of interface and dry wall areas. To do this, the
Surface Evolver program41 can be used. In this program, the interface and dry wall areas as a func-
tion of Bond number, tank geometry, and liquid fill level are calculated. The liquid-ullage interface
shapes calculated by Surface Evolver for hydrogen in the MHTB tank are presented in figure 27.

Bond=0 Bond=0.1 Bond=10 Bond=200 Bond=1,000

Figure 27. Liquid-ullage interface in MHTB tank with 90% fill level.
F27_1627

52
 Interface and dry wall areas in dimensionless form (as a ratio of current value to the value
at flat interface, or at high Bond numbers) are calculated with the Surface Evolver and can be
approximated by the function:

m
⎛ ⎞
1− ⎜ Bo Bo ⎟
⎝ 3⎠
A= a+b m
, (124)
⎛ ⎞
1+ ka ⎜ Bo Bo ⎟
⎝ 3⎠

where coefficients a, b, and ka are calculated from three points on the graph, and coefficient m can
be found by the least squares method or graphically for best matching the Surface Evolver data
(for details, see app. D.5). Figure 28 represents results of Surface Evolver calculations and their
approximations.

3.5
1.2
A2, Bo2
3 1
Interface Area, SE
2.5 Interface Area, SE
Relative Areas

Relative Areas

A3, Bo3 0.8 Interface Area, Approximation

Interface Area, Approximation
A1, Bo1 2 Dry Wall Area, SE
Dry Wall Area, SE 0.6 Dry Wall Area, Approximation A2, Bo2
1.5 Dry Wall Area, Approximation
0.4
1
0.5 0.2
0 0
1×10–5 1×10–3 1×10–1 1×10 1×103 1×10–5 1×10–3 1×10–1 1×10 1×103 1×105 1×107
(a) Bond Number (b) Bond Number

Figure 28. Interface and dry wall areas for hydrogen in MHTB tank: (a) 25% fill level
F28_1627
and (b) 90% fill level.

It has to be noted that for different tank shapes and different fill levels, areas are different,
so for each particular case it needs to find these functions using the Surface Evolver program.

7.4.1 Heat Transfer in Nonmixing Cases

For unsettled liquids during the nonmixing (self-pressurization) period for calculating heat
transfer, the concentric spherical enclosure model is used. It includes an outer spherical tank shell
and ullage sphere inside of it, divided by a liquid layer. The volume of the outer sphere is equal
to  the tank volume, the volume of the internal sphere is equal to the ullage volume, and the volume
of the liquid layer is equal to the bulk liquid volume (fig. 3). Applying the Teertstra et al.42 method,
the heat transfer to the ullage is calculated using Au as a linear scale and conduction shape
factor, S*:
(
qmg = k Au S * Tw − Tu . (125) )

53
In the case of ‘sphere in sphere,’ the shape factor is:

1/3 Au
⎛ 4π ⎞
S* = ⎜ ⎟ + 2 π . (126)
⎝ 3 ⎠
( 1/3
Vtot −Vu1/3 )
Heat flow rate from the tank wall to the ullage is:
⎡ 1/3 Au ⎤
⎛ 4π ⎞
( ⎢
qmg = kl Au Tw − Tu ⎜ ⎟ ) ⎥
+ 2 π . (127)
⎢⎝ 3 ⎠
⎣
1/3
Vtot (
−Vu1/3 )
⎥
⎦

7.4.2 Heat Transfer During Mixing

To calculate heat transfer during mixing, a homogeneous fluid model was applied. In that
model, all fluid inside the tank is assumed to be a homogeneous mixture of vapor, noncondens-
able gas (if any) and liquid. Thermodynamic and physical properties (density, viscosity, conduction
coefficient, specific heat capacity, enthalpy, and other) are calculated as mass average values:

C m + Cv mv + Cl ml
C mix = he he . (128)
mhe + mv + ml

54
 8. COMPUTER MODEL DESCRIPTION

8.1 TankSIM Program Base

The code of the TankSIM program is based on FORTRAN 90 language and is released on
the Intel Visual FORTRAN compiler for Windows, but it can be used at any platforms with minor
adjustments just in external files operations.

The program is structured and consists of the main program and 75 subroutines. All
subroutines can be divided into groups as follows:

• Service and error handling.

• Pressure and mission phases control.
• Initialization and input-output operations.
• Heat and mass transfer in the tank.
• Calculation by NIST procedures or tables interpolation to obtain fluids thermodynamic properties.
• Ullage venting and pressurization.
• Spray bar and axial jet TVS.
• Calculation of tank geometry.

8.2 Program Possibilities

The program TankSIM is created for predicting cryogenic tank temperatures and pressures
under different external conditions. It is possible to use:

• Cylindrical tanks with different upper and lower domes geometries (flat, elliptical, and spherical)
and different lengths of the cylindrical part.
• Different cryogenic fluids and incondensable pressurization gases.
• Different materials for a tank wall (properties of material are input parameters).
• Unbounded number of six mission phase combinations, which are:
– Regular (‘Slf’)—self-pressurization with ullage venting, spray bar or axial jet pressure control
systems.
– Engine and feed system chilldown (‘Chl’).
– Noncondensable gas (‘Hep’) prepressurization.
– Autogenous (‘Vap’) prepressurization.
– Engines firing (‘Fir’) with noncondensable gas or autogenous pressurization.
– Vapor liquefaction (‘Lqf’).
• Different types of back-pressure devices—critical orifice or Venturi and pressure regulation.
• Three different types of pressure control logic.

55
 8.3 Program Error Handling

Program error handling and control system includes two subroutines and more than 75 code
segments. There are two types of code segments—custom error handling and these, which use the
RefProp error procedures. When an error happens, the program pauses and an error message appears
on the interactive screen. This message includes the name of a code block (main program or subrou-
tine), serial number of error handling segment, time that the error happened (calculation time, not
a  system time), and error description. In addition, it includes four possibilities to continue calcula-
tions: just continue; continue without option massage (error is printed on the screen, but program
does not pause); continue without error message (does not print any messages); and stop (stop the
program and error message is printed in the output file).

User needs to choose and input one option. Figure 29 is an example of an error message
that appears on the screen.

Figure 29. Error message example.

It should be noted that the option chosen by the user (excluding ‘Stop’) will be applied only
to a particular serial number. Each first appearance of the error handling segment will give the full
set of options. See figure 30 for a typical custom error handling code segment.

Figure 30. Custom error handling code segment.

56
 In this code segment, subname is the code block name (main program or subroutine), reason
is a custom error description, ier is a unique serial number of the error handling code segment, and
Error_Action is an error handling subroutine.

The typical error handling segment for the standard procedure is different from the custom
sequence. In figure 31, ierr is a RefProp error number and herr is a RefProp error description.

Figure 31. Standard error handling code segment.

8.4 TankSIM Functionality

8.4.1 Program Flowcharts

Figures 32–34 represent working flowcharts for the regular mission phase, which include the
self-pressurization and pressure control system for the cost time period, initialization subroutine,
and heat exchanger transfer subroutine for the spray bar TVS system, respectively.

57
 Care other than Begin
regular phases
No

Regular Mission Phase

Yes Initialization
Phase Control

regime=self No pui>max Yes figpu=1

Yes
Time No
Cycle

figpu=0
Printing
Option

Pressure
Control

regime=self or regime=jetm or
No regime=barm or No
(figpu=0 & jeth
barh
figpl=0)

Yes Yes
Yes
Spray Bar Axial Jet
Calculations Calculations
Low G Heat
Yes Low G case
Transfer

No
Update
Heat & Mass
Print Output Files Parameters for
Transfer Rates Next Time Step

No Final Time

Yes
End

Figure 32. Program flowchart for regular mission phase.

F32_1627

58
 Read Read
Enter
Address File Input Files

thei > 0 Yes phei > 0

No Yes No

Pressurization Bond
Calculation

Helium Ullage-Liquid
Parameter Interface Area
Calculations Calculation

Homogenous
Fluid Calculation No Bond> 0.1

Yes

Spray Bar regime= barm or regime= jetm or

Yes No
Initialization barh jeth
Yes

No

Initial Axial Jet

Gesses Initialization

Initial Energy Time Integration

Calculation of Variables

Return

F33_1627
Figure 33. Initialization subroutine flowchart.

59
 Enter

Initial J-T Valve

Assignment Calculations

Heat Exchanger
Calculation Loop

Non boiling Boiling

No Xei<Xe<1 Yes
nb=1 nb=2

Heat Exchanger
Heat

Heat Transfer Rate Heat Transfer Rate Heat Transfer Rate

From Manifold to From Tank to Inside
Heat Exchanger Heat Exchanger Heat Exchanger

Heat Exchanger Heat Exchanger

Pressure Loss Pressure

Back Pressure
No Loop End Yes
Device

Yes delm>0.0001

No

Return

Figure 34. Heat exchanger transfer subroutine flowchart.

F34_1627

60
8.4.2 File System

The TankSIM file system includes program files, input files, service files, and several output
files. The file structure and file names are hardcoded, excluding output files where names include
date, time of creation, and short description of the case. The overall structure of the TankSIM file
system is represented in figure 35. ‘TankSIM Main Folder’ can have any name and address, because
all internal folders have relative addresses.

Debug

TankSIM
Fluids

TankSIM_InOut
TankSIM
Main Folder

TankSIM_Serv

Figure 35. TankSIM file structure.F35_1627

Folder ‘TankSIM’ contains the program source file, TankSIM FORTRAN project file, and
several files for RefProp subroutines. In addition, the folder includes an address input file and two
folders: ‘Debug,’ which is a compiler service folder (contains executable TankSIM file), and ‘Fluids,’
which contains the RefProp fluids description.

Folder ‘TankSIM_InOut’ contains working fluid and noncondensable gas property files,
initial conditions input file, tank, spray bar, or axial jet measurements and properties, and mission
profile input files. It also can include three optional files with external venting, filling, and pressur-
ization line measurements, and properties for additional heat leak calculations. All input files are
written in ‘text’ format.

Program output files are also included in the ‘TankSIM_InOut’ folder. The main output file,
which includes input and calculated information in table form, is made in ‘comma separated values’
(CSV) format. That format allows an easy conversion to an Excel file to use all its possibilities for
data analysis. If needed, the user can print several additional output files with more detailed infor-
mation in ‘text’ format. Folder ‘TankSIM_Serv’ is used to write and read intermediate information
during working with input files. All service files are automatically deleted after using and the folder
remains empty, although it should still be there.

61
 There is the possibility to use TankSIM executable file. In this case, the structure is the same
as in figure 35, but the folder ‘TankSIM’ contains only a ‘Fluids’ subfolder and two files, executable
and address.

A detailed description is provided in appendix E.

62
 9. PROGRAM VALIDATION AND MULTIPHASE MISSION EXAMPLE

The main objective of this validation effort is to verify the TankSIM code using MHTB
hydrogen tests and compare program prediction and test data and to demonstrate the program
ability for modeling multiphase flight missions.

9.1 Multipurpose Hydrogen Test Bed Data Validation

In this section, TankSIM predictions are compared to MHTB test data. Validation was
made for three different tank fill levels, 25%, 50%, and 90%. For calculations and analysis, two
types of external heat leaks were used—heat leak estimated from boiloff tests and different heat
leaks for self-pressurization and TVS operation, which were found for the best matching test data.

9.1.1 Hydrogen 25% Fill Level Test, File p263968kl, 1996

Results of calculations with constant values for heat leaks and test data for 25% fill level
are presented in figures 36–40.

140
138
136
134
Pressure (kPa)

132
130
128
126
Pu, 25% Fill
124
122 Pu, Calc
120
0 20,000 40,000 60,000 80,000 100,000 120,000 140,000 160,000 180,000 200,000
Time (s)

Figure 36. Ullage pressure, 25% fill level.

F36_1627

63
 140
138
136
134
Pressure (kPa)

132
130
128
126
Plsat, 25% Fill
124
122 Plsat, Calc
120
0 20,000 40,000 60,000 80,000 100,000 120,000 140,000 160,000 180,000 200,000
Time (s)

Figure 37. Liquid saturation pressure, 25% fill level. F37_1627

28

27

26
Temperature (K)

25

24

23
Tu, 25% Fill
22 Tu, Calc
21
0 20,000 40,000 60,000 80,000 100,000 120,000 140,000 160,000 180,000 200,000
Time (s)

Figure 38. Ullage temperature, 25% fill level. F38_1627

21.5

21.4
Pressure (kPa)

21.3

21.2

21.1
Tl, 25% Fill
21
Tl, Calc
20.9
0 20,000 40,000 60,000 80,000 100,000 120,000 140,000 160,000 180,000 200,000
Time (s)

Figure 39. Liquid temperature, 25% fill level. F39_1627

64
 28

27

26
Temperature (K)

25

24

23
Tuw, 25% Fill
22 Tuw, Calc
21
0 20,000 40,000 60,000 80,000 100,000 120,000 140,000 160,000 180,000 200,000
Time (s)

Figure 40. Ullage wall temperature, 25% fill level. F40_1627

Comparison of predicted ullage pressure and liquid saturation pressure with test data are
shown in figures 36 and 37. The ullage and liquid saturation pressure rising rates predicted by
TankSIM during the tank lockup are in a good agreement with that of measured data. The ullage
pressure and liquid saturation pressure drop rates during TVS operation also are in good agree-
ment with those of the test data. However, ullage pressure predicted by TankSIM is slower than
test data during the TVS operation. The number of cycles predicted by TankSIM and that of the
test data are 11.5 and 18, respectively. The average cycling duration for prediction and test data
are  23,050 and 14,821 s, respectively. The predicted liquid saturation pressure range and that of test
data are 120.2 to 135.8 kPa and 120.2 to 134 kPa, respectively.

Figures 38 and 39 show a comparison of present predicted ullage and liquid temperatures
with that of measured data. The predicted ullage temperature is within the 23.3 to 27.7 K range,
while measured values are within the 21.3 to 26.2 K range. Predicted liquid temperature is in good
agreement with that of the test data during the self-pressurization regime but the agreement is
slightly different during the TVS operation. Predicted liquid temperature changes are within the
21 to 21.4 K range, while the measured values are within the 21 to 21.35 K range. From figure 40,
it  can be seen that the ullage predicted wall and test data temperatures are in very close agreement.

Figures 41 to 45 show calculating results with the best matching heat leaks, which are 10 W
(ullage) and 10.1 W (liquid) for initial self-pressurization and 14.5 W (ullage) and 15 W (liquid)
during the TVS operation period. By comparing predicted values with the test data one can
conclude good agreement with the test data during the TVS operation.

65
 140

135

130
Pressure (kPa)

125

120
Pu, 50% Fill
115
Pu, Calc
110
0 10,000 20,000 30,000 40,000 50,000 60,000 70,000 80,000 90,000
Time (s)

Figure 41. Ullage pressure, 50% fill level. F41_1627

140

135

130
Pressure (kPa)

125

120
Plsat, 50% Fill
115
Plsat, Calc
110
0 10,000 20,000 30,000 40,000 50,000 60,000 70,000 80,000 90,000
Time (s)

Figure 42. Liquid saturation pressure, 50% fill level. F42_1627

30
29
28
27
Temperature (K)

26
25
24
23
Tu, 50% Fill
22
21
Tu, Calc
20
0 10,000 20,000 30,000 40,000 50,000 60,000 70,000 80,000 90,000
Time (s)

Figure 43. Ullage temperature, 50% fill level. F43_1627

66
 21.4
21.3
21.2
Temperature (K)

21.1
21
20.9
20.8 Tl, 50% Fill
20.7 Tl, Calc
20.6
0 10,000 20,000 30,000 40,000 50,000 60,000 70,000 80,000 90,000
Time (s)

Figure 44. Liquid temperature, 50% fill level. F44_1627

30
29
28
27
Temperatures (K)

26
25
24
23
22 Tuw, 50% Fill
21 Tw, Calc
20
0 10,000 20,000 30,000 40,000 50,000 60,000 70,000 80,000 90,000
Time (s)

Figure 45. Ullage wall temperature, 50% fill level. F45_1627

9.1.2 Hydrogen 50% Fill Level Test, File p263981t, 1998

Figure 41 presents a comparison of ullage pressure predicted by TankSIM and that of

test data. Pressure rise predicted by TankSIM during the tank lockup increases faster than that
of measured data. However, ullage pressure predicted by TankSIM is somewhat slower than test
data during TVS operation. The ullage pressure drop rate predicted by TankSIM is in a reasonable
agreement with that of test data. The number of cycles predicted by TankSIM and that of test data
are 9 and 10, respectively, which are in close agreement. The average cycling duration for prediction
and test data are 6,770 and 6,053 s, respectively.

TankSIM predicted liquid saturation pressure is in a good agreement with that of test data
as shown in figure 42. The predicted liquid saturation pressure range and that of test data are 112
to 130.6 kPa and 112 to 133.1 kPa, respectively.

67
 Comparison between TankSIM predicted ullage temperature with that of test data is shown
in figure 43. The TankSIM predicted ullage temperature range is 29.5 to 23.8 K, while measured
ullage temperature is 27.9 to 22.2 K.

Figure 44 compares predicted liquid temperature with that of measured data. The predicted
liquid temperature is 20.7 to 21.26 K, while the measure liquid temperature range is 20.7 to 21.35 K.

A comparison between the predicted ullage wall temperature and that of test data is pro-
vided in figure 45. The predicted TankSIM ullage wall temperature range is 24.3 to 29.5 K, while
the measured ullage wall temperature range is 21.6 to 27.1 K.

9.1.3 Hydrogen 90% Fill Level Test, File p263981d, 1998

Figure 46 presents a comparison of ullage pressure predicted by TankSIM and that of test
data. Pressure rise predicted by the tank during the tank lockup is faster than that of measured
data. However, ullage pressure rate predicted by TankSIM is somewhat slower than test data
during TVS operation. The ullage pressure drop rate predicted by TankSIM is in a reasonable
agreement with that of test data. The number of cycles predicted by TankSIM and that of test data
are 8 and 14, respectively. The average cycling durations for prediction and test data are 3,800 and
2,150 s, respectively. TankSIM predicted liquid saturation pressure is in good agreement with that
of test data as shown in figure 47. The predicted liquid saturation pressure range and that of test
data are 110.4 to 117.4 kPa and 110.4 to 121.1 kPa, respectively.

140

135

130
Pressure (kPa)

125

120
Pu, 90% Fill
115
Pu, Calc
110
0 5,000 10,000 15,000 20,000 25,000 30,000 35,000 40,000 45,000
Time (s)

Figure 46. Ullage pressure, 90% fill level. F46_1627

68
 140

135

130
Pressure (kPa)

125

120
Plsat, 90% Fill
115
Plsat, Calc
110
0 5,000 10,000 15,000 20,000 25,000 30,000 35,000 40,000 45,000
Time (s)

Figure 47. Liquid saturation pressure, 90% fill level. F47_1627

Comparison between TankSIM predicted ullage temperature with that of test data is shown
in figure 48. The TankSIM predicted ullage temperature range is 20.9 to 25 K, while the measured
ullage temperature is 20.9 to 24.2 K.

26
Tu, 90% Fill
25
Tu, Calc
24
Temperature (K)

23

22

21

20
0 5,000 10,000 15,000 20,000 25,000 30,000 35,000 40,000 45,000
Time (s)

Figure 48. Ullage temperature, 90% fill level. F48_1627

69
 Figure 49 presents predicted liquid temperature with that of measured data. The predicted
liquid temperature is 20.7 to 21 K while the measured liquid temperature range is 20.7 to 20.9 K.

21
20.95
20.9
Temperature (K)

20.85
20.8
20.75
20.7 Tl, 90% Fill
20.65 Tl, Calc
20.6
0 5,000 10,000 15,000 20,000 25,000 30,000 35,000 40,000 45,000
Time (s)

Figure 49. Liquid temperature, 90% fill level. F49_1627

9.1.4 Hydrogen 25% Fill Level Test, File p263968kl, 1996, Adjusted Heat Loads

The example of adjusted heat loads is provided in this section. Three sets of heat loads
are shown in table 2. The first set is estimated from the boiloff test, which is usually made before
tank lockup. The second set is the best matched heat leaks during self-pressurization. As shown
in  table  2, these values are close to the estimated values. The third set is values for the best matched
during the TVS operation.

Table 2. Heat leaks in 25% fill level test, adjusted heat loads.
Estimated Best Matching Best Matching
Fill Level 25% From Boiloff Lock-up Values TVS Values
Ullage 9.64 10.0 14.5
Liquid 9.06 10.1 15.0
Total 18.70 20.1 29.5

Results of modeling compared to test data are presented in figures 50–54. By comparing
figures 36–40 with figures 50–54, one can conclude that the calculated parameters for the adjusted
heat leaks lead to much closer agreements with the test data.

70
 140
138
136
134
Pressure (kPa)

132
130
128
126
126 Pu, 25% Fill
122 Pu, Calc
120
0 20,000 40,000 60,000 80,000 100,000 120,000 140,000 160,000 180,000 200,000
Time (s)

Figure 50. Ullage pressure, adjusted heat loads, 25% fill level.
F50_1627

139
137
135
133
Pressure (kPa)

131
129
127
125
123 Plsat, 25% Fill
121 Plsat, Calc
119
0 20,000 40,000 60,000 80,000 100,000 120,000 140,000 160,000 180,000 200,000

Time (s)

Figure 51. Liquid saturation pressure, adjusted heat loads, 25% fill level.

F51_1627

71
 30
29
28
27
Temperature (K)

26
25
24
23 Tu, 25% Fill
22
Tu, Calc
21
20
0 20,000 40,000 60,000 80,000 100,000 120,000 140,000 160,000 180,000 200,000

Time (s)

Figure 52. Ullage temperature, adjusted heat loads, 25% fill level.
F52_1627

21.4
21.35
21.3
21.25
Temperature (K)

21.2
21.15
21.1
21.05
Tl, 25% Fill
21
Tl, Calc
20.95
20.9
0 20,000 40,000 60,000 80,000 100,000 120,000 140,000 160,000 180,000 200,000

Time (s)

Figure 53. Liquid temperature, adjusted heat loads, 25% fill level.
F53_1627

28

27

26
Temperature (K)

25

24

23
Tuw, 25% Fill
22 Tuw, Calc

21
0 20,000 40,000 60,000 80,000 100,000 120,000 140,000 160,000 180,000 200,000

Time (s)

Figure 54. Ullage wall temperature, adjusted heat loads, 25% fill level.
F54_1627

72
 9.2 Multiphase Mission Example With Ullage Venting

TankSIM allows combining the number of mission phases, such as:

• Self-pressurization (or tank lockup) with different pressure control systems.

• Autogenous or noncondensable gas prepressurization.
• Feed lines chilldown.
• Engine firing with tank pressurization.
• Propellant liquefaction.

Each of these phases could have its own duration, heat leaks, acceleration, mass flow rates,
and other parameters. Usually, each mission is a combination of the first four phases. Figures 55–57
present the model predictions for results for ullage and liquid saturation pressures as well as ullage
and liquid temperatures. The model simulates a mission, which consist of the following nine phases:

(1) First self-pressurization with cycling ullage venting.

(2) Second self-pressurization with different heat leaks and acceleration.
(3) First prepressurization with noncondensable gas.
(4) First feed lines chilldown.
(5) First engine firing.
(6) Third self-pressurization with cycling ullage venting.
(7) Second prepressurization with noncondensable gas.
(8) Second feed lines chilldown.
(9) Second engine firing.

225
4 8
Pu
205
Plsat
5
185 3 7
Pressure (kPa)

165 9
1 2 6
145

125

105
0 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 4,500 5,000 5,500 6,000
Time (s)

F55_1627
Figure 55. Ullage and liquid saturation pressures (cycling ullage venting).

73
 120
Tu
100 Tl
Temperature (K)

80

60

40

20
0 500 1,000 1,500 2,000 2,500 3,000 3,500 4,000 4,500 5,000 5,500 6,000

Time (s)
F56_1627
Figure 56. Ullage and liquid temperatures (cycling ullage venting).

45
Tu 21.39
40 Tl 21.37

Liquid Temperature (K)

Ullage Temperature (K)

21.35
35
21.33
30 21.31
21.29
25
21.27
20 21.25
2,800 3,000 3,200 3,400 3,600 3,800 4,000 4,200 4,400 4,600 4,800

Time (s)
F57_1627
Figure 57. Ullage and liquid temperatures (cycling ullage venting, large scale).

74
 10. SUMMARY AND FUTURE PROGRAM IMPROVEMENT

10.1 Summary

A full and detailed description of the TankSIM program is presented in this TM. At this
time, TankSIM is one of the most robust non-computational fluid dynamics (CFD) programs for
predicting heat and mass transfer inside cryogenic propellant tanks with different external and
internal conditions. Program features include the following:

• Use of different tank shapes, working fluids, and noncondensable gases.

• Use of different types of pressure control systems—ullage venting in permanent and cycling
regimes, spray bar and axial jet TVS, and their combinations in different mission phases.

• Simulation of full mission profiles by combining different mission phases in needed sequences
with the ability to add new phases as required.

• Computations with relatively high precision and acceptable run time.

10.2 Future TankSIM Improvements

For future improvement of TankSIM, the following tasks are recommended:

• Develop new correlation equations for tank wall-to-ullage, tank wall-to‐liquid, and liquid-to-
ullage interface heat transfer under normal and low gravity, dependent on tank geometry, endcaps
ratio, barrel section height, and fill level. This can be done using CFD software and regression
analysis.

• Improve calculation of the ullage-to-liquid and tank wall-to-ullage interface areas, dependent
on Bond number, fill level, and shape of the tank. This can be done using ‘Surface Evolver’
software and regression analysis.

• Develop calculation of the concentration-dependent coefficient for condensation mass flow rate
when noncondensable gas is present in the tank.

• Replace nonsettled liquid heat transfer calculations during mixing from homogeneous fluid model to
gas ullage and set of liquid spheres; diameters can be found from the mixer jet spray Weber number.

• Create a graphical user interface with Visual Studio tools.

• Continue to validate the program paying special attention to preparation of test data and develop
an algorithm for estimating experimental heat leaks. Find the sources of discrepancies between
predictions and test data at TVS operation and improve the program to reduce differences.

75
 APPENDIX A—CONSERVATION EQUATIONS FOR MULTINODE ANALYSIS
OF CRYOGENIC STORAGE TANKS

A.1 Leibnitz’s Theorem

The following theorem is useful in expressing the mass and energy conservation equations
in integral form. From Panton43, the Leibnitz’s Theorem for the time derivative of a volume
integral of a scalar, f, where the control volume, V, changes with time, t, may be written as:

⎧ ⎫ !
d ⎪ ⎪ ∂f

dt ⎪ ∫
⎨ f dV ⎬= ∫ ( )
dV + ∫ f VS i n̂S dS . (129)
⎩V (t) ⎪⎭ V (t) ∂t S(t)

!
In equation (129), the scalar function,
! f, may be a function of time and space, f = f ( t, r );
!
S is the surface bounding the volume, V ; VS is the surface velocity, and nS is the outward unit
normal to surface S. The above equation may also be rearranged to express the volume integral
of  ∂f / ∂t as

⎧ ⎫ !
∂f d ⎪ ⎪
∫ ∂t dV =
dt ⎪ ∫
⎨ f dV ⎬ ∫
− (f V )
S i n̂S dS . (130)
V (t) ⎩V (t) ⎪⎭ S(t)

By applying equation (129) to the scalar function f = 1, a useful relationship between the rate
of volume change and surface velocity is found:

⎧ ⎫ !
dV d ⎪ ⎪

dt dt ⎪ ∫
= ⎨ f dV ( )
⎬ − ∫ VS i n̂S dS . (131)
⎩V (t) ⎪⎭ S(t)

If the control volume, V, is constant (i.e., V is not changing with time), then

d ⎧⎪ ⎫⎪ ∂f dV
dt ⎪ ∫
⎨ f dV ⎬= ∫ dV , = 0 . (132)
⎩V ⎪⎭ V ∂t dt

A.2. Conservation of Mass

The conservation of mass for no creation or destruction sources within the volume, V,
follow Panton40 and may be written in integral form as

76
 ⎧ ⎫ ! !
d ⎪ ⎪
⎨ ∫
dt ⎪
ρ dV ( )
⎬ + ∫ ρ V − VS i n̂S dS = 0 . (133)
⎩V (t) ⎪⎭ S(t)

This also follows directly from integrating the differential form of conservation of mass:

∂ρ !

∂t
( )
+ ∇ i ρV = 0 . (134)

!
( )
Over the control volume, V, while applying the Gauss’s theorem to the ∇ i ρV term,
! !
∫ ∇ i( )
ρV dV = ∫ S i n̂S dS , (135)
ρ V
V S

∂ρ
and applying equation (130) to the term ∫ dV where the scalar function, f, in equation (130)
corresponds to the fluid density, r. ∂t
V (t)

Using the following definitions for the terms of equation (132),

! !
m= ∫ ρ dV ; m! k = ∫ ρ (V −VS ) i n̂S dS; S(t) = ∑ Sk (t) , (136)
V (t) S(t) k

where
m – fluid mass within the volume
V – volume
m! k – mass flow rate out of the volume through the surface element
S – total surface
Sk – surface element.

Using the above definitions, the integral form for the conservation of mass previously given by
equation (132) as
d
{m} + ∑ m! k = 0 . (137)
dt k

1
V ∫
A volume-average density also can be defined as m = ρV, where ρ = ρ dV , and the conserva-
tion of mass written as V (t)

d
{ ρV} + ∑ m! k = 0 . (138)
dt k

77
 A.3 Conservation of Thermal Energy

The conservation of thermal energy in differential form written in terms of specific internal
energy, e, is

∂( ρe) ! !

∂t
( )
!
(
+ ∇i ρeV = − ( ∇i q ) − P ∇i V + φ , (139) )
or written in terms of specific enthalpy, h, is

∂( ρ h) ! ! ∂P !

∂t
( )
+ ∇i ρ hV = − ( ∇i q ) +
∂t
(
+ V i ∇P + φ , (140) )
where φ is the viscous dissipation function (which is always positive). Integrating over a control
volume, V(t), and applying Leibnitz’s rule given by equation (130) to the volumes integrals,
∂  (ρe) /∂t and ∂  (ρh) /∂t, the above two thermal energy equations become

⎧ ⎫ ! ! !
d ⎪ ⎪

dt ⎪
( )
!
⎨ ∫ ρe dV ⎬ + ∫ ρe V − VS i n̂S dS = − ∫ q i n̂S dS − ∫ P ∇ iV dV + ∫ φ dV (141) ( )
⎩V (t) ⎪⎭ S(t) S(t) V (t) V (t)

and

⎧ ⎫ ! ! !
d ⎪ ⎪ ⎡ ∂P
⎨ ∫
dt ⎪
ρ h dV ( ) !
⎬ + ∫ ρ h V − VS i n̂S dS = − ∫ q i n̂S dS − ∫ ⎢
⎣ ∂t
⎤
(
+ V i ∇P ⎥ dV + ∫ φ dV . (142)
⎦
)
⎩V (t) ⎪⎭ S(t) S(t) V (t) V (t)

For a constant! control volume, V (i.e., V is not changing with time), the control volume
surface velocity is VS = 0.

Volumes integrals in equations (140) and (141) can be simplified by ignoring the spatial
!
variation of pressure within the control volume, P ( r,t ) ≈ P(t), which is a reasonable approximation
for low speed flows.

For the thermal energy equation in terms of specific internal energy, the above described
approximation for P is used to simplify the following integral involving pressure:
! ! !
∫ ( )
P ∇ iV dV ≈ P ∫ (∇ iV ) dV = P ∫ V i n̂S dS
V (t) V (t) S(t)

! ! ! ! !
∫ (V −VS ) i n̂S dS + P ∫ ∫ (V −VS ) i n̂S dS + P
dV
=P VS i n̂S dS = P , (143)
dt
S(t) S(t) S(t)

78
 !
( )
where equation (131) was used to relate the surface integral of VS i n̂S to the rate of volume change
!
dV/dt. For the thermal energy equation in terms of specific enthalpy, h, P ( r,t ) ≈ P(t) is used to
simplify the following integral involving pressure:

⎡ ∂P !

∫ ⎢⎣ ∂t + V( i ∇P
⎤
⎥⎦ dV )
= ∫
dP
dt
dV =
dP
dt ∫ dV =
dP (144)
dt
V.
V (t) V (t) V (t)

!
Using the simplifications given above which are valid when P ( r,t ) ≈ P(t), the integral form
of the thermal energy equations in terms of specific internal energy, e, is

⎧ ⎫ ! !
d ⎪ ⎪
( ) !
⎨ ∫ ρe dV ⎬ + ∫ ρ h V − VS i n̂S dS = − ∫ q i n̂S dS − P
dt ⎪
dV
dt
+ ∫ φ dV , (145)
⎩V (t) ⎪
⎭ S(t) S(t) V (t)

and the integral form of the thermal energy equations in terms of specific enthalpy, h, is

⎧ ⎫ ! !
d ⎪ ⎪
( ) !
⎨ ∫ ρ h dV ⎬ + ∫ ρ h V − VS i n̂S dS = − ∫ q i n̂S dS +V
dt ⎪
dP
dt
+ ∫ φ dV . (146)
⎩V (t) ⎪
⎭ S(t) S(t) V (t)

In equation (143), for the internal energy, the relation ρe + P = ρh was employed.
! For the
(
case in which the control volume is constant in size and shape, dV / dt = 0 and VS = 0 , )
the above thermal energy equations simplify to

d ⎧⎪ ⎫⎪ ! !
⎨ ∫ ρe dV ⎬ + ∫ ρ hV i n̂S dS = − ∫ q i n̂S dS + ∫ φ dV (147)
dt ⎪ ⎪⎭ S
⎩V S V

and

d ⎧⎪ ⎫⎪ ! ! dP
⎨∫ ρ h dV ⎬ ∫
+ ρ h V i n̂S dS = − ∫ q i n̂S dS +V + ∫ φ dV . (148)
dt ⎪ ⎪⎭ S dt
⎩V S V

Defining mass-averaged internal energy, e , enthalpy, h , and a total fluid mass, m, within
the control volume as

1 1
e=
m ∫ ρ edV , h =
m ∫ ρ h dV , m = ∫ ρ dV , (149)
V (t) V (t) V (t)

79
The equation can be written as
⎧ ⎫ ⎧ ⎫
d ⎪ ⎪ d d ⎪ ⎪ d

⎨ ∫
dt ⎪
ρ edV ⎬ = { me } ; ⎨ ∫ ρ hdV ⎬= { }
mh . (150)
⎩V (t) ⎪⎭ dt dt ⎪
⎩V (t) ⎪⎭ dt

Using mass conservation equation (137), the time derivative of mass, the equation can be
written as

d dm de de
{me } = e +m =m − e ∑ m! k (151)
dt dt dt dt k
and

d dm de dh

dt
{ }
mh = h
dt
+m
dt
=m
dt
− h ∑ m! k . (152)
k

After replacing integrals in equations (147) and (148) by the definition from equation (150),
the thermal energy equations in integral form become

!
m
de
dt k
( )
− ∑ m! k e + ∫ ρ h V −VS i n̂S dS = − ∫ q i n̂S dS − P
dV
dt
+ ∫ φ dV (153)
S(t) S(t) V (t)

and
!
m
dh
dt
( )
− ∑ m! k h + ∫ ρ h V −VS i n̂S dS = − ∫ q i n̂S dS +V
dP
dt
+ ∫ φ dV . (154)
k S(t) S(t) V (t)

The only assumption in deriving the above equations was that pressure was spatially uniform.

The viscous dissipation contribution is expected to be small, although when the fluid mixer
device is operating, the kinetic energy introduced into the fluid will eventually dissipate and lead
to an increase in thermal energy. For now, the following term will be carried:

S!visc = ∫ φ dV . (155)
V (t)

One further simplification is to treat the enthalpy along the various surface elements, Sk,
( !
)
as spatially constant hk ( r,t ) ≈ hk (t) for the integration over the control surface, S, so that

80
 ! ! ! ! ! !
∫ ( )
ph V − VS i n̂S dS = ∑ ∫ ( )
ph V − VS i n̂S dS ≈ ∑ hk ∫ ( )
h V − VS i n̂S dS = ∑ m! k hk . (156)
S(t) k Sk k Sk k

For consistency, the heat flux surface integral will be split into the same surface elements,
Sk, and be written as

! !
− ∫ q i n̂S dS = − ∑ ∫ q i n̂S dS =∑ qk , (157)
S(t) k Sk k

where Qk is the heat rate into the fluid through surface element, Sk.

!
Using the above relations consistent with the approximation hk ( r,t ) ≈ hk (t) , the thermal
energy equations in integral form become

d dV !
{me } + ∑ m! k hk = ∑ qk − P + S (158)
dt k k
dt
and

d dP !

dt
{ }
mh + ∑ m! k hk = ∑ qk +V
dt
+ Sk , (159)
k k

or

de dV !
m
dt k
( )
+ ∑ m! k hk − e = ∑ qk − P
dt
+ Sk (160)
k
and

dh dP !
m
dt k
( )
+ ∑ m! k hk − h = ∑ qk +V
dt
+ Sk . (161)
k

81
 APPENDIX B—HEAT AND MASS TRANSFER INSIDE THE TANK

B.1 Ullage-Bulk Liquid Interface Temperature

Temperature of the interface is calculated by Alabovskii’s equation, suggested in reference 1:

0.5 0.5
Θ u kl ⎛ µu ρl ⎞ ⎛ µl ⎞
= (1+ m ) , (162)
Θ l ku ⎜⎝ µl ρu ⎟⎠ ⎜ρ a ⎟
⎝ l l⎠
where

Θ u = Tu − Tint = ullage temperature difference, K

Θ l = Tl − Tint = liquid temperature difference, K

µu , µl = viscosity of ullage and liquid respectively, kg/m·s

C pl ρl
al = = liquid thermal diffusivity, m2/s
kl

qevap
m= = evaporation to total heat flows through interface ratio.
qintl

Because interface temperature, in most cases, is higher than liquid temperature and lower than
ullage temperature, and the right side of equation (162) is positive, the temperature difference
for  the liquid needs to be reassigned to the positive value as:

Θ l = Tint − Tl . (163)

After rearranging, equation (162) can be written as:

0.5
Tu − Tint kl ⎛ µu ρl ⎞
= ⎜ Prl ⎟ (1+ m ) . (164)
Tint − Tl ku ⎝ µl ρu ⎠
Suppose
0.5
k ⎛µ ρ ⎞
K = l ⎜ u l Prl ⎟ (1+ m ) , (165)
ku ⎝ µl ρu ⎠

82
and equation (164) is resolved for the interface temperature, Tint  ,

KTl + Tu Tl + Tu / K
Tint = = . (166)
K +1 1+ 1/ K

Table 3 presents coefficient K for some cryogenic fluids; hydrogen, oxygen, methane, and
nitrogen are calculated with assuming (1+m) is close to unity, which implies evaporation from the
interface is relatively low. As can be seen from equations (165) and (166), if evaporation is not neg-
ligible, the interface temperature will be even closer to liquid temperature.

Table 3. Coefficient K from equation (165).

Saturation Saturation
Fluid Temperature (K) Coefficient K Fluid Temperature (K) Coefficient K
Hydrogen Oxygen
18 20.16 85 116.42
19 17.24 88 95.21
20 14.88 90 83.79
21 12.94 95 62.01
Methane Nitrogen
105 99.60 75 87.80
110 76.31 78 70.45
115 59.51 80 61.24
120 47.10 82 53.51

As can be seen from table 3 and equation (166), the interface temperature is always closer
to the liquid temperature than to the ullage temperature.

B.2 Requirements for Using Flat Plate Correlations to the Cylindrical

Segment of the Tank

Conditions when the flat plate model can be used follow. For thermal conductivity (the heat
flow rate through the cylindrical wall) (see fig. 58), the tank is:

qcyl =
(
2π hcyl k T2 − T1 ),
(167)
(
ln r2 / r1 )
where hcyl  is the height of the wall and T2 > T1.

83
 L
r2
r1
T1 T2
r

T1 T2
hcyl
hcyl

Figure 58. Cylindrical to flat tank wall transformation.

F58_1627

For a flat wall, the heat flow rate can be written as:

q flat =
(
Ak T2 − T1 ) , (168)
δ

where A is a wall surface area, d is a thickness of the wall, and d = r2 – r1. It can be found
using the length of the cylindrical wall calculated with the average of the radius of the cylinder,
r = (r1 + r2) / 2, and height, hcyl , A = 2prhcyl , or

π ( r1 + r2 ) hcyl k (T2 + T1 )
q flat = . (169)
δ

Suppose temperatures, wall lengths and thicknesses, and wall heights are the same for both
cylindrical and flat walls. In this case, the relative difference between heat flow rates as a function
of the ratio of wall thickness to radius can be written as:

q flat − qcyl
Δ= , or Δ = 1− (0.5 + 1/ x)ln(1+ x), (170)
q flat

where x = d / r1.

84
 Figure 59 presents this function in the ratio range, which takes place in propellant tank
applications.

0.6

0.5
Relative Difference (%)

0.4

0.3

0.2

0.1

0
0 0.002 0.004 0.006 0.008 0.01
Wall Thickness-to-Radius Ratio, δ /r1

Figure 59. Relative difference between heat flow rates through circular and flat tank walls.
F59_1627

B.3 Free Convection Heat Transfer From the Tank Wall Inside the Tank

Following Elenbaas44, the Nusselt number can be accepted for the internal convection
in  a circular cylinder as:

C3
⎧ ⎡ ⎛ h 20 ⎞ C2 ⎤ ⎫
αr 1 r ⎪ ⎥ ⎪⎬
Nur = = Rar ⎨1− exp ⎢ − ⎜ . (171)
( )
kl Tw 16 h ⎪⎩ ⎣
⎟
⎢ ⎝ r Rar ⎠ ⎥ ⎪
⎦⎭

Equation (171) was developed, assuming a radius of cylinder r → ∞. Equation (171) should
match the Schmidt-Pohlhausen-Beckmann45 solution with the experimental adjustment for the flat
vertical plate:

αh
Nuh = = 0.6Rah 1/4 . (172)
( )
ku Tw
Using these equations,

gβCpu ρu ΔTh3 gβCpu ρu ΔTr 3

C2 i C3 = 3 / 4, Rah = , and Rar = . (173)
vku vku

and from equations (172), it can be found that:

Rar = Rah r 3 / h3 , Nur = Nuh r / h . (174)

85
Using equations (170)–(172), the cylindrical Nusselt to flat plate Nusselt ratio can be created:

C3
Nur 1 r ⎧ ⎡ ⎛ h 20 ⎞ C2 ⎤ ⎫
⎪ ⎥ ⎪⎬
= Rar ⎨1− exp ⎢ − ⎜ ⎟ Ra1/4
h , (175)
Nuh c1 h ⎪⎩ ⎢ ⎝ r Rar ⎠ ⎥ ⎪
⎣ ⎦⎭

where c1 = 16 × 0.6 = 9.6.

Simplifying equation (175),

C3
Nur 1 ⎛ r ⎞ 7/4 3/4 ⎧⎪ ⎡ ⎛ h 20 ⎞ C2 ⎤ ⎫
= ⎜ ⎟ Rar ⎨1− exp ⎢ − ⎜ ⎟ ⎥ ⎪⎬ . (176)
⎝
Nuh c1 h ⎠ ⎢ ⎝ r Rar ⎠ ⎥ ⎪
⎪⎩ ⎣ ⎦⎭

For cylindrical vertical channels, Elenbaas44 suggested C2 = 3/4, C3 = 1. With this
assumption, equation (176) becomes:

Nur 1 ⎛ r ⎞ 7/4 3/4 ⎧⎪ ⎡ ⎛ h 20 ⎞ 3/4 ⎤ ⎫

= ⎜ ⎟ Rar ⎨1− exp ⎢ − ⎜ ⎟ ⎥ ⎪⎬ . (177)
Nuh c1 ⎝ h ⎠ ⎪⎩ ⎢ ⎝ r Rar ⎠ ⎥ ⎪
⎣ ⎦⎭

For a small argument, the exponent in equation (177) can be expended into the Taylor
series. Keeping linear and quadratic terms only, gives:

⎡ ⎛ h 20 ⎞ 3/4 ⎤ ⎛ h 20 ⎞
3/4
1 ⎛ h 20 ⎞
3/2
exp ⎢ − ⎜ ⎥ ≈ 1− ⎜ ⎟ + 2 ⎜ r Ra ⎟ . (178)
⎢ ⎝ r Rar ⎟⎠ ⎥ ⎝ r Rar ⎠ ⎝ r ⎠
⎣ ⎦

3/4
⎛ h 20 ⎞
This series converges when ⎜ ⎟ < 1. Using equation (174), it can be found that
⎝ r Rar ⎠
1/4
r ⎛ 20 ⎞
> . (179)
h ⎜⎝ Rah ⎟⎠

1/4
r ⎛ 20 ⎞
The lowest ratio value, = ⎜ , is depicted in figure 60.
h ⎝ Rah ⎟⎠

86
 0.25

0.2

0.15
r/h

0.1

0.05

0
1×104 1×105 1×106 1×107 1×108 1×109 1×1010 1×1011
Rah

F60_1627
Figure 60. Minimal radius-to-height ratio value allows using the flat wall
convective heat flow model for cylindrical walls.

By substituting equation (178) into equation (177),

3/4 ⎡ 3/4 ⎤ 3/4 ⎡ 1 ⎛ h 20 ⎞ 3/4 ⎤

Nur 1 ⎛ r ⎞ 7/4 3/4 ⎛ h 20 ⎞ ⎛ ⎞
= ⎜ ⎟ Rar ⎜ ⎟ ⎢1− 1 ⎜ h 20 ⎟ ⎥= 20 r ⎢1− ⎥ . (180)
Nuh c1 ⎝ h ⎠ ⎝ r Rar ⎠ ⎢ 2 ⎝ r Rar ⎠ ⎥ c1 h ⎢ 2 ⎜⎝ r Rar ⎟⎠ ⎥
⎣ ⎦ ⎣ ⎦

Using equation (174),

r Nur Nur→h r
Nur = Nur→h , or = . (181)
h Nuh Nuh h

Now, equation (180) can be rewritten as:

3/4 3/4
Nur→h r 203/4 r ⎡ 1 ⎛ h 20 ⎞ ⎤ Nur→h 203/4 ⎡ 1 ⎛ h 20 ⎞ ⎤
= ⎢1− ⎥ , or = ⎢1− ⎥ . (182)
Nuh h c1 h ⎢ 2 ⎜⎝ r Rar ⎟⎠ ⎥ Nuh c1 ⎢ 2 ⎜⎝ r Rar ⎟⎠ ⎥
⎣ ⎦ ⎣ ⎦

203/4
Because this ratio should be equal to unity for large radius-to-height ratios, = 1.
c1
Calculating Nu relative difference and using equation (182),

87
 Nur→h 203/4 ⎛ h⎞
3
−3/4 ⎛ h⎞
3
−3/4
δ Nu = 1− = ⎜⎝ ⎟⎠ Ra h ≈ 4.7287 ⎜⎝ ⎟⎠ Rah . (183)
Nuh 2 r r

The relative difference, dNu, calculated by equation (182) for 104 ≤ Rar ≤ 1011 is presented in figure 61.

0.9
Relative Difference, 1–Nur→h /Nuh (%)

0.8 1×104
0.7 1×105

0.6 1×106
1×107
0.5
1×109
0.4 1×1011
0.3

0.2

0.1

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0
Cylinder Radius-to-Height Ratio (r/h)

Figure 61. Relative Nusselt numbers difference versus radius-to-height ratio for different
Rayleigh numbers.

B.4 Ullage Wall Condensation Model

B.4.1 Vertical Cylindrical Wall

Condensation on the internal side of the cylindrical wall segment is shown in figure 62.

For this system with the boundary layer approximation, the momentum equation can be
expressed as reported by Incropera and DeWitt46:

∂2 u g

∂y 2
=−
µl l
(
ρ − ρv . (184) )

Integrating twice with boundary conditions ∂u = 0; u(0) = 0, results in:

∂y δ = 0

u( y) =
(
g ρl − ρv ⎡
⎢ yδ −
y2 ⎤ )
⎥ . (185)
µl ⎢⎣ 2 ⎦⎥
88
 δ 0(x)
y
v
⋅
m(x)
u
dq=h fg dm⋅
dx dq
dm⋅
g
⋅ ⋅
m+dm
Wall, Tw δ (x)

Vapor, Tsat

Figure 62. Condensation on the internal side of cylindrical wall segment.

F62_1627

Integrating this velocity across the film, results in the following relation for the mass flow rate per
unit width of surface:

δ
( )
ρl ρl − ρv g ⎛ y2δ y3 ⎞
δ
( )
ρl ρl − ρv g ⎡ 3 δ 02
(
⎤
)⎥ . (186)
m! ′ = ρl ∫ udy =
µl ⎜ 2 − 6⎟
⎝ ⎠
=
3µl
⎢δ −
⎢ 2
3δ − δ 0
⎥⎦
δ0 δ0 ⎣

Differentiating equation (185) with respect to δ yields:

=
( ⎢
)
−
2
d m! ′ ρl ρl − ρv g ⎡ 2 δ 0 ⎤
δ ⎥ . (187)
dδ µl ⎢⎣ 2 ⎥⎦

Assuming that heat flux across the film is mainly due to conduction, and using conservation
of  energy from the condensation, one can get:

kl ΔT
dx = h′fg d m! ′, (188)
δ

where referring to Rohsenow47, h′fg = h fg (1+ 0.68Ja) and Ja =

(
C pl Tsat − Tw ).
h fg

Combining equations (187) and (188), the following differential equation is obtained for δ:

dx =
(
ρl ρl − ρv gh′fg ) ⎡ δ 2δ ⎤
dδ ⎢δ 3 − 0 ⎥ . (189)
µl kl ΔT ⎢⎣ 2 ⎥⎦

89
Integrating equation (189) and solving for δ, the following can be obtained:

1/2 ⎤ 1/2

⎡ 2 4
⎢ δ 0 + δ 0 + 4kx
δ (x) = ⎢
( ⎥ )
2 ⎥ , (190)
⎢⎣ ⎥⎦

4 µl kl ΔT
where k = .
( )
ρl ρl − ρv gh′fg

Because heat transfer across the film is by conduction alone, the heat transfer coefficient
is  given by
1/2
⎛ ⎞
2
hl (x) = kl / δ (x) = kl ⎜ ⎟ . (191)
⎜⎝ δ 2 + δ 4 + 4kx ⎟⎠
0 0

The mean heat transfer coefficient is:

1/2
L ⎛ ⎞
1 k L 2 k L
1
hl = ∫ h l (x) dx = l ∫ ⎜ ⎟ dx = 1/4l ∫ dx , (192)
L
0
L ⎜ δ 2 + δ 4 + 4kx ⎟
0⎝ 0 0 ⎠ k L 0 ( a + a+x )1/2

where a = δ 04 4k.

Integrating equation (192), one gets:

1/4
( 3
4 ⎡ g ρl ρl − ρv kl h′fg ⎤)
hl = ⎢
3 2⎢ µl ΔTL
⎥
⎥
⎡
⎣⎢
( b +1 − 2 b )( b + b +1 )1/2 + 2b3/4 ⎤ , (193)
⎦⎥
⎣ ⎦

where b = δ 04 4kL.

Supposing d0 = 0 results in the Nusselt equation for the mean laminar heat transfer
coefficient on the vertical wall,

1/4

hl = 0.9428 ⎢
(
⎡ g ρ ρ − ρ k 3h ′ ⎤
l l v l fg
⎥
) . (194)
⎢ µl ΔTL ⎥
⎣ ⎦

90
The average Nusselt number can be expressed in the form:

1/4
hlL (
⎡ g ρ ρ − ρ h′ L3 ⎤)
Nul =
kl
= 0.9428 ⎢
l l
(
v fg
⎥
⎢ kl µl Tsat − Tw ⎥ )
⎡
⎢⎣ ( b +1 − 2 b )( b + b +1 )1/2 + 2b3/4 ⎤ . (195)
⎥⎦
⎣ ⎦

The thickness of the film will be:

1/4
⎛ x⎞
δ = (kL)1/4 ⎜ b + b + ⎟
1/2
= ⎜ l l sat
(
⎛ 4µ k T − T L ⎞
w ) ⎛ x⎞
1/2
⎟ ⎜⎝ b + b + L ⎟⎠ . (196)
⎝ L⎠ (
⎜⎝ ρl ρl − ρv g h′fg ⎟⎠ )

The dimensionless parameter b that depends on the initial liquid film thickness is:

b= =
(
δ 04 δ 04 ρl ρl − ρv gh′fg
. (197)
)
4kL 16 µl kl Tsat − Tw L ( )

B.4.2 Upper Dome

The condensation appears in the underside of inclined surfaces for this case. Following
Gerstmann and Griffith7, the following assumptions were made:

• Flow is laminar.
• Angle of inclination is small.
• Pressure in the film is hydrostatic.
• Vapor exerts negligible shear on the interface.
• Wall temperature and the vapor temperature are uniform.
• Group kl  DT / ml   hfg is much less than unity.
• Group Ja = cpl DT/ hfg is much less than unity.
• Capillary length scale is accepted in the form:

1/2
⎡ σ ⎤
ψ =⎢ ⎥ . (198)
( )
⎢⎣ g ρl − ρv cos(α ) ⎥⎦

91
Using previous assumptions and figure 63, the governing equations may be written as:

∂p ∂2 vx
= µl , (199)
∂x ∂x 2

δ
α
y

z g

Figure 63. Condensation on the upper dome.

F63_1627

d 2δ
( )
p = psat − ρl − ρv g cos(α )(δ − z) − σ
dx 2
, (200)

∂vx ∂vz
+ = 0 , (201)
∂x ∂z
and

=
( .
)
Q kl Tsat − Tw (202)
A δ

By differentiating equation (200) and substituting in equation (199), the following equation can be
obtained:

∂2 vx dδ d 3δ
µl
∂x 2
= − ρ l(− ρ v g cos )
(α )
dx
− σ
dx 3
. (203)

Integrating equation (203) with boundary conditions of zero velocity at the wall and zero
shear stress at the interface gives:

1 ⎡ dδ d 3δ ⎤ ⎛ z2 ⎞

µl ⎢⎣
( )
vx = ⎢ ρl − ρv g cos (α )
dx
+ σ 3 ⎥ ⎜ zδ − ⎟ . (204)
dx ⎥⎦ ⎝ 2⎠

Integrating equation (200) over the film thickness, using equations (201) and (203),
equation (205) can be obtained:

92
 d ⎧⎪ 1 ⎡ ⎤ ⎫⎪
2
dδ d 3δ ⎤ ⎡ 3 δ 0 kl ΔT
δ ⎨
dx ⎪ 3µl ⎢⎣
( )
⎢ ρl − ρv g cos (α )
dx
+ σ 3 ⎥ ⎢δ −
dx ⎥⎦ ⎢⎣ 2
(
3δ − δ 0 ⎥ ⎬ = ) ρl h fg
. (205)
⎩ ⎥⎦ ⎪⎭

⌢ δ
Equation (204) can be made dimensionless by defining the following new variables: δ =
⌢ x ψ
and x = .
ψ
Using primes to signify differentiation with respect to x, equation (204) may be restated as
⌢ ⌢
⌢ ⎡ ⌢ 3 δ 02 ⌢ ⌢ ⎤ ⌢ ⌢ ⌢ ⎛ ⌢ 2 δ 02 ⎞ ⎡δ⌢′′′δ⌢′ + δ⌢′ 2 ⎤ = Θ , (206)
δ ⎢δ −
⎢⎣ 2
(
3δ − δ 0 )(
⎥⎦
)
⎥ δ ′′ + δ ′′′′ + 3δ ⎜ δ − ⎟
⎝ 2 ⎠ ⎢⎣ ( ) ⎥⎦

where
3kl ΔT µl
Θ= .
ψ ρl h fg σ

⌢ ⌢
Assuming that δ 02 ≪ δ 2 in equation (206) results in the same equation as obtained by Gerstmann7:

⌢ ⌢ ⌢ ⌢ ⌢ ⌢ ⌢ 2
( ) ( )
δ 4 δ ′′ + δ ′′′′ + 3δ 3 ⎡⎢δ ′′′δ ′ + δ ′ ⎤⎥ = Θ . (207)
⎣ ⎦

The solution of this equation is given in reference 6 for Ra > 106:

0.9
Nu = ,

(
Ra1/6 1+ 1.1Ra −1/6 ) (208)

where Ra is a modified Rayleigh number:

1/2
ρlσ h fg ⎡ σ ⎤
Ra = ⎢ ⎥ . (209)
( ) ( )
kl µl Tsat − Tw ⎢⎣ g ρl − ρv cos(α ) ⎥⎦

It was recommended by Gerstmann and Griffith7 that equations (208) and (209) could be used
for small inclinations, a < 20o.

For angles 20o < a < 90o, the modified Nusselt model was recommended. At this model,
the  free-fall acceleration g in equations (195)–(197) should be replaced by gsin(a). In this case,
equations (195)–(197) result in:

93
 1/4
hl L ( )
⎡ g sin(b) ρ ρ − ρ h′ L3 ⎤
Nul =
kl
= 0.9428 ⎢
⎢ (
l l v fg
kl µl Tsat − Tw
⎥
⎥ )
⎡
⎣⎢
( b +1 − 2 b )( b + b +1 )1/2 + 2b3/4 ⎤ , (210)
⎦⎥
⎣ ⎦

1/4

δ =⎜ l l sat(
⎛ 4µ k T − T L ⎞
w
⎟
) ⎛
b + b +
x⎞
1/2
, (211)
( )
⎜⎝ ρl ρl − ρv g sin(α )h′f g ⎟⎠ ⎜⎝ L ⎟⎠

and

b= 0
(
δ 4 ρl ρl − ρv g sin(α )h′f g
. (212)
)
16 µl kl Tsat − Tw L ( )
The tangent line equation for the ellipse is:

xx0 yy0
+ = 1. (213)
hb2 r2

From equation (213), points (0, y*) and (x*, 0) can be found:

hd2 r2
(0, y*) : y* = ; (x*,0) : x* = . (214)
y0 x0

From figure 64, coordinates x0 and y0 can be found:

h ⎛ h ⎞
y0 = hb − h1; x0 = r 1 ⎜ 2 − 1 ⎟ . (215)
hd ⎝ hd ⎠

y
(0, y*)

h1 α (x0, y0)
b, hb
α
(x*, 0)
x
a, r g

Figure 64. Angle a calculation schematic.

F64_1627

94
From figure 64, equations (214) and (215), it can be found that:

tan(α ) =
y * hd
=
(
h1 2 − h1 ) , (216)
x* r (1− h1)2
h
where h1 = 1 ,
hd

and

tan(α ) 1
sin(α ) = , cos(α ) = . (217)
2 2
1+ tan (α ) 1+ tan (α )

B.4.3 Lower Dome

As shown by Oosthuizen8, the condensation on the upper side of the inclined plate can be
described by the Nusselt equations, where the free-falling acceleration is replaced by the product
gsin(b). In this case, b is an angle between the tangent line and free-fall acceleration. For the con-
densation on the lower dome, equations (210)–(212) convert to:

1/4
hl L ⎡ g sin( β ) ρ ρ − ρ h′ L3 ⎤ ( )
Nul =
kl
= 0.9428 ⎢
⎢
l l v fg
kl µl Tsat − Tw (
⎥
⎥ )
⎡
⎢⎣ ( b +1 − 2 b )( b + b +1 )1/2 + 2b3/4 ⎤ (218)
⎥⎦
⎣ ⎦

and
1/4
⎛ x⎞
δ = ( kL )1/4 ⎜ b + b + ⎟
1/2 ⎛ 4µ k T − T L ⎞
=⎜ l l sat w ( ) ⎛ x⎞
1/2
⎟ ⎜⎝ b + b + L ⎟⎠ , (219)
⎝ L⎠ ( )
⎜⎝ ρl ρl − ρv g sin( β )h′f g ⎟⎠

where

b= =
( )
δ 04 δ 04 ρl ρl − ρv g sin( β )h′fg
,

4kL 16 µl kl Tsat − Tw L( )
4 µl kl ΔT
k= ,
( )
ρl ρl − ρv g sin( β )h′fg

h′fg = h fg (1+ 0.68Ja ),

Ja =
(
cpl Tsat − Tw ).
h fg

95
Figure 65 shows the condensation on the lower dome.

(x*, 0)

x
h1

(x0, y0)
β
α
β
(0, y*)
g

F65_1627
Figure 65. Condensation on the lower dome.

B.4.4 Angle b Calculation.

Because a + b = 90o, tan( β ) = cot(α ) =

1
, or tan( β ) =
r 1− hl
,
( )2
tan(α ) hb hl 2 − hl ( )
tan( β )
and sin( β ) = .
2
1+ tan ( β )

96
 APPENDIX C—TANK THERMODYNAMIC VENTING PRESSURE CONTROL SYSTEMS

C.1 Liquid Film Heat Transfer

To make calculations, the wall is divided into layers and the number of these layers are equal
to the number of sprays receiving the wall. For each layer, the mass flow rate, velocity of the liquid,
and thickness of the liquid film can be found. Assume that droplets in the ullage and liquid in the
film on the wall have no friction losses. In this case, the distance between two sequential droplets
reaching the wall can be found using the initial distance between the orifices and orthogonal com-
ponent of the droplet velocity (fig. 66).

h1
Δ h1
m⋅ 1
h2 Vdr1
vg1
Δ h2 H1
m⋅ 2
Vdr2
vg2
H2
m⋅ 3
Vdr3
vg3

vgf

Figure 66. Heat and mass transfer in theF66_1627

liquid film.

2 2
g⎛ l ⎞ g⎛ l ⎞
H1 = h1 + Δ h2 − Δ h1; Δ h1 = ⎜ ⎟ ; Δ h2 = ⎜ ⎟ . (220)
2 ⎜⎝ Vdr1 ⎟⎠ 2 ⎜⎝ Vdr 2 ⎟⎠

From equation (220), one can get:

gl 2 ⎛ 1 1 ⎞ gl 2 ⎛ 1 1 ⎞
H1 = h1 + ⎜ 2 − 2 ⎟ , H2 = h2 + ⎜ 2 − 2 ⎟ . (221)
2 ⎝ Vdr 2 Vdr1 ⎠ 2 ⎝ Vdr 3 Vdr2 ⎠

97
 In general, for each orifice Hi can be written as:

gl 2 ⎛ 1 1 ⎞
H i = hi + ⎜ 2 − 2 ⎟⎟
. (222)
2 ⎜⎝ Vdr(i+1) Vdri ⎠

In the beginning and in the end of each layer, the mass flow rates and velocities of the liquid
can be found. In the beginning of the layer,

m! b(i + 1) = m! fi + m! i + 1 − Δ m! imp(i + 1) (223)

and

m! fi vfi + ⎡ m! i + 1 − Δ m! imp(i + 1) ⎤ vi + 1
vb(i + 1) = ⎣ ⎦ . (224)
m! fi + m! i + 1 − Δ m! imp(i + 1) ⎤
⎡
⎣ ⎦

In the end of the layer,

m! f (i + 1) = m!b (i + 1) − Δ m! evap (i + 1)
(225)
and
vf (i + 1) = vb(i + 1) + 2gHi + 1 . (226)

98
 APPENDIX D—TANK AXIAL JET THERMODYNAMIC VENTING
PRESSURE CONTROL SYSTEMS

D.1 Heat and Mass Transfer on the Ullage-Liquid Interface During Axial Jet Mixing

To calculate heat and mass transfer on the interface, the following assumptions were made:

β g ΔTL
• Buoyancy effects are insignificant, Ri ≡ b 2 b < 1.
vb

µbc pb
• Bulk liquid Prandtl number is lower than 6, 1 < Prb ≡ < 6.
kb

cpb ΔT
• Bulk liquid Jacob number is small, Ja ≡ < 0.2.
h fg

ρ v L
• Ullage-liquid interface is turbulent, 350 < Re ≡ b b b < 11,000.
µb

ρj vj d
• Jet is fully turbulent, Re j ≡ > 2,300.
µj

• Jet nozzle diameter is small compared to the nozzle submergence, d << s.

• Lb is a turbulent linear scale.

From figure 67, the jet diameter at any height, x, can be calculated as D(x) = d + 2tan(a/2)x.
When a jet reaches the interface, its diameter is Dj = d + 2tan(a /2)s. If a diameter of the jet is larger
than the interface diameter, the result is high nozzle submergence; if smaller, low nozzle submer-
gence. The boundary between these two cases can be found when diameters are equal: Dj = Dint,
or Dint = d + 2tan(a /2)s.

99
 Dint

Dj

D(x)
s

α /2
x

Figure 67. Geometrical parameters of axial jet.

F67_1627

The boundary submergence, s/Dint, follows:

s 1 ⎛ d ⎞
= 1− . (227)
Dint 2tan(α / 2) ⎜⎝ Dint ⎟⎠

s 1 ⎛ d ⎞
From equation (227), it follows that < 1− is a low nozzle submergence
Dint 2tan(α / 2) ⎝ Dint ⎟⎠
⎜
s 1 ⎛ d ⎞
and ≥ 1− is a high nozzle submergence.
Dint 2tan(α / 2) ⎝ Dint ⎟⎠
⎜

As shown by Tollmien48, angle α of the jet expansion in the same liquid is about 11o ÷ 12o.
Dominick49 suggested correlation for the expansion angle, which for SI units is as follows:
⎛α⎞
tan ⎜ ⎟ = 0.8252(v)0.135 , (228)
⎝ 2⎠

where v is a kinematic viscosity of liquid, m2/s.

100
 D.2 Helicoidally Coiled Tube-in-Shell Heat Exchanger Transfer

Hydraulic diameter for the flow inside a coiled pipe can be found by the standard equation
Dhvp = 4A/P.

From figure 68, the area of ellipse is A = pDd/4. The first Ramanujan’s formula for the
perimeter of the ellipse follows:

Dec

Dc

Dic

Dis

Des

Figure 68. Geometrical parmeters of helicoidally

F68_1627
coiled tube-in-shell heat exchanger.

π
P= ⎡3(D + d ) − (3a + b)(a + 3b) ⎤⎦ (229)
2⎣
can be modified to

P=
π
2
( )
(D + d ) 3 − 4 − hin , (230)

where hin = [( D − d ) / ( D + d )] .(from fig. 26(a)).

2

101
 From equation (230), the hydraulic diameter for the elliptical tube can be found:

4A 2Dd
Dhv ρ = = . (231)
(
P ( D + d ) (3 − 4 − h in )
As can be seen in figure 68, the open area on the shell can be calculated as:

As =
π⎡ 2
4 ⎣
Des − Dis2 − Dec
2
( 2 ⎤
− Dic
⎦
. (232) )
From figure 26(a), Dec = Dc + d + 2d and Dic = Dc – d – 2d. In this case, the open area is:

As =
π⎡ 2
4 ⎣ (
Des − Dis2 − 4Dc d + 2δ vp ⎤ . (233)
⎦ )
The mass flux in the shell side can be found by:

4 m! shell
Gshell = − . (234)
π ⎡ Des
⎣
2
(
– Dis2 − 4Dc d + 2δ vp ⎤
⎦ )
Hydraulic diameter for the shell side can be found by Dhs = 4V/S, as suggested by Patil50,
which is equivalent to the standard formula Dh = 4A/P for the cylindrical geometry. The free vol-
ume for the shell side can be calculated as follows (assume the heat exchanger takes n coils with
total length, L):

2
⎛ π Dc ⎞
L=n π 2
Dc2 2
+ p = np 1+ ⎜ ⎟ = np 1+ γ −2 , (235)
⎝ p ⎠

where

p  is pitch of a helicoidally coiled pipe, m

γ  = p/pDc is a dimensionless pitch.

With this assumption, the open (free) volume of the shell side is:

Vs = np
π 2
4
( π
) (
Des − Dis2 − np ⎡ D + 2δ vp (d + 2δ ) 1+ γ −2 ⎤
4 ⎣⎢ ⎦⎥
)
= np
π
4 {( 2
Des
⎣ ) ( ⎦)
− Dis2 − ⎡ Dd + 2δ vp (D + d ) ⎤ 1+ γ −2 . (236) }

102
Total area of the shell side can be calculated as:

( ) ( )
As = npπ Des + Dis + npPs 1+ γ −2 = np ⎡π Des + Dis + Ps 1+ γ −2 ⎤ . (237)
⎢⎣ ⎥⎦

Replacing the external perimeter of the coiled pipe by Ramanujan’s formula gives:

( π
) (
As = np ⎡⎢π Des + Dis + np D + d + 4δ vp 3 − 4 − hs
⎣ 2
)( ) 1+ γ −2 ⎤⎥ . (238)
⎦

In this case, the shell side hydraulic diameter is:

Dis2 − Des
2 ⎡
− Dd + 2δ vp (D + d ) ⎤ 1+ γ −2
Dhs = 2 ⎣ ⎦ , (239)
( ) (
2 Dis + Des + D + d + 4δ vp 3 − 4 − his)( ) 1+ γ −2

( )
2
where his = ⎡( D − d ) D + d + 4δ vp ⎤ .
⎣ ⎦

D.3 Heat Transfer at Low Acceleration Conditions

Assuming that the plot of the approximation function would go through three points
of the graph calculated by Surface Evolver, from figure 28 the following can be found:

m m m
⎛ Bo ⎞ ⎛ Bo ⎞ ⎛ Bo ⎞
1− ⎜ 1 Bo ⎟ 1− ⎜ 2 Bo ⎟ 1− ⎜ 3 Bo ⎟
⎝ 3⎠ ⎝ 3⎠ ⎝ 3⎠
A1 = a + b m
, A2 = a + b m
, and A3 = a + b m
. (240)
⎛ Bo ⎞ ⎛ Bo ⎞ ⎛ Bo ⎞
1+ ka ⎜ 1 Bo ⎟ 1+ ka ⎜ 2 Bo ⎟ 1+ ka ⎜ 3 Bo ⎟
⎝ 3⎠ ⎝ 3⎠ ⎝ 3⎠

Solving the system of equations (240), it can be found that

a = A3 , (241)

ka =
( A1 − A3)(1− Bo23 ) − ( A2 − A3)(1− Bo13 ) , (242)
( A2 − A3)(1− Bo13 ) Bo23 − ( A1 − A3)(1− Bo23 ) Bo13
and
⎛ 1+ ka Bo23 ⎞
(
b = A1 − A3 ⎜ )
⎝ 1− Bo23 ⎠
⎟ , (243)

103
where
m m
Bo
Bo13 = ⎛⎜ 1 Bo ⎞⎟ and Bo23 = ⎛⎜ Bo2 ⎞
,
⎝ 3⎠ ⎝ Bo3 ⎟⎠

A1, A2, A3 – interface areas or dry wall areas corresponding to three chosen points
on the graph
Bo1, Bo2, Bo3 – Bond numbers corresponding to the same three chosen points on the graph
m – a coefficient that can be found with the least squares method or graphically
for the best match with the Surface Evolver data.

Point 1 on figure 28 corresponds to very low Bond numbers (actually, Bo = 0), when all
gas volume is inside the liquid phase and the dry wall area disappears. Point 3 corresponds to
high Bond numbers, where gas-liquid interface is flat and both relative areas are equal to the unit.
Point 3 can be taken as inflection points of the graph (fig. 28(a)) or as a minimum point on the
interface area graph and an inflection point on the dry wall area (fig. 28(b)). On both interface
and dry area lines, point 2 can be different.

104
 APPENDIX E—COMPUTER MODEL DESCRIPTION

E.1 Input Files Description

The following eight items describe the input files:

(1) Address file (fig. 69)

– ‘TankSIM_Addr_Input.txt’
– Includes Input-Output folder address and file number assignments for TankSIM input
and output open file procedures.
– Should be in the same folder as the TankSIM (‘TankSIM.f 90’) source file.
All files in figure 69 should be in the In-Out folder.

Figure 69. TankSIM address input file.

(2) TankSIM input file (fig. 70)

– ‘TankSIM_Input.txt’
– Includes fluid, initial conditions—pressures, temperatures, flow rates, external heat
flows, calculation time, and time steps for regular and mixing regimes.

105
 Figure 70. TankSIM input file.

(3) Mission profile input file (fig. 71)

– ‘TankSIM_Mission_Profile.txt’
– Includes mission profile information—number of mission phases, regimes informa-
tion, mission time, time steps, heat loads, pressure regulation ranges, flow rates, pump parameters,
and other information.

Figure 71. TankSIM mission profile input file.

106
Figure 71. TankSIM mission profile input file (Continued).

107
 (4) Tank input file (fig. 72)
– ‘TankSIM_Tank.txt’
– Includes tank information—material, shape, tank measurements, and tank material
specific heat and heat conduction table as a function of temperature.

Figure 72. TankSIM tank input file.

108
 (5) Spray bar, heat exchanger, and pump input file (fig. 73)
– ‘TankSIM_SprayBar.txt’
– Includes spray manifold, injection pipes, and heat exchanger measurements
and properties.

Figure 73. TankSIM spray bar input file.

109
 (6) Axial jet input file (fig. 74)
– ‘TankSIM_AxialJet.txt’
– Includes axial jet nozzle and helicoidally coiled heat exchanger measurements
and pressure drop coefficients.

Figure 74. TankSIM axial jet input file.

(7) Optional external lines input (fig. 75)

– ‘TankSIM_Tank_Add_Feed.txt’
– ‘TankSIM_Tank_Add_Pres.txt’
– ‘TankSIM_Tank_Add_Vent.txt’ (example of this file presented in fig. 75)
– These files have the same structure and include internal volume, connected
to the ullage, environmental and approximate pipe temperatures, line segment measurements,
and constant thermophysical properties.

110
 Figure 75. TankSIM external venting line input file.

(8) Main output file (fig. 76)

– ‘TankSIM_20160830_110721_ 1_Out-H2_90%d3_9-50_2W_81d.csv’
– Includes all initial and calculated information. The length of the record (line) may
consist of 55 columns of data depending on TVS type.

Figure 76. TankSIM main output file.

111
 E.2 Variables Description for Input and Output Files

Table 4 lists the description of the variables for the input and output files.

Table 4. Variables description for input and output files.

Description of TankSIM Input and Output Files Variables

Address File (TankSIM_Addr_Input.txt)

Path Path to the InOut folder 40 pos.

PathServ Path to the folder for service files 40 pos.
FN1-FN17 Working file numbers –

Mission Profile File (TankSIM_Mission_Profile.txt)

Common
flgmis Mission phase flag: 'Slf' - self pressurization; 'Chl' - feed line chilldown; 3 pos.
'Hep' - non-condensable gas pressurization; 'Fir' - engine firing; 'Lqf' - wall liquefaction
dtmiss Mission duration sec
tmstep Time step for calculation sec
prtmiss Number of records, skipping during printing –

'Slf'
qext Uniformly distributed heat loads W
qheatu Constant flow rate to the ullage (usually from supporting elements, manholes, etc.) W
qheatl Constant flow rate to the liquid (usually from supporting elements, manholes, etc.) W
gr Vehicle acceleration-to-g ratio –
teta Angle between thrust vector and tank longitudinal axis degree
pumin Minimum required ullage pressure kPa
pumax Maximum required ullage pressure kPa
plmin Minimum required liquid saturation pressure kPa
plmax Maximum required ullage liquid saturation pressure kPa
flgpcl Ullage pressure control logic flag: -
1 - control by ullage and liquid saturation pressures independently
2 - control by ullage pressure only
3 - control by liquid saturation pressure only
regime Pressure control regime: 4 pos.
'barm' - spray bar mixing, no venting; 'barh' - spray bar TVS;
'jetm' - axial jet mixing, no venting; 'jeth' - axial jet TVS;
'ullg' - ullage venting; 'sbul' - ullage venting and spray bar mixing;
'preg' - ullage pressure regulation;
'prbm' ---“--- with spray bar mixing;
'prbh' ---“--- spray bar TVS;
'prjm' ---“--- axial jet mixing;
'prjh' ---“--- axial jet TVS;
'prul' ---“--- ullage venting;
'boff' - boil-off
flgback Back pressure device: 'orifice' - critical orifice; 'venturi' - critical Venturi; ‘prconst’ - constant backpressure 7 pos.

112
 Description of TankSIM Input and Output Files Variables (Continued)

kfri Friction pressure lost coefficient from the pump to the heat exchanger (vent line)
kdevi Device pressure lost coefficient from the pump to the vent line (valves, flowmeters)
mdvent Initial vent flow rate kg/sec
Lohm Joule-Thomson device resistance (for Lee equation) 1/ft2
pbacki Initial backpressure kPa
dback Diameter of backpressure device mm
cdor Back pressure orifice discharge coefficient –
deltmix Time step for calculation during TVS functioning sec
iprtmix Number of records, skipping during printing (TVS functioning)
vdgpm Pump volumetric flow rate gpm
p0 Constant pump Delta P (all other coefficients = 0) or free term in Delta P vs. Flow Rate curve kPa
p1…p5 Coefficients for Delta P (kPa) vs. Flow Rate (gpm) curve ( 5-th order maximum) –
e0 Constant pump Power (all other coefficients = 0) free term in Power vs. Flow Rate curve W
e1…e5 Coefficients for Power (W) vs. Flow Rate (gpm) curve ( 5-th order maximum) –

'Chl'
mdch Feed line chilling down mass flow rate kg/sec

'Hep'
puhp Nominal start box pressure kPa
dpuhp Start box pressure margin (ratio) –

'Fir'
mdlf Propellant mass flow rate during firing kg/sec
pufir Nominal run box pressure kPa
dpufir Run box pressure margin (ratio) –
flgprs Pressurization flag: 0 - autogenous; 1 - non-condensable –
pvp Pressure of pressurization vapor or gas entering tank kPa
tvp Temperature of pressurization vapor or gas entering tank K
hvp Enthalpy of pressurization vapor or gas entering tank kJ/kg
gr Vehicle acceleration-to-g ratio –
teta Angle between thrust vector and tank longitudinal axis degree

'Lqf'
qext Uniformly distributed heat loads W
qcrc Uniformly distributed heat eliminating by cryocooler W
qheatu Constant flow rate to the ullage (usually from supporting elements, manholes, etc.) W
qheatl Constant flow rate to the liquid (usually from supporting elements, manholes, etc.) W
gr Vehicle acceleration-to-g ratio –
teta Angle between thrust vector and tank longitudinal axis degree
mdlqf Propellant vapor mass flow rate liquefaction kg/sec
pulqf Nominal run box pressure kPa
dpulqf Run box pressure margin (ratio) –
pvp Pressure of vapor entering tank kPa
tvp Temperature of vapor entering tank K
hvp Enthalpy of vapor entering tank kJ/kg
Input File (TankSIM_Input.txt)

fluid Working fluid: ''Methane_', 'Hydrogen', 'Oxygen__', 'Nitrogen', 'Parahydr' 8 pos.

113
 Description of TankSIM Input and Output Files Variables (Continued)
pheini Initial pressurization gas pressure kPa
pvini Initial vapor pressure kPa
tuini Initial ullage temperature K
tlini Initial liquid temperature K
twini Initial ullage wall temperature K
tlwini Initial liquid wall temperature K
tkwlini Initial thickness of ullage wall liquid film m
full Liquid fill level in volume percent %
phest Bottle helium pressure in the beginning of the mission kPa
thesti Bottle helium temperature in the beginning of the mission K
dhest Helium regulator orifice diameter mm
tkhest Thickness of regulator orifice diameter mm
vhest Helium bottles total volume m3
brktime Program breaking time, used for debugging sec
iplotreg Number of time steps between output plotting, pump is off NOT USED NOW
iplotmix Number of time steps between output plotting, pump is on NOT USED NOW
itesti Initial printing flag (usually use 5) –
flgprt Print flag: –
0 - only inputs and summary outputs, no details;
1 - basic printing with details;
2 - Ullage venting printing;
3 - Spray bar TVS printing;
4 - Axial jet TVS printing;
5 - Combination of ullage venting & spray bar TVS;
6 - Combination of ullage venting & axial jet TVS;
7 - Liquefaction;
tent Environmental temperature K
penv Environmental pressure kPa
casnum Any text, used for file name and naming a case 20 pos

Tank File (TankSIM_Tank.txt)

tmater Tank wall material: ' SS 304 ', 'Aluminum', '  Titanium ', etc. 8 pos.
tshape Tank bulkhead shape: ' Spherical  ', ' Elliptical ', ' Flat ' 10 pos.
flgadd ' 1 ' indicates additional external lines attached to the tank wall: feed, vent, pressure. –
dtank Tank outside diameter m
hcyl Tank cylinder part (barrel) height m
hbulk Tank bulkhead height m
tkw Tank wall thickness m
vaddt Additional volume added to the Ullage from additional lines m3
ptmax Tank design maximum pressure kPa
rhow Tank material average density kg/m3
ksi Surface average heat absorption coeff. (used in Stephan-Boltzmann equation) –
nmet Number of points in tank material properties table –
tmet Temperature in tank material properties table K
cpmet Tank material specific heat table J/(kg•K)
kmet Tank material thermal conductivity table W/(m•K)

114
 Description of TankSIM Input and Output Files Variables (Continued)
Spray Bar File (TankSIM_Spray_Bar.txt)

dmno spray manifold external diameter m

thmn spray manifold wall thickness m
lmn spray manifold length m
kmn spray manifold material average thermal conductivity W/(m•K)
ltop distance from top of the tank to top part of spray manifold (flgtop = 1) or m
distance from top of the tank to the top of the injection line manifold (flgtop = 2)
nbar number of spray bars in tank –
kmni total pressure lost coefficient from the pump to the spray manifold –
kmnf total pressure lost coefficient from the spray manifold to injection lines –
dito injection line external diameter m
thit injection line wall thickness m
lit injection line length m
dzsi distance between orifices on the injection lines m
norsi number of orifice levels on the injection line –
norit number of orifices on the each orifice level –
ninj number of injection lines per spray bar –
ditom distance from the top to the first top orifice (flgtop = 1), or injection line manifold external diameter (flgtop = 2) m

115
 Description of TankSIM Input and Output Files Variables (Continued)
thitm Manifold height (flgtop = 1), or injection line manifold wall thickness (flgtop = 2) m
litm Spray manifold top part length (flgtop = 1), or injection line manifold length (flgtop = 2) m
dzsim Distance between orifices on the top part (flgtop = 1), or injection lines manifold (flgtop = 2) m
nortop Number of orifice levels on the top part of the spray bar manifold (flgtop = 1), or –
on the injection line manifold (flgtop = 2)
norimn Number of orifices on the each orifice level on the top part of spray bar manifold –
(flgtop = 1), or on the injection lines manifold (flgtop = 2)
lspr Distance from the injection line spray orifice to the tank wall (used for calculating m
amounts of liquid spraying to the wall and falling to the bulk liquid)
flgtop Type of the spray bar top part: 1 - MHTB with perforated spray manifold top part; –
2 - SLS type with perforated injection line manifold (large tanks)
ks Spray orifice loss coefficient –
cds Spray orifice discharge coefficient –
nlim Maximum number of iterations –
nsec Number of sections with the same orifice sizes –
node Node number (number of different orifices in injection line table –
(if all same orifices - node = 1)
dorf Orifice nodes diameter table m
dexo Heat exchanger external diameter m
thex Heat exchanger wall thickness m
lex Heat exchanger length m
xei Heat exchanger initial gas quality –
nex Number of points heat exchanger divided –
alhex Angle between heat exchanger and tank longitudinal axis degree
kvent Kvent = 1, vent line begins after pump (flow rate subtracts from pump flow rate) –

Axial Jet File (TankSIM_AxialJet.txt)

bet1 Coefficient in axial jet equation (0.34) –

bet2 Coefficient in axial jet equation (0.24) –
hj Axial jet nozzle exit height m
dj Axial jet internal nozzle diameter m
nzj Axial jet number of nozzles –
dexol Heat exchanger pipe minor axis m

TankSIM Output File

File Name of the current file –

Fluid Working fluid –
Thank Measurement
Same as in tank input file
Mission Profile
Same as in mission profile input file
Regular printing for 'Slf' and common part for all other regimes
Regime Working regime during phase: 3 pos.
'Slf' - self pressurization (no mixing and TVS working);
'Chl' - lines chill-down );

116
 Description of TankSIM Input and Output Files Variables (Continued)
'Chp' - chill-down with autogenuous pressurization;
'Hep' - non-condensable gas pressurization;
'Hec' - pressurization phase with zero non-condensable flow rate (no pressurization
during pressurizaton phase, ullage pressure higher then required max);
'Fir' - engine firing without pressurization;
'Fip' firing with autogenous pressurization;
'Bar' - spray bar TVS is switched on (cycling);
'MiB' spray bar only mixing (cycling);
'PhB' - spray bar TVS working permanently;
'PmB' spray bar only mixing permanently;
'Jet' - axial jet TVS is switched on (cycling);
'MiJ' - axial jet only mixing (cycling);
'PhJ' - axial jet TVS working permanently;
'PmJ' - axial jet only mixing permanently;
'Ull' - ullage venting (open valve).
'Lqf' - liquefaction
Time Time from the beginning of the mission sec
Time h Time from the beginning of the mission hour
Pu Ullage pressure kPa
Plsat Liquid saturation pressure corresponds to liquid temperature kPa
Tu Ullage temperature K
Tl Liquid temperature K
Tlw Temperature of the wall interfacing with liquid K
Tw Temperature of the wall interfacing with ullage K
Tint Liquid - ullage interface temperature K
Pv Vapor pressure in ullage kPa
Phe Non-condensable gas pressure in ullage kPa
Bond Bond number by tank diameter –
Ml Bulk liquid mass kg
Ml lost Total mass lost by chill-down, liquid venting and through ullage-liquid interface kg
Fill Volumetric tank fill level %
Del Fill Sw Difference between initial liquid volumetric fill level and liquid level swell %
Liq Lev Bulk liquid level m
Del Lev Sw Difference between initial liquid level and liquid level swell m
Mv Mass of vapor in ullage kg
Mhe Ul Mass of non-condensable in ullage kg
Mhe Lost Non-condensable gas lost mass kg
Md Evap Evaporation-condensation mass flow rate through ullage liquid interface kg/sec
Md Boff Bulk liquid boil-off mass flow rate kg/sec
Md Cond Ullage condensation mass flow rate kg/sec
Mlost int Liquid total mass lost through ullage-liquid interface (evaporation and boil-off) kg
M Drop U Mass of liquid droplets in ullage (from ullage condensation and mixing kg
Err_Etot Total energy balance error compare to the external heat loads %
Err_Mtot Total fluid mass balance error compare to the total fluid mass %
Err_Mds Spayed mass flow rate error compare to the flow rate in manifold %

117
 Description of TankSIM Input and Output Files Variables (Continued)
Addition to Common print for all regimes excluding 'Slf'
Mix.Cycl Number of cycles during mixing (no TVS) or ullage wenting –
Hex.Cycl Number of cycles during TVS working or ullage venting by liquid saturation pressure –
Dutypu Duty cycle for ullage pressure controlling %
Dutypl Duty cycle for liquid saturation pressure controlling addition to all regimes for Ullage Venting %
Md Uvent Ullage venting mass flow rate kg/sec
M Uv tot Mass venting from the ullage kg
Time Uven Total ullage venting time sec
Addition to all regimes for Spray bar or Axial jet
Md Pump Pump mass flow rate kg/sec
Qd Pump Pump volumetric flow rate GPM
Dp Pump Pressure increase at the pump kPa
Md Lvent Liquid venting mass flow rate kg/sec
Qman tot Power eliminating from the manifold W
Tman End Temperature at the manifold end K
Pman End Pressure at the manifold end kPa
Qex tot Power eliminating by heat exchanger W
Lefex/Xsef Relative length where liquid fully evaporated in the heat exchanger, or vapor quality if it does not fully evaporated –
Tex End Temperature at the heat exchanger end K
Pex End Pressure at the heat exchanger end kPa
Dt J-T Temperature drop at the Joule-Thomson device K
Dp J-T Pressure drop at the Joule-Thomson device kPa
X J-T Vapor quality after the Joule-Thomson device –
Pback Back pressure at venting device kPa
Addition to all regimes for Liquefaction
Md Ull Ullage condensation flow rate kg/sec
Md Uwall Ullage wall condensation flow rate kg/sec
Md Vapor External vapor added to the tank flow rate kg/sec
Qwu Heat eliminated from the ullage W
Tsat-Twu Ullage wall subcooling compare to the vapor saturation temperature K

118
 REFERENCES

1. Alabovskii, N.: “Thermal Interaction at the Gas-Liquid Interface in an Apparatus with an

Immersed Burner,” Kiev Polytechnical Institute. Translated from Inzhenerno-Fizicheskii
Zhurnal, Vol. 22, No. 1, pp. 117–122, January 1972.

2. Delhaye, J.M.: “Jump conditions and entropy sources in two-phase systems. Local instant for-
mulation,” Int. J. Multiphase Flow, Vol. 1, Issue 3, pp. 395–409, 1974.

3. Meserole, J.S.; Jones, 0.S.; Brennan, S.M.; and Fortini, A.: “Mixing-Induced Ullage Condensa-
tion and Fluid Destratification,” AIAA-87-2018, 23rd Joint Propulsion Conference, San Diego,
CA, doi.org/10.2514/6.1987-2018, 1987.

4. Rohsenow, W.M.; Hartnett, J.P.; and Cho, Y.I.: “Handbook of Heat Transfer,” 3rd ed.,
McGraw-Hill, NY, 1998.

5. Churchill, S.W.; and Usagi, R.: “A General Expression for the Correlation of Rate of Transfer
and Other Phenomena,” AIChE Journal, Vol. 18, pp. 1121–1128, 1972.

6. Gerstmann, J.; and Griffith, P.: “The Effects of Surface Instabilities on Laminar Film Conden-
sation,” Technical Report No. 5050-36, MIT, 1965.

7. Gerstmann, J.; and Griffith, P.: “Laminar Film Condensation on the Underside of Horizontal
and Inclined Surfaces” Int. J. Heat Mass Transfer. Vol. 10, pp. 567–580, 1967.

8. Oosthuizen, P.H.; and Naylor, D.: An Introduction to Convective Heat Transfer Analysis,
McGraw-Hill, NY, 1999.

9. Goldstein, R.J.; Sparrow, E.M.; and Jones, D.C.: “Natural Convection Mass Transfer Adjacent
to Horizontal Plates,” Int. J. Heat Mass Transfer, Vol. 16, Issue 5, pp. 1025–1035, 1973.

10. Rich, B.R.: “An Investigation of Heat Transfer from an Inclined Flat Plate in Free Convec-
tion,” Trans. ASME, Vol. 75, pp. 489–499, 1953.

11.

12. NIST Standard Reference Data 23, “Reference Fluid Thermodynamic and Transport Proper-
ties–REFPROP,” Version 9.0, 2010.

13. Miller, R.W.: “Flow Measurement Engineering Handbook,” 3rd ed., McGraw-Hill, NY, 1996.

119
14. Spalding, D.B.; and Cole, E.H.: Engineering Thermodynamics, 3rd ed., SI Units, Edward
Arnold Ltd, London, 1973.

15. Ward-Smith, A.J.: “Critical Flowmetering: The Characteristics of Cylindrical Nozzles With
Sharp Upstream eEdges,” Int. J. Heat & Fluid Flow, Vol. 1, No. 3, pp. 123–132, 1979.

16. Simmons, F.S.: “Analytical Determination of the Discharge Coefficient of Flow Nozzles,”
NACA TN 3447, April 1955.

17. Cole, R.; and Rohsenow, W.M.: “Correlation of bubble departure diameters for boiling of
saturated liquids,” Chem. Eng. Prog. Symp. Ser., Vol. 65, No. 92, pp. 211–213, 1968.

18. Levich, V.G.: Physicochemical Hydrodynamics, Prentice-Hall Inc., Englewood Cliffs, NJ, 1962.

19. Lockhart, R.W.; and Martinelli, R.C.: “Proposed Correlation of Data for Isothermal Two-
Phase, Two-Component Flow in Pipes,” Chem. Eng. Progress, Vol. 45, Issue 1, pp. 39–48, 1949.

20. Martinelli, R.C.; and Nelson, D.B.: “Prediction of Pressure Drop During Forced-Circulation
Boiling of Water,” Trans. ASME, Vol. 70, pp. 695–702, 1948.

21. Changhong, P.; Yun, G.; Suizheng, Q.; et al.: “Two-phase low and boiling heat transfer in two
vertical narrow annuli,” Nuclear Engineering and Design, Vol. 235, pp. 1737–1747, 2005.

22. Carey, V.P.: Liquid Vapor Phase Change Phenomena, Tailor & Francis, Hebron, KY, 1992.

23. Ghabari, A.; Farshad, F.F.; and Rieke, H.H.: “Newly developed friction factor correlation for
pipe flow,” J. Chem. Eng. Mater. Sci., Vol. 2, Issue 6, pp. 83–86, 2011.

24. Chen, J.C.: “Correlation for boiling heat transfer to saturated fluids in convective flow,” Ind.
Eng. Chem. Proc. Design and Dev, Vol. 5, No. 5, pp. 322–339, 1996.

25. Collier, J.G.: “Forced convective boiling,” Two-Phase Flow and Heat Transfer in the Power
and Process Industries, A.E. Bergles, J.G. Collier, J.M. Delhaye, G.F. Hewitt, and F. Mainger,
eds., Hemisphere, NY, 1981.

26. Idelchik, I.E.: Handbook of Hydraulic Resistance, 3rd ed., Jaico Publishing House, Mumbai,
2008.

27. Pappel, S.S.; Saiyed, N.H.; and Nyland, T.W.: “Acquisition and Correlation of Cryogenic
Nitrogen Mass Flow Data through a Multiple Orifice Joule-Thomson Device,” NASA
TM-103121, 1990.

28. Pappel, S.S.; Nyland, T.W.; and Saiyed, N.H.: “Liquid Hydrogen Mass Flow Through
a Multiple Orifice Joule-Thomson Device,” NASA TM–105583, AIAA-92-2881, 1992.

120
29. Jurns, J.M.: “Visco Jet Joule-Thomson Device Characterization Tests in Liquid Methane,”
NASA/CR—2009–215497, NASA Glenn Research Center, Cleveland, OH, March 2009.

30. The Lee Company, Technical Hydraulic Handbook, 11th ed., 2009.

31.

32. Ranz, W.E.; and Marshall, Jr., W.R.: “Evaporation From Drops,” Chemical Engineering
Progress, Vol. 48, No. 3, pp. 141–146, pp. 173–180, 1952.

33. Pasandideh-Fard, M.; Aziz, S.D.; Chandra, S.; and Mostaghimi, J.: “Cooling Effectiveness
of a Water Drop Impinging on a Hot Surface,” Int. J. Heat and Fluid Flow, Vol. 22, Issue 2,
pp. 201–210, 2001.

34. Castanet, G.; Liénart, T.; and Lemoine, F.: ”Dynamics and temperature of droplets impacting
onto a heated wall,” Int. J. Heat and Mass Transfer, Vol. 52, Issues 3–4 , pp. 670–679, 2009.

35. Baumeister, K.J.: and Simon, F.F.: “Leidenfrost Temperature—Its Correlation for Liquid
Metals, Cryogens, Hydrocarbons, and Water,” ASME J. Heat Trans., Vol. 96, No. 2,
pp. 166–173, 1973.

36. Chun, K.R.; and Seban, R.A.: “Heat Transfer to Evaporating Liquid Films,” J. Heat Transfer,
Vol. 93, No. 4, pp. 391–396, 1971.

37. Sonin, A.A.; Shimko, M.A.; and Chun, J.H.: “Vapor Condensation onto a Turbulent Liquid-
I. The Steady State Condensation Rate as a Function of Liquid-Side Turbulence,” Int. J. Heat
Mass Trans., Vol. 29, No. 9, pp. 1319–1332,1986.

38. Thomas, R.M.: “Condensation of Steam on Water in Turbulent Motion,” Int. J. Multiphase
Flow, Vol. 5, No. 1, pp. 1–15, 1979.

39. Brown, J.S.; Khoo, B.C.; and Sonin, A.A.: “Rate Correlation for Condensation of Pure Vapor
on Turbulent, Subcooled Liquid,” Int. J. Heat Mass Trans., Vol. 33, No. 9, pp. 2001–2018,
1990.

40. Lin, C.S.; and Hasan, M.M.: “Effect of Liquid Surface Turbulent Motion on the Vapor
Condensation in a Mixing Tank,” Proceedings of the 4th International Symposium on Transport
Phenomena in Heat and Mass Transfer, Sidney, Australia, July 14–19, 1991, Transport Phenom-
ena in Heat and Mass Transfer, J.A. Reizes, ed., pp. 1526–1537, 1992.

41 Chen, Y.; Schaeffer, B.M.; and Weislogel, M.M.: “Surface Evolver - Fluid Interface Tool, Ver-
sion 1.060. User’s Manual,” Department of Mechanical and Materials Engineering, Portland
State University, May 20, 2011.

121
42. Teertstra, P.; Yovanovich, M.M.; and Culham, J.R.: “Conduction Shape Factor Models for
Three-Dimensional Enclosures,” J. Thermophysics and Heat Transfer, Vol. 19, No. 4, October–
December, 2005.

43. Panton, R.L.: Incompressible Flow, 4th ed., John Wiley & Sons, NY, 2013.

44. Elenbaas, W.: “The Dissipation of Heat by Free Convection. The Inner Surface of Vertical
Tubes of Different Shapes of Gross-Section,” Physica IX, No. 8, pp. 865–874, September 1942.

45. Schmidt, E.; and Beckmann, W.: Das Temperatur- und Geschwindigkeitfeld vor einer Warme
abgebenden senkrechten Plate bei natúrlicher Konvektion,” Tech. Mech. und Thermodynamik,
Bd.1, Nr. 10, Okt.1930, pp. 341–349 and Bd. 1, Nr.11, pp. 391–406 (in German), November
1930.

46. Incropera, F.P.; and DeWitt, D.P.: Fundamental of Heat and Mass Transfer, 4th ed., John Wiley
& Sons, NY, 1996.

47. Rohsenow, W.M.: “Heat Transfer and Temperature Distribution in Laminar Film Condensa-
tion,” Trans. ASME, Vol. 78, pp. 1645–1648, 1956.

48. Tollmien, W.: “Calculation of Turbulent Expansion Processes,” NACA TM–1085, Langley,
1945.

49. Dominick, S.M.: “Mixing Induced Condensation Inside Propellant Tanks,” AIAA-84-
0514, 22nd Aerospace Sciences Meeting, Aerospace Sciences Meetings, Reno, NV, doi.
org/10.2514/6.1984-514 1984.

50. Patil, R.K.; Shende, B.W.; and Ghosh, P.K.: “Designing a Helical-Coil Heat Exchanger,”
Chemical Engineering, pp. 85–88, December 13, 1982.

122
123
 Form Approved
REPORT DOCUMENTATION PAGE OMB No. 0704-0188
The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining
the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing
this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operation and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302.
Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid
OMB control number.
PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY) 2. REPORT TYPE 3. DATES COVERED (From - To)

01–04–2017 Technical Memorandum
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER

Tank System Integrated Model: A Cryogenic Tank Performance 5b. GRANT NUMBER

Prediction Program
5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S) 5d. PROJECT NUMBER

L.G. Bolshinskiy,* A. Hedayat, L.J. Hastings (retired), S.G. Sutherlin, 5e. TASK NUMBER

A.R. Schnell, and J.P Moder**

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION

REPORT NUMBER
George C. Marshall Space Flight Center
Huntsville, AL 35812 M–1431
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING/MONITOR’S ACRONYM(S)

National Aeronautics and Space Administration NASA

Washington, DC 20546–0001 11. SPONSORING/MONITORING REPORT NUMBER

NASA/TM—2017–218239
12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified-Unlimited
Subject Category 20
Availability: NASA STI Information Desk (757–864–9658)
13. SUPPLEMENTARY NOTES

Prepared by the Propulsion Systems Department, Engineering Directorate, *Jacobs ESSSA Group
(contract NNM12AA41C), Huntsville, AL, **NASA Glenn Research Center, Cleveland, OH
14. ABSTRACT

Accurate predictions of the thermodynamic state of the cryogenic propellants, pressurization rate, and performance
of pressure control techniques in cryogenic tanks are required for development of cryogenic fluid long-duration
storage technology and planning for future space exploration missions. This Technical Memorandum (TM) presents
the analytical tool, Tank System Integrated Model (TankSIM), which can be used for modeling pressure control and
predicting the behavior of cryogenic propellant for long-term storage for future space missions. Utilizing TankSIM,
the following processes can be modeled: tank self-pressurization, boiloff, ullage venting, mixing, and condensa-
tion on the tank wall. This TM also includes comparisons of TankSIM program predictions with the test data and
examples of multiphase mission calculations.

15. SUBJECT TERMS

modeling on-orbit liquid cryogens storage and pressure control, cryogenic propellants, in-space cryogenic
propulsion, heat and mass transfer in cryogenic liquids
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF ABSTRACT 18. NUMBER OF 19a. NAME OF RESPONSIBLE PERSON
a. REPORT b. ABSTRACT c. THIS PAGE PAGES STI Help Desk at email: help@sti.nasa.gov

U U U UU 144 19b. TELEPHONE NUMBER (Include area code)

STI Help Desk at: 757–864–9658
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std. Z39-18

124
National Aeronautics and NASA/TM—2017–
Space Administration
IS02
George C. Marshall Space Flight Center
Huntsville, Alabama 35812

Tank System Integrated Model: A Cryogenic

Tank Performance Prediction Program
L.G. Bolshinskiy
Jacobs ESSSA Group, Huntsville, Alabama

A. Hedayat, L.J. Hastings (retired), S.G. Sutherlin, and A.R. Schnell

Marshall Space Flight Center, Huntsville, Alabama

J.P. Moder
NASA Glenn Research Center, Cleveland, Ohio

January 2017