Open access peer-reviewed chapter

Two-Phase Flow Modeling of Cryogenic Loading Operations

Written By

Ekaterina Ponizovskaya-Devine, Dmitry Luchinsky, Michael Foygel, Vasil Hafiychuk, Michae Khasin, Jared Sass and Barbara Brown

Submitted: 30 October 2018 Reviewed: 23 January 2019 Published: 15 March 2019

DOI: 10.5772/intechopen.84662

From the Edited Volume

Heat and Mass Transfer - Advances in Science and Technology Applications

Edited by Alfredo Iranzo

Chapter metrics overview

1,246 Chapter Downloads

View Full Metrics


We consider problem of modeling and controlling two-phase cryogenic flows during ground loading operations. We introduce homogeneous moving front and separated two-phase flow solvers that are capable of fast and accurate online predictions of flow dynamics during chilldown and transfer under nominal conditions and in the presence of faults. Concise sets of cryogenic correlations are proposed in each case. We present results of application of proposed solvers to the analysis of chilldown in large-scale experimental cryogenic transfer line build in Kennedy Space Center. We discuss optimization of parameters of cryogenic models obtained using general inferential framework and an application of the solvers to the fault detection and evaluation based on D-matrix approach. It is shown that solver’s predictions are in good agreement with experimental data obtained for liquid nitrogen flow in nominal regime and in the presence of faults.


  • cryogenic flow
  • chilldown
  • moving front
  • two-phase flow
  • optimization
  • flow boiling
  • correlations
  • heat transfer

1. Introduction

Deep space exploration requires development of autonomous management system [1] that can recognize, predict, and control two-phase cryogenic flows online without human interaction.

This long-standing problem is of great complexity because flowing fluids are often far away from thermal equilibrium (e.g., during chilldown), and their dynamics involves heat and mass transfer [2, 3, 4]. A sketch of the two-phase flow is shown in the Figure 1. Within the accepted approximation [3, 4], the pressure p in the system is common. The temperature of the wall T w , vapor T g , and liquid T l can all be different. Gas and liquid velocities can have different values and directions. The heat and mass exchange and the drag take place both at the wall and at the interface. The shape of the interface can be complex [5].

Figure 1.

Sketch of two-phase flow. α is the void fraction, A is the cross-section of the flow, l g and l l are dry and wetted perimeters, and l i is interface perimeter. Other parameters are explained in the text.

During the past decades, a number of efficient algorithms have been developed for analysis of multiphase flows [3, 4, 6, 7, 8, 9]. However, the problem still presents a major challenge to the scientific community due to the lack of a general agreement regarding the fundamental physical models describing two-phase flows [10] and relatively poor knowledge of heat and mass transfer correlations in boiling cryogenic flows [11, 12, 13, 14, 15].

To address these issues critical issues, NASA has implemented an impressive program of research [1, 16, 17, 18] that resulted in emergence of space-based fluid management technologies; see, e.g. [14, 15, 19, 20, 21, 22, 23, 24] and references therein. Specifically, two-phase separated flow models were developed for some of the flow regimes [25, 26]. A number of optimization techniques have become commercially available for analysis of the model parameters and data correlations [27].

Despite these efforts small-time steps and instabilities [26] or implicitness of numerical scheme [28, 29] impose substantial limitations on the speed of the solution and possibility of online application. And an extensive research is currently under way to improve accuracy of predictions of transient cryogenic flows [15, 19, 24, 27].

Some of the grand challenges of this research include development of fast algorithms, choice of cryogenic correlations, and applications of the algorithms to the full-scale practical systems.

In this paper, we report on the development of fast two-fluid models suitable for online analysis of cryogenic flows and discuss applications of these models to the large-scale experimental prototypes of cryogenic transfer lines for ground launch systems.

The paper is organized as follows. In the next section, we describe the moving front model and algorithm of integration of homogeneous cryogenic flow. In Section 3, we discuss the separated two-phase flow model. In Section 4 we briefly discuss applications of these models to the inference of cryogenic correlations and fault detection in cryogenic transfer line. Finally, in the Conclusions, we summarize the obtained results and discuss directions of future work.


2. Homogeneous moving front model

2.1 Model equations

One of the simplest models that can recognize and predict cryogenic two-phase fluid dynamics online in nominal and off-nominal (i.e., in the presence of faults) flow regimes is moving front homogeneous model [30, 31]. The model was developed for the integrated health management system for cryogenic loading operations at the Kennedy Space Center (KSC). It is similar to the model [32] that describes transient phenomena in evaporators and condensers in refrigeration. It accounts for the heat exchange with the pipe walls and for the moving liquid front in chilldown regime. It reproduces with good accuracy the time evolution of the pressure, temperature, and flow rate in spatially distributed cryogenic systems in the nominal loading regime. It can detect and isolate major faults of cryogenic loading operation in the real time. It works much faster than commercially available full-scale fluid dynamics solvers and can be used for preliminary analysis of the loading protocols on the ground and under microgravity conditions [33].

Within this model the cryogenic transfer line is divided into a set of control volumes with constant cross-section shown in Figure 2. The long pipes can consist of several control volumes. Following [30] we will assume that within each control volume flow is one-dimensional along the tube’s axis in z-direction. The change in pipe diameter considered only at the control volume boundaries. We further neglect the viscous dissipation and spatial variations of the pressure. The model is based on the conservation equations for the mass, momentum, and specific energy e in one dimension:

Figure 2.

(a) Three main control volumes of the model. Subindexes 1, 2, and 3 refer to (1) subcooled liquid region; (2) saturated two-phase region; and (3) superheated vapor region. (b) The convention for the nodes and junction notations: K , L and M refer to the centers of control volumes and index j enumerates junctions between control volumes.

ρ , t + ρu , x = 0 ρu , t + u ρu , x + p x + 1 A τ w l w + ρg sin θ = 0 ρe , t + ρuH , x + 1 A q ̇ w l w = 0 E1

Here H is enthalpy, ρ is density, and u is velocity of the flow. τ w is the friction losses, l w is the pipe length, g is the gravity, θ is the angle of the pipe with horizontal plane, q ̇ w is the heat transfer from the pipe wall to the fluid mixture, and A is the cross-section area of the pipe.

The wall temperature is described by the energy conservation equation:

c w ρ w A w T w , t = hπD T T w + h amb πD T amd T w + κ A w T , xx E2

where A w is the cross-sectional area of the pipe wall with diameter D , which is in thermal contact with the ambient fluid of temperature T amb and with the fluid of temperature T flowing in the pipe. c w is the specific heat, ρ w is density, and κ is thermal conductivity of the wall material. h is the heat transfer coefficient at the inner wall surface, and h amb is the heat transfer coefficient at the ambient side.

Two-phase flow is considered to be homogeneous with void fraction of the vapor equal α and 1 α for the liquid. Each volume can contain pure liquid with α = 0 , pure gas with α = 1 , and two-phase flow, and there are two volumes that contain the moving front from liquid to two-phase and from two-phase to gas. The local void fraction and the vapor mass fraction (fluid quality x ) in the two-phase region are

