## Abstract

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.

### Keywords

- 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

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

Here

The wall temperature is described by the energy conservation equation:

where

Two-phase flow is considered to be homogeneous with void fraction of the vapor equal

where

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:

The model equations are closed by adding equation of state

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]

where the factor

and

The estimated average void fraction

For the fully developed laminar flow,

for

for

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

where laminar and turbulent contributions to the Nusselt number

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

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

Here subindexes

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

Here

The mass fraction

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

Next the mass

The subindex

For each control volume

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

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 (

## 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]

and liquid phases

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

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

where

The heat transfer coefficients

For relatively low mass fluxes considered in this work (

**Convective heat transfer** coefficient

are Prandtl and Grashof numbers, respectively,

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

**The onset of nucleate boiling** correlations is based on the analysis of the balance between mechanical and thermodynamical equilibrium [3]. The corresponding

where

When the wall superheat exceeds

**The critical heat flux**

For the pool boiling value of

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].

where

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

The temperature

where

When wall superheat exceeds

**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

Iloeje [55, 66] has corrected Berenson equation to take into account the dependence of the

where

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

Typical values of the parameters used in simulations are the following: (i)

The minimum film boiling heat flux can now be defined as

In the region of single-phase gas flow, the heat transfer coefficient is given by

**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

where

Similarly, for the regime where the wall superheat is increasing from

where

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

The two-phase pressure drop

where

Friction factors

The interface drag is given by

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

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.

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

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.

### 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

### 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].

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

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.

## Acknowledgments

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 |

Subscripts | |

e | equilibrium |

g | gas |

l | liquid |

s | saturated |

sub | subcooled |

w | wall |