α = A v A = A v A v + A l and x = M v M v + M l = 1 + 1 α ρ v α / ρ l S 1 E3

where S = u v / u l is the slip ratio. It can be shown that for the steady flow without slip ( S = 1 ), the local density and specific enthalpy of the mixture are [34]:

ρ = αρ v + 1 α ρ l ; H = xH v + 1 x H l E4

The momentum equation is reduced to a set of quasi-static algebraic equations assuming that the flow is relatively slow and there are no shock waves:

u 2 , x + 1 A τ w l w tp = p , x ρgsin θ E5

The model equations are closed by adding equation of state ρ g , l = ρ g , l p H v , l in the form of tables and by providing a set of cryogenic correlations discussed below.

There is a multitude of data on void fraction-quality correlations for different flow patterns in horizontal and upward inclined pipes [30]. Here we use correlation of the form [30]

α x p = 1 + a p 1 x 1 1 E6

where the factor a is

a = ρ v ρ l n 1 μ l μ v n 2 E7

and μ is the viscosity of the fluid. Coefficients n 1 and n 2 can have different values, e.g., according to Thome [35] n 1 = 0.89 , n 2 = 0.18 , while it is suggested to use n 1 = 0.67 , n 2 = 0 in [34].

The estimated average void fraction α can be used to evaluate the heat transfer coefficient from the tube wall to the TP mixture in the TP region n = 2 as h 2 = α h v + 1 α h l , where h v , l = κ D h Nu v , l , Nu is Nusselt number, and Pr is Prandtl number.

For the fully developed laminar flow, Nu = 3.66 . For the turbulent flow Re > 2300 , one can use the modified Dittus-Boelter approximation with the Nusselt number calculated separately for each of the vapor and liquid components. For example, according to [34]

Nu = 0.0214 Re 4 / 5 100 Pr 2 / 5 1 + D h / L 2 2 / 3 E8

for 0.5 < Pr < 1.5 and

Nu = 0.012 Re 0.87 280 Pr 2 / 5 1 + D h / L 2 2 / 3 E9

for 1.5 < Pr < 500 . Here the effective Reynolds number for each of the components, Re v , l = m ̇ D h μA , is estimated by means of the average hydraulic diameters D h , v = αD and D h , l = 1 α D .

For the forced convection, the heat transfer coefficient [9] at the environmental side is

h amb = 2 κ πD Nu amb , Nu amb = Nu lam + Nu turb , E10

where laminar and turbulent contributions to the Nusselt number

Nu lam = 0.0664 Re 1 / 2 Pr 1 / 3 and Nu turb = 0.037 Re 0.8 Pr 1 + 2.443 Re 0.1 Pr 2 / 3 1 E11

are calculated for the ambient fluid.

The specific form of the linearized conservation equations depends on the type of the flow. Using Taylor expansion for density ρ in terms of pressure p and enthalpy H , we can rewrite the mass and energy conservation equations in terms of pressure p and enthalpy H or pressure p and void fraction α . In the case of two-phase control volume, the equations are (with β = 1 α ):

α ρ v p + β ρ l p dp dt + ρ v ρ l dt = m ̇ L m ̇ R LA α ρH v p + β ρH l p 1 dp dt + ρ v H v ρ l H l dt = m ̇ L H L m ̇ R H R LA + 4 h T w T D E12

and in the case of single-phase control volume (e.g., vapor), we have

ρ p dp dt + ρ H dH dt = m ̇ L m ̇ R LA ρH p 1 dp dt + ρH H dH dt = m ̇ L H L m ̇ R H R LA + 4 h T w T D E13

Here subindexes R and L refer to the right and left boundaries; ρ v and ρ l are the saturated density for vapor and liquid in Eq. (12), while in Eq. (13) for the single-phase, ρ is the density of vapor.

The case when the saturated two-phase flow and gas flow coexist in one control volume integration of the conservation equations that accounts for the moving front position L 2 results in the following set of equations:

α ρ v P + β ρ l p dp dt + ρ v ρ l dt β dln L 2 dt = m ̇ L m ̇ R L 2 A α ρ v H v p + β ρ l H l p 1 dp dt + ρ v H v ρ l H l dt β dln L 2 dt = m ̇ L H L m ̇ R H v L 2 A + 4 h D T w T ρ p dp dt + ρ H dH dt + ρ ρ v dln L 3 dt ) = m ̇ L m ̇ R L 3 A ρH p 1 dp dt + ρH H dH dt + ρH ρ v H v dln L 3 dt = m ̇ L H l m ̇ R H R L 3 A + 4 h D T w T E14

Here L 3 = L L 2 and L is the length of the control volume; see Figure 2. These equations are solved simultaneously to find pressure, enthalpy, and density in each region and the position of the moving front [30]. The wall temperature is calculated for each region according to Eq. (2).

The mass fraction x is used to calculate the enthalpy flux at the boundaries, and it is assumed to change linearly with coordinate in the phase transition region. The new enthalpy for two-phase flow depends on x , and the new density is defined by the void fraction α in Eq. (4), while the relation between x and α is defined by Eq. (6).

The model was applied to consider two-phase flow of liquid nitrogen in the KSC experimental cryogenic transfer line. It takes into consideration compressibility of the vapor and condensation-evaporation effects in the flowing cryogenic fluid through the initially hot transfer line. Both the transient and steady-state regimes of flow can be described by accounting for the viscous and thermal interaction between the two-phase fluid and the pipe wall. Similarly, the fluid dynamics in the storage and vehicle tanks includes evaporation at the liquid surface and the heat exchange at the tanks walls as described in [36].

2.2 Algorithm

The conservation equations (1) are solved for the set of control volumes that form a staggered grid. The solution proceeds by iterating the following steps. For each volume we calculate pressure P , position of the front L , void fraction α , vapor ( v ) and liquid ( l ) enthalpies H v l , and two wall temperatures before and after the front T w 1 and T w 2 . These variables are calculated at the center of control volume.

Next the mass m ̇ and enthalpy fluxes H m ̇ are calculated at the boundaries of each control volume; see Figure 2, using the following equation:

m ̇ j 2 A n 2 1 ρ n f n L n D n + K n + x 2 γρ g + 1 x 2 1 γ ρ l n = Δ p j Δ z j ρ j gsinθ E15

The subindex n for frictional factor f and minor losses K refers to the control volume on the left- or right-hand side of the junction. The stability of this approximation was enforced by the donor-like formulation of the frictional losses for the mass fluxes [35].

For each control volume k , the continuity condition is added as follows: m ̇ k , in = m ̇ k 1 , out for the serial connections of the k 1 , k , and m ̇ k 1 , out = m ̇ k 1 , in + m ̇ k 2 , in if the k -th control volume consists of two parallel control volumes k 1 and k 2 .

Updated values of the thermodynamic variables, enthalpy, and mass fluxes are used to calculate new parameters of friction and heat transfer correlations.

The algorithm can be briefly summarized as follows: (i) calculate variables P , H , L , α , and T w at the center of the control volume using mass and heat fluxes from the previous step, (ii) update location of the moving fronts, (iii) calculate new mass and heat fluxes at the boundaries of each control volume, and (iv) update friction and heat transfer correlations and return to the first step.

The moving front model allows to account for multiple faults in cryogenic transfer lines including, e.g., valve clogging, valve stuck open/closed, and heat and mass leaks [37, 38]. However, the homogeneous moving front model has some important limitations. It does not include shock waves and cannot distinguish nonhomogeneous ( u g u l ) and non-equilibrium ( T g T l ) flows such as e.g. counter flows and stratified flows. To overcome these limitations we developed two-phase separated solver for cryogenic flows discussed below.


3. Two-phase flow model

Having in mind fast online applications of the solver we have tested a number of algorithms for the separated two-phase flow [30, 31, 39, 40, 41, 42, 43, 44]. It was shown that the nearly implicit algorithm (NIA), similar to one developed in [7], can be applied successfully for online predictions of nonhomogeneous and non-equilibrium flows. Below we describe briefly the NIA and the corresponding model equations. The details can be found in [39, 40, 42, 43, 44, 45].

In the nearly implicit algorithm, the six-equation model consists of a set of conservation laws for the mass, momentum, and energy of the gas [6, 7, 39]

A αρ g , t + A αρ g u g , x = A Γ g A αρ g u g , t + A αρ g u g 2 , x + p , x = A Γ g u ig Aαg ρ g z , x τ gw l wg τ gi l i A αρ g E g , t + A αρ g E g u g , x = A Γ g H g Ap α , t pAα u g , x + q ̇ gw l wg + q ̇ gi l i E16

and liquid phases

A βρ l , t + A βρ l u l , x = A Γ g A βρ l u l , t + A βρ l u l 2 , x + p , x = Aβg ρ l z , x τ lw l wl τ li l i A Γ g u il E l ρ l , t + E l ρ l u l , x = Ap β , t pAβ u l , x + q ̇ lw l wl + q ̇ li l i A Γ g H l . E17

Variables in Eqs. (16) and (17) are defined following convention of Eq. (1); the details are provided in the nomenclature. Importantly, velocities and temperatures of the gas and liquid phases are now independent variables; the mass flux and the drag are now defined both at the wall and at the interface.

The equations above are coupled to the equation for the wall temperature T w :

ρ w c w d w T w t = h wg T g T w + h wl T l T w + h amb T amb T w . E18

In the model (16), (17) is non-hyperbolic [46]. It does not have a complete set of real eigenvalues and does not represent a well-posed system of equations [39, 47]. It also lacks positivity and possesses instabilities due to phase appearance/disappearance process [8, 48].

Despite these difficulties it was successfully employed to predict two-phase flows of boiling water in large-scale system [6, 7]. In our developments we were following the guidelines of this research.

3.1 Algorithm

The algorithm originates from [49, 50] all-speed implicit continuous-fluid Eulerian algorithm. In our work we closely follow recent modifications [6, 7] that enhance stability and eliminate material CFL restrictions using two-step nearly implicit schemes. Within this approach Eqs. (16) and (17) are expanded and solved [39, 42] for velocities and pressure in the first step, and the unexpanded conservation equations are solved for densities and energies in the second step of the algorithm.

The near implicitness of the algorithm is determined by the fact that nearly all terms in the numerical scheme are taken at the new time step. Yet the algorithm is very efficient because linearized form of equations is used, which allows to use block tridiagonal solver for velocity in the first step and to reduce the second step to the solution of four independent tridiagonal matrices.

The stability of the code was further enhanced by using upwinding and staggered grid methods as well as ad hoc smoothing and multiple time step control techniques [39]. The non-hyperbolicity was suppressed using so-called virtual mass term [7]. The time step was controlled to keep all the thermodynamic variables within the predetermined limits, to ensure that the change of these variables at any given time step does not exceed 25% of their values obtained at the previous time step and to enforce mass conservation in each control volume and in the system as whole. In addition, smooth transition between phases [39, 51] was used to adjust temperature, velocity, and density.

The models (15)(17) have to be completed with the equations of state and the constitutive relations. The equations of state were implemented using tables for each phase. The constitutive relations are briefly summarized below.

3.2 Constitutive relations

Following the algorithm outlined above, we calculate new velocities, pressure, densities, and energies for both phases using Eqs. (15)(17). Once these variables are calculated, we have to determine the flow pattern and boiling regime in the system and update drag, heat, and mass fluxes using appropriate correlations. First, the flow patterns have to be recognized by the code. There are numerous flow pattern correlations, but their validation for cryogenic flows remains an open problem [11, 12, 45]. Here we considered only a few flow regimes following earlier work by Wojtan et al. [52] including stratified-Wavy-to-stratified transition, stratified-Wavy-to-annular-intermittent transition, and dry-out transition. We approximated the location of these boundaries by low-dimensional polynomials and used polynomial coefficients as fitting parameters. The details of the transformation can be found in [45].

Next, the heat transfer and pressure losses for a given flow pattern have to be calculated. There are literally hundreds of heat transfer correlations proposed in the literature [53]. A concise subset of these correlations briefly discussed below was selected as a result of extensive numerical analysis and validation using two experimental setups; see Section 4 and Refs. [44, 45, 54].

3.2.1 Heat and mass transfer

As was mentioned above, the total mass transfer in the two-phase model is the sum Γ g = Γ wg + Γ ig of the mass transfer at the wall and at the interface:

Γ wg = q ̇ wl H g H l ; Γ ig = q ̇ li + q ̇ gi H g H l ; E19

where H g H l is the difference between saturated gas and liquid enthalpies H g , s H l for positive flux Γ > 0 while it is H g H l , s for negative flux.

The heat transfer coefficients h for the gas (subscript g ) and liquid ( l ) at the wall ( w ) and at the interface ( i ) are defined by the following relations:

q ̇ wg = h wg T w T g ; q ̇ ig = h ig T l , s T g ; q ̇ wl = h wl T w T l ; q ̇ il = h il T l , s T l . E20

For relatively low mass fluxes considered in this work ( < 1 metric ton per square meter per second), the heat transfer correlations can be obtained as the mass flux and void fraction dependent corrections to the pool boiling correlations [55, 56, 57, 58, 59]. Accordingly, the following pool boiling mechanisms of the heat transfer were modeled: convective heat transfer, nucleate, transition, and film boiling.

Convective heat transfer coefficient h cb was calculated for four flow regimes [6, 60]. The value of D h h cb / κ was taken to be 4.36 in the forced laminar regime, while in forced turbulent regime, it was calculated as 0.023 Re 0.8 Pr 0.4 [53]. The coefficient D h h cb / κ for natural convection was calculated as 0.1 Gr Pr 1 / 3 for laminar and as 0.59 Gr Pr 1 / 4 for turbulent flows [61]. Here

Pr = μ C p κ and Gr = ρ 2 g β T T w T l g D 3 μ 2 E21

are Prandtl and Grashof numbers, respectively, β T is the coefficient of thermal expansion, and D h is the hydraulic diameter. The maximum value of h cb calculated for four regimes is taken as the value for the convective heat transfer.

The next critical temperature that defines the shape of the boiling curve and influences the chilldown corresponds to the onset of nucleation boiling T onb .

The onset of nucleate boiling correlations is based on the analysis of the balance between mechanical and thermodynamical equilibrium [3]. The corresponding T onb and heat flux q ̇ onb are [62, 63, 64].

T onb = T s + F 1 + 1 + 2 Δ T sub F , q ̇ onb = B Pr 2 Δ T sat 2 = h cb T onb T l E22

where B = ρ g h lg κ l 8 σ T s , F = h cb Pr l 2 2 B , Δ T sat = T onb T s is the wall superheat, and Δ T sub = T s T l is liquid subcooling temperature.

When the wall superheat exceeds Δ T sat = T onb T s , the nucleation boiling begins, and the heat flux to the wall may increase by more than an order of magnitude. This increase continues until the heat flux approaches its critical value q ̇ chf .

The critical heat flux q ̇ chf and the critical wall superheat T chf correlations are crucial for predicting chilldown and dry-out phenomena in cryogenic flows. Yet the corresponding experimental data remains sparse, and values of q ̇ chf and T chf are often estimated using mechanistic models [3, 15, 56, 65].

For the pool boiling value of q ̇ chf , 0 , we used Zuber [58] correlation

q ̇ chf , 0 = π 24 h lg ρ g σg ρ l ρ g ρ g 2 1 / 4 ρ l ρ l + ρ g 1 / 2 . E23

In boiling flows further corrections have to be introduced to take into account the dependence of the heat flux on the void fraction, velocity, and subcooling of the flow [55, 56].

q ̇ chf = q ̇ chf , 0 α cr α 1 + a 1 ρ l c l Δ T sub ρ g h lg + a 2 Re l + a 3 Re l ρ l c l Δ T sub ρ g h lg , E24

where α cr is the critical value of the void fraction and a i are constants, e.g., a 1 = 0.0144 , a 2 = 10 6 , a 3 = 0.5 × 10 3 [57], and α cr = 0.96 [55] for water.

In practice, we often used a simpler expression; see for the details [55, 66].

The temperature T chf for the critical heat flux was estimated in this work using the approach proposed by Theler and Freis [67].

T chf = T s 1 T s R g h lg log 2 k g + 1 , E25

where k g is the isoentropic expansion factor that for ideal diatomic gases is 7/2 and R g is the specific gas constant. T chf was further corrected for boiling flows using Iloeje-type correlations [66] similar to the one applied Δ T mfb , see below.

When wall superheat exceeds Δ T chf = T chf T s , the transition boiling begins, and the heat flux to the wall decreases sharply as a function of the wall temperature until the latter reaches minimum film boiling temperature T mfb .

The minimum film boiling regime signifies a complete separation of the fluid flow from the wall by the vapor film. The corresponding value of the wall superheat Δ T mfb = T mfb T s was estimated by Berenson as [68, 69].

Δ T mfb , 0 = 0.127 ρ g h lg κ g × g ρ l ρ g ρ l + ρ g 2 / 3 σ g ρ l ρ g 1 / 2 μ g ρ l ρ g 1 / 3 E26

Iloeje [55, 66] has corrected Berenson equation to take into account the dependence of the Δ T mfb on the quality and mass flux of the boiling flows in the form

Δ T mfb = c 1 Δ T mfb , 0 1 c 2 χ e c 3 1 + c 4 G c 5 , E27

where χ e is the equilibrium quality, G is liquid mass flux, and a i are constants, e.g., a 1 = 0.0144 , a 2 = 10 6 , a 3 = 0.5 × 10 3 [57], and α cr = 0.96 [55] for water.

In this work the heat flux to the wall in the film boiling regime was taken in the form of Bromley correlations:

h f b = c 1 g ρ g κ g 2 ρ l ρ g h ˜ lg c pg D T w T spt Pr g 0.25 1 c 2 χ e c 3 1 + c 4 G c 5 E28

Typical values of the parameters used in simulations are the following: (i) c 1 = 2.0 , (ii) c 2 = 1.04 , (iii) c 3 = 2.0 , (iv) c 4 = 0.2 , and (v) c 5 = 0.1 .

The minimum film boiling heat flux can now be defined as

q ̇ mfb = h f b Δ T mf b . E29

In the region of single-phase gas flow, the heat transfer coefficient is given by h cb discussed above with appropriately modified parameters. Transition to the single-phase heat transfer is initiated when dry-out transition is detected.

Interpolation: The transitions between various flow boiling regimes are characterized by a number of critical points including (i) onset of nucleate boiling, (ii) critical heat flux, (iii) minimum film boiling, and (iv) onset of dry-out. The values of the heat flux between the first three critical points were interpolated as follows.

In the regime where the wall superheat is increasing from Δ T onb to Δ T chf , the heat flux was defined using the following interpolation [6]:

q ̇ nb = y n q ̇ onb + 1 y n q ̇ chf , E30

where n is constant, y is defined as T w T onb / T chf T onb , while q ̇ chf and q ̇ onb are given by Eqs. (24) and (22), respectively.

Similarly, for the regime where the wall superheat is increasing from Δ T chf to Δ T mfb , the transition boiling heat flux was interpolated as [6].

q ̇ tb = f tb q ̇ chf + 1 f tb q ̇ mfb , E31

where f tb = T w T mfb T chf T mfb 2 , while T chf , T mfb , and q ̇ mfb are given by Eqs. (25), (27), and (29), respectively.

Using parameterization of the heat transfer correlations described in Section 3.2, we could obtain smooth transition from the pool boiling correlations to the flow boiling regime as the mass flux and void fraction of the flow increase [44].

3.2.2 Pressure drop

The constitutive relations are completed by providing pressure drop correlations as follows. For the single-phase flow, the wall drag was calculated using standard relations τ wl = 0.5 f wl ρ l u l 2 and τ wg = 0.5 f wg ρ g u g 2 with the friction factors f wl g for turbulent and laminar flow given by Churchill approximation [70].

The two-phase pressure drop dp dz 2 ϕ is given by the Lockhart-Martinelli correlations [71]. The pressure losses are partitioned between the phases as follows [7]:

τ wg l wg = α g dp dz 2 ϕ 1 α g + α l Z 2 , τ wl l wl = α l dp dz 2 ϕ Z 2 α g + α l Z 2 E32

where Z 2 is given by

Z 2 = f wl Re l ρ l u l 2 α wl α l / f wg Re g ρ g u g 2 α wg α g , E33

Friction factors f wg l are the same as in above, while coefficients α wl and α wg depend on the flow pattern [7].

The interface drag is given by τ ig = τ il = 1 2 C D ρ g u g u l u g u l , where interfacial drag coefficient C D depends on the flow pattern [6].

We note that the functional form of the correlations adopted in this work is not unique and a number of alternative presentations can be used [2, 6, 45, 60].


4. Applications

The algorithms discussed above were used to develop a number of applications related to the management of cryogenic flows including, e.g., fault detection [72], diagnostics [38], evaluation [37], parameter inference, and optimization [44, 73].

To develop these applications and to validate our approach, we used experimental data obtained for cryogenic chilldown in two systems: (i) horizontal transfer line at National Bureau of Standards (currently NIST) [74] and (ii) large-scale experimental transfer line build at KSC [18].

The first system [74] is a 61 m long vacuum jacketed line with internal diameter of the copper pipe 3/4 inches. Four measurement stations were located along the pipe as shown in Figure 3. In this case the cold front continuously propagates through the initially hot pipe, and chilldown is achieved at approximately 30, 80, 100, and 130 s at each station. The model predictions for the fluid temperature are compared with the experimental data in the inset of the Figure 3. These results were obtained for subcooled flow and tank pressure 4.2 atm. A detailed comparison for various chilldown regimes is given in [44].

Figure 3.

Sketch of the cryogenic transfer line built at NIST. It includes storage tank (ST) and horizontal transfer line with four measurement stations. The model predictions for the temperature at four stations are compared with experimental data in the inset [44]. The experimental data are shown by solid black lines and open symbols. The model predictions are shown by colored solid lines.

A more complex cryogenic test-bed (CTB) system with multiple control and bleed valves was built at KSC; see Figure 4. It consists of a 6000 gallon storage tank (ST) connected to a 2000-gallon vehicle tank (VT), control (CV) and bleed (BV) valves, and pump that control the conditions of cryogen flow. The total length of the transfer line is about 45 m. The diameter of the stainless steel pipe varies from 0.1524 to 0.0254 m. A set of the temperature (TT) and pressure (PT) sensors allows for accurate detection of the flow conditions and gradual chilldown of the system. The liquid motion through the transfer line is driven by an elevated pressure in the storage tank.

Figure 4.

Sketch of the cryogenic transfer line built at KSC. It includes storage tank (ST) and vehicle tank (VT); the in-line control valves, CV1 through CV8; remotely controlled bleed valves, BV1 through BV8; and 10 in-line temperature sensors (TT) and 9 in-line pressure sensors (PT).

We note that both algorithms can be used for online fault detection and control of the flow. The homogeneous algorithm is about 10 times faster as compared to the NIA and can integrate 2 h of real-time loading in the CTB in less than 1 s. However, the NIA provides deeper insight into mechanism of two-phase flow and can distinguish a number of important non-equilibrium regimes including, e.g., counter flow and vertical stratification. We therefore proposed to use for online applications both algorithms in parallel to take advantage of speed of homogeneous model and greater precision of the two-phase flow model. Below we discuss some of the applications of these algorithms.

The validation and convergence of the algorithms using first horizontal cryogenic transfer line build in NIST were considered in details in our earlier work [44, 45, 54]. Here we briefly discuss some of the applications of the algorithms to the analysis of the flow in the CTB system.

4.1 Parameter inference and optimization

Practical application of the algorithms to management of specific large-scale cryogenic system involves adjustment of parameters of cryogenic flow boiling and system component correlations. In addition, parameter inference and optimization are required for continuous learning of multiple faults in the system (e.g., [45, 54] clogging, valves stuck open/closed, heat and mass leaks) and development of the recovery strategy. Parameter inference is important because not only functional form of the correlations is not unique, but also the values of numerical “constants” in these correlations should be considered as fitting parameters [44].

To address this problem, we developed inferential framework [44] that allows for simultaneous estimation and optimization of a large number of parameters of cryogenic models.

The application of this framework to the analysis of cryogenic correlations in NIST system was discussed in details in our earlier work [44]. Here we provide a brief example of parameter optimization during chilldown by applying two-phase flow model to CTB system [45, 54]; see Figure 4. In the first example we use direct search for preliminary estimation of some of the model parameters. The search is performed simultaneously for six parameters, which are scales for flow rate of input valve CV1, in-line valve CV2, dump valve BV3, dump valve BV2, dump valve BV2, and minimum film boiling heat transfer. The cost function is the sum of the square distances between experimental data and model predictions for four temperature sensors (TT202, TT105, TT162, and TT174) and four pressure sensors (PT102, PT104, PT161, PT173). It can be seen from Figure 5 that model predictions nearly converge to the sensor’s readings after 100 iterations.

Figure 5.

Convergence of the direct search algorithm for the model predictions (dashed lines) toward experimental data (solid black lines) measured for three temperature sensors: TT105, TT162, and TT174. The cost values for three dashed lines are highlighted by squares in the right panel.

To optimize and infer the full set of nearly 100 parameters of the KSC system, we used a number of steps. In the first step, we performed extensive sensitivity analysis to find parameters of the model that can significantly modify the dynamics of the system. At this step the number of optimization parameters can be reduced by the factor of ~2 by selecting only significant parameters for further optimization. Next, we inferred the obtained subset of sensitive parameters of the heat transfer and pressure correlations using NIST system [44] that does not have system components other than input valve. Next, we inferred parameters of the components of more complex CTB system using heat transfer and pressure correlations obtained at the previous step. The optimization itself was performed in two steps. First, we used direct search to roughly estimate parameter values close to global minimum, see example above and Figure 5. Then we refined the results using global stochastic optimization methods. The details of the approach are provided in [44, 54].

The comparison of the model predictions with the temperature sensor measurements obtained within this inferential framework is shown in Figure 6. It can be seen from the figure that by using this approach, we were able to accurately fit all the sensor’s measurements of the CTB system. We note that similar results were obtained using homogeneous model but for a fewer sensor’s measurements.

Figure 6.

Separated model predictions of the fluid temperature (green) during complete sequence of the CTB chilldown as compared to the experimental data (black lines) for the following temperature sensors: (a) TT202, (b) TT105, (c) TT162, (d) TT174, (e) TT156, (f) TT221.

4.2 Fault detection

One of the major applications of the algorithms in the cryogenic flow management is detection and evaluation of the system faults [37, 38]. Our approach to the fault management is based on the ability of the models to accurately predict fluid dynamics in the system, which was discussed above. To efficiently manage the faults, we developed [37, 38] approach based on the application of D-matrix, which is a causal representation of the relationship between faults and tests with 1 representing the relationship that the test can detect corresponding failure in the component and 0, otherwise.

Within this approach each fault is represented by a binary array of sensor readings, and each distinguishable failure mode has its own “signature” binary array. To account for strong variability of the fluid dynamics during chilldown, the original approach was extended [38] by introducing time-dependent D-matrix. The system failure modes were defined (as binary arrays) in several time windows corresponding to the different stages of chilldown.

The extended approach was verified by performing numerical analysis of 7 faults using 12 sensors readings [37, 38]. Here we provide an example of fault detection using this method. In this example the fault (valve BV3 is stuck open) is injected at t = 300 s. The system is simulated over 1600 s using homogeneous moving front model. At each time step, the model checks if the sensor readings stay within predefined maximum and minimum threshold limits corresponding to the nominal regime. Once the sensor reading crosses the thresholds, the fault is detected as shown in Figure 7. The failure mode is determined by comparing the binary array of the fault with the set of signature arrays.

Figure 7.

Fault detection in the transfer line when the valve BV3 is stuck opened. Top (red) and bottom (green) lines indicate margins of the nominal regime. Middle (blue) line correspond sensor readings during loading operation. The fault (BV3 stuck open) is injected at 300 s. Fault detection (crossing the margins) is shown by open circles.

4.3 Fault evaluation and model predictions

To successfully apply the fault detection and evaluation algorithm in real cryogenic systems, we need to validate the predictions of the models in the presence of faults. The validation was performed by conducting experiment with injected valve stuck open/closed faults in the CTB. In this experiment three chilldown regimes are explored, in which the opening of valve CV2 had values 15% (stuck closed), 25% (nominal), and 40% (stuck open). All other controls remain at the fixed values corresponding to the nominal flow considered in Section 4.1. During simulations the model parameters are kept constant at the values found during optimization described above. Only the value of the model parameter corresponding to the valve cv112 opening is modified according to experimental data.

The results of the comparison of the homogeneous model predictions of the unseen experimental data are shown in Figure 8. It can be seen from the figure that the model predicts reasonably well the main trends of the fault including the shift in the chilldown initiation and the change in the slopes of the time traces of temperature. Similar results were obtained using two-phase flow model [54].

Figure 8.

Comparison of the model predictions for the temperature time-series data (colored dashed lines) with the experimental results (colored solid lines) for the different openings of the valve cv112: (i) 15% blue line, (ii) 25% black line (nominal), and (iii) 40% red line.

We note, however, that the accuracy of the predictions remains limited for both models. To improve the accuracy, we proposed to embed the optimization algorithms discussed above within the machine learning framework. Such an approach would allow for semiautomated continuous learning of the model parameters using multiple databases and specifically databases obtained during fault modeling.


5. Conclusion

To summarize, we developed fast and reliable homogeneous moving front and separated two-fluid solvers for cryogenic loading operations. We proposed concise sets of cryogenic flow boiling correlations capable of reproducing a wide range of experimental time-series data obtained for two cryogenic transfer lines using both solvers.

Both solvers offer similar performance. Homogeneous solver is 10 times faster and can integrate 2 h of real-time loading in the CTB in less than 1 s. The two-phase flow algorithm provides access to important correlation parameters of two-phase flow, describes more accurately sensor readings, and can distinguish a number of non-equilibrium regimes unavailable in homogeneous flow including counterflow and vertical stratification.

We validated the performance of the models using horizontal transfer line at National Bureau of Standards [74] and large-scale experimental transfer line built at KSC [18]. We concluded that both solvers can be used for online control of cryogenic loading operations, and by applying two models in parallel, it is possible to resolve the trade-off issue between speed and accuracy of each solver.

We demonstrated an example of applications of these solvers to the analysis of cryogenic flow in the CTB transfer line including inference of the parameters of cryogenic correlations and chilldown modeling. We discussed an approach to the fault detection in cryogenic loading system-based D-matrix approach and demonstrated a capability of the solvers to resolve faults observed experimentally in the CTB system during chilldown.

We demonstrated that solvers allow for efficient optimization of the model parameters and discussed an extended approach to the inference of cryogenic correlations developed in our earlier work [44, 54]. We proposed that optimization algorithms should be included into machine learning framework for efficient offline learning of cryogenic correlations in the future work.

We emphasize that the machine learning approach will most likely underlie autonomous control and fault management of two-phase flows in the future space missions.



This work was supported by the Advanced Exploration Systems and Game Changing Development programs at NASA HQ.


A pipe cross-section
D pipe diameter
D h hydraulic pipe diameter
E specific energy
Gr Grashof number
H specific enthalpy
K n minor losses
N u Nusselt number
Pr Prandtl number
Re Reynolds number
T temperature
χ mass fraction
m ̇ mass flux
q ̇ heat flux
τ w friction losses
c specific heat
c p specific heat for constant pressure
g gravity
h heat transfer coefficient
h lg latent heat of evaporation
p pressure
u fluid velocity
Greek symbols
α void fraction
β liquid void fraction
Γ mass flow rate
k thermal conductivity
μ viscosity
ρ density
σ surface tension
τ shear stress
h heat transfer coefficient
e equilibrium
g gas
l liquid
s saturated
sub subcooled
w wall


  1. 1. Chato DJ. Cryogenic fluid transfer for exploration. Cryogenics. 2008;48(5–6):206-209
  2. 2. Collier JG, Thome JR. Convective Boiling and Condensation. Oxford, Clarendon Press; 1994
  3. 3. Ghiaasiaan SM. Two-phase flow, boiling, and condensation. In: Conventional and Miniature Systems. Cambridge: Cambridge University Press; 2007
  4. 4. Ishii M, Hibiki T. Thermo-Fluid Dynamics of Two-Phase Flow. Bücher: Springer; 2010
  5. 5. Two-Phase Flow Regimes. 2010. Available from:
  6. 6. TRACE V5.0 Theory Manual - Volume 1: Field Equations, Solution Methods, and Physical Models. U. S. Nuclear Regulatory Commission, Washington: Nuclear Regulatory Commission; 2008. ML120060218
  7. 7. RELAP5-3D Code Manual Volume I: Code Structure, System Models and Solution Methods, Revision 4.2, Idaho National Laboratory, Idaho Falls; 2014. INL-EXT-98-00834-V1
  8. 8. Nourgaliev R. Solution Algorithms for Averaged Equations. Idaho National Laboratory, Idaho Falls; 2012. INL/EXT-12-27187
  9. 9. Cheng L, Ribatski G, Wojtan L, Thome JR. New flow boiling heat transfer model and flow pattern map for carbon dioxide evaporating inside horizontal tubes. International Journal of Heat and Mass Transfer. 2006;49(21–22):4082-4094
  10. 10. Prosperetti A, Tryggvason G. Computational Methods for Multiphase Flow. Cambridge: Cambridge University Press; 2007
  11. 11. Van Dresar NT, Siegwarth JD, Hasan MM. Convective heat transfer coefficients for near-horizontal two-phase flow of nitrogen and hydrogen at low mass and heat flux. Cryogenics. 2001;41(11–12):805-811
  12. 12. Jackson JK. Cryogenic Two-Phase Flow During Chilldown: Flow Transition and Nucleate Boiling Heat Transfer. Gainesville: University of Florida; 2006
  13. 13. Wang S, Wen J, Li Y, Yang H, Li Y, Tu J. Numerical prediction for subcooled boiling flow of liquid nitrogen in a vertical tube with MUSIG model. Chinese Journal of Chemical Engineering. 2013;21(11):1195-1205
  14. 14. Darr SR, Hu H, Shaeffer R, Chung J, Hartwig JW, Majumdar AK. Numerical Simulation of the Liquid Nitrogen Chilldown of a Vertical Tube. AIAA SciTech. Kissimmee, FL: American Institute of Aeronautics and Astronautics; 2015
  15. 15. Konishi C, Mudawar I. Review of flow boiling and critical heat flux in microgravity. International Journal of Heat and Mass Transfer. 2015;80:469-493
  16. 16. Notardonato W. Active control of cryogenic propellants in space. Cryogenics. 2012;52(4–6):236-242
  17. 17. Kim J. Review of reduced gravity boiling heat transfer: US research. Journal of The Japan Society of Microgravity Application. 2003;20(4):264-271
  18. 18. Johnson R, Notardonato W, Currin K, Orozco-Smith E. Integrated ground operations demonstration units testing plans and status. In: AIAA SPACE 2012 Conference & Exposition. SPACE Conferences & Exposition. American Institute of Aeronautics and Astronautics; 2012
  19. 19. Hartwig JW, Vera J. Numerical modeling of the transient chilldown process of a cryogenic propellant transfer line. In: 53rd AIAA Aerospace Sciences Meeting. AIAA SciTech. American Institute of Aeronautics and Astronautics; 2015
  20. 20. Kawaji M. Boiling heat transfer during quenching under microgravity. In: Fluid Dynamics Conference. Fluid Dynamics and Co-located Conferences. American Institute of Aeronautics and Astronautics; 1996
  21. 21. Majumdar A, Steadman T. Numerical modeling of thermofluid transients during chilldown of cryogenic transfer lines. In: 33rd International Conference on Environmental Systems, Vol. 3; Vancouver, B.C.: Society of Automotive Engineers; 2003
  22. 22. Majumdar AK, Ravindran SS. Fast, nonlinear network flow solvers for fluid and thermal transient analysis. International Journal of Numerical Methods for Heat and Fluid Flow. 2010;20(6–7):617-637
  23. 23. Chung JN, Yuan K. Recent progress on experimental research of cryogenic transport line chilldown process. Frontiers in Heat and Mass Transfer. 2015;6:1-17
  24. 24. Darr S, Dong J, Glikin N, Hartwig J, Majumdar A, Leclair A, et al. The effect of reduced gravity on cryogenic nitrogen boiling and pipe chilldown. NPJ Microgravity. 2016;2:16033
  25. 25. Hedayatpour A, Antar BN, Kawaji M. Cool-down of a vertical line with liquid-nitrogen. Journal of Thermophysics and Heat Transfer. 1993;7(3):426-434
  26. 26. Yuan K, Ji Y, Chung JN. Numerical modeling of cryogenic chilldown process in terrestrial gravity and microgravity. International Journal of Heat and Fluid Flow. 2009;30(1):44-53
  27. 27. Cullimore BA. Optimization and automated data correlation in the nasa standard thermal/fluid system analyzer. In: Proceedings 33rd Intersociety Engineering Conference on Energy Conversion; 1998
  28. 28. Majumdar A, Ravindran SS. Numerical prediction of conjugate heat transfer in fluid network. Journal of Propulsion and Power. 2011;27(3):620-630
  29. 29. SINDA/FLUINT, General Purpose Thermal/Fluid Network Analyzer, Version 5.8. Boulder, Colorado: C&R Technologies, Inc; 2015
  30. 30. Hafiychuk V, Foygel M, Ponizovskaya-Devine E, Smelyanskiy V, Watson MD, Brown B, et al. Moving-boundary model of cryogenic fuel loading, I: Two-phase flow in a pipe. Journal of Thermophysics and Heat Transfer. 2015;29:533-544
  31. 31. Hafiychuk V, Foygel M, Ponizovskaya-Devine E, Smelyanskiy V, Watson MD, Brown B, et al. Moving-boundary model of cryogenic fuel loading, II: Theory versus experiments. Journal of Thermophysics and Heat Transfer. 2015;29:545-550
  32. 32. Grald EW, MacArthur JW. A moving-boundary formulation for modeling time-dependent two-phase flow. International Journal of Heat and Fluid Flow. 1992;13:266-272
  33. 33. Doherty MP, Gaby JD, Salerno LJ, Sutherlin SG. Cryogenic Fluid Management Technology for Moon and Mars Missions. AIAA SPACE 2009 Conference & Exposition, Pasadena, California; 2009
  34. 34. Schmidt FW, Henderson RE, Wolgemuth CH. Introduction to Thermal Sciences; New York: Wiley, John \& Sons, Incorporated; 1993
  35. 35. Thome JR, The heat transfer engineering data book III, Wieland-Werke AG, Edition 3, Ulm; 2016
  36. 36. Osipov VV, Hafiychuk H, Ponizovskaya-Devine E, Khasin M, Smelyanskiy V. Risk Assessment and Scaling for the SLS LOx ET. NASA, ARC; 2012. NASA/TP—2012–216485
  37. 37. Ponizovskaya-Devine E, Hafiychuk VV, Luchinsky DG, Khasin M, Perotti J, Brown B. Fault diagnostics and evaluation in cryogenic loading system using optimization algorithm. In: Annual Conference of the Prognostics and Health Management Society; 2015
  38. 38. Kodali A, Ponizovskaya-Devine E, Robinson P, Luchinsky D, Bajwa A, Khasin M, et al. D-matrix based fault diagnostics for cryogenic loading systems. In: Annual Conference of the Prognostics and Health Management Society; 2015
  39. 39. Luchinsky DG, Smelyanskiy VN, Brown B. Physics based model for cryogenic chilldown and loading. Part IV: Code structure. 2014. NASA, ARC. NASA/TP-2014-218399
  40. 40. Luchinsky DG, Smelyanskiy VN, Brown B. Physics based model for cryogenic chilldown and loading. Part II: Verification and validation. 2014. NASA, ARC. NASA/TP-2014
  41. 41. Ponizovskaya-Devine E, Luchinsky DG, Khasin M, Perotti J, Sass J, Brown B. Towards physics based autonomous control of the cryogenic propellant loading system. In: 51st AIAA/SAE/ASEE Joint Propulsion Conference. Propulsion and Energy Forum. American Institute of Aeronautics and Astronautics; 2015
  42. 42. Luchinskiy DG, Ponizovskaya-Devine E, Hafiychuk V, Kashani A, Khasin M, Timucin D, et al. Hierarchy of two-phase flow models for autonomous control of cryogenic loading operation. IOP Conference Series: Materials Science and Engineering. 2015;101(1):012069
  43. 43. Luchinsky D, Ponizovskaya-Devine E, Khasin M, Kodali A, Perotti J, Sass J, et al. Two-phase flow modelling of the cryogenic propellant loading system. In: 51st AIAA/SAE/ASEE Joint Propulsion Conference. Propulsion and Energy Forum. American Institute of Aeronautics and Astronautics; 2015. p. 4214
  44. 44. Luchinsky DG, Khasin M, Timucin D, Sass J, Brown B. Inferential framework for two-fluid model of cryogenic chilldown. International Journal of Heat and Mass Transfer. 2017;114:796-808
  45. 45. Luchinsky DG, Khasin M, Timucin D, Sass J, Johnson RG, Perotti J, et al. Physics based model for cryogenic chilldown and loading. Part III: Correlations. 2016. NASA, ARC. NASA/TP-2016
  46. 46. Dinh TN, Nourgaliev RR, Theofanous TG. Understanding of the Ill-posed two-fluid model. In: Proceedings of the Tenth International Topical Meeting on Nuclear Reactor Thermal Hydraulics. KNS; 2003. p. 1CD-ROM
  47. 47. Liou MS, Nguyen L, Chang CH, Sushchikh S, Nourgaliev R, Theofanous T. Hyperbolicity, discontinuities, and numerics of two-fluid models. In: Deconinck H, Dick E, editors. Springer, Berlin: Computational Fluid Dynamics 2006. 2009. p. 625
  48. 48. Cordier F, Degond P, Kumbaro A. Phase appearance or disappearance in two-phase flows. Journal of Scientific Computing. 2014;58(1):115-148
  49. 49. Harlow FH, Amsden AA. A numerical fluid dynamics calculation method for all flow speeds. Journal of Computational Physics. 1971;8(2):197-213
  50. 50. Liles DR, Reed WH. A semi-implicit method for two-phase fluid dynamics. Journal of Computational Physics. 1978;26(3):390-407
  51. 51. Liou MS, Chang CH, Nguyen L, Theofanous TG. How to solve compressible multifluid equations: A simple, robust, and accurate method. AIAA Journal. 2008;46(9):2345-2356
  52. 52. Wojtan L, Ursenbacher T, Thome JR. Investigation of flow boiling in horizontal tubes: Part I—A new diabetic two-phase flow pattern map. International Journal of Heat and Mass Transfer. 2005;48(14):2955-2969
  53. 53. Nellis G, Klein S. Heat Transfer. Cambridge: Cambridge University Press; 2009
  54. 54. Luchinsky DG, Khasin M, Timucin D, Sass J, Johnson RG, Perotti J, et al. Physics based model for cryogenic chilldown and loading. Part V: Optimization techniques. 2016. NASA, ARC. NASA/TP-2016
  55. 55. Franchello G. Development of a Heat Transfer Package Applicable to a Large Variety of Fluids. Luxembourg: European Commission; 1993. EUR 14985 EN
  56. 56. Seader JD, Miller WS, Kalvinskas LA. Boiling Heat Transfer for Cryogenics. NASA, Washington: NASA CR-243; 1965
  57. 57. Griffith P. The Correlation of Nucleate Boiling Burnout Data. Massachusetts Institute of Technology; 1957. NS-035-267
  58. 58. Zuber N, Tribus M. Further Remarks on the Stability of Boiling. Oak Ridge, AECU-3631: Technical Information Service Extension; 1958. 68 pages
  59. 59. Kutateladze SS. Critical heat flux to flowing, wetting, subcooled liquids. Energetika. 1959;2:229-239
  60. 60. RELAP5-3D Code Manual Volume IV: Models and Correlations, Revision 4.2, Idaho National Laboratory, Idaho Falls; 2014. INL-EXT-98-00834-V4
  61. 61. Holman JP. Heat Transfer. Mechanical Engineering Series. Boston: McGraw-Hill; 1989
  62. 62. Sato T, Matsumura H. On the conditions of incipient subcooled-boiling with forced convection. Bulletin of JSME. 1964;7(26):392-398
  63. 63. Frost W, Dzakowic GS. An Extension of the Method for Predicting Incipient Boiling on Commercially Finished Surfaces. New York: ASME; 1967. p. 67-HT-61
  64. 64. Huang L. Evaluation of onset of nucleate boiling models. In: ECI International Conference on Boiling Heat Transfer, Vol. 40. INUS; 2009. p. 40079216
  65. 65. Tong LS, Tang YS. Boiling Heat Transfer and Two-Phase Flow. Series in Chemical and Mechanical Engineering. New York: Taylor & Francis; 1997
  66. 66. Iloeje OC, Plummer DN, Rohsenow WM, Griffith P. Effects of mass flux, flow quality, thermal and surface properties of materials on rewet of dispersed flow film boiling. Journal of Heat Transfer. 1982;104(2):304-308. DOI: 10.1115/1.3245088
  67. 67. Theler G, Freis D. Theoretical critical heat flux prediction based on non-equilibrium thermodynamics considerations of the subcooled boiling phenomenon. Mecánica Computacional. 2011;XXX:1713-1732. Möller O, Signorelli JW, Storti MA, editors
  68. 68. Carbajo JJ. A study on the rewetting temperature. Nuclear Engineering and Design. 1985;84(1):21-52
  69. 69. Berenson PJ. Film-boiling heat transfer from a horizontal surface. Journal of Heat Transfer. 1961;83:351-358
  70. 70. Churchill SW. Friction-factor equation spans all fluid-flow regimes. Chemical Engineering. 1977;84(24):91-92
  71. 71. Chisholm D. A theoretical basis for the Lockhart-Martinelli correlation for two-phase flow. International Journal of Heat and Mass Transfer. 1967;10(12):1767-1778
  72. 72. Kashani A, Ponizhovskaya E, Luchinsky D, Smelyanskiy V, Sass J, Brown B, et al. Physics based model for online fault detection in autonomous cryogenic loading system. AIP Conference Proceedings. 2014;1573(1):1305-1310
  73. 73. Kashani A, Luchinskiy DG, Ponizovskaya-Devine E, Khasin M, Timucin D, Sass J, et al. Optimization of cryogenic chilldown and loading operation using SINDA/FLUINT. IOP Conference Series: Materials Science and Engineering. 2015;101(1):012115
  74. 74. Brennan JA, Brentari EG, Smith RV, Steward WG. Cooldown of Cryogenic Transfer Lines. Boulder, Colorado: National Bureau of Standards; 1966. NBS Report 9264

Written By

Ekaterina Ponizovskaya-Devine, Dmitry Luchinsky, Michael Foygel, Vasil Hafiychuk, Michae Khasin, Jared Sass and Barbara Brown

Submitted: 30 October 2018 Reviewed: 23 January 2019 Published: 15 March 2019