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
- moving front
- two-phase flow
- flow boiling
- heat transfer
Deep space exploration requires development of autonomous management system  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 in the system is common. The temperature of the wall , vapor , and liquid 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 .
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  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 .
Despite these efforts small-time steps and instabilities  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  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 .
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  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 in one dimension:
Here is enthalpy, is density, and is velocity of the flow. is the friction losses, is the pipe length, is the gravity, is the angle of the pipe with horizontal plane, is the heat transfer from the pipe wall to the fluid mixture, and is the cross-section area of the pipe.
The wall temperature is described by the energy conservation equation:
where is the cross-sectional area of the pipe wall with diameter , which is in thermal contact with the ambient fluid of temperature and with the fluid of temperature flowing in the pipe. is the specific heat, is density, and is thermal conductivity of the wall material. is the heat transfer coefficient at the inner wall surface, and 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 for the liquid. Each volume can contain pure liquid with , pure gas with , 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 ) in the two-phase region are
where is the slip ratio. It can be shown that for the steady flow without slip (), the local density and specific enthalpy of the mixture are :
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 in the form of tables and by providing a set of cryogenic correlations discussed below.
where the factor is
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 as , where , is Nusselt number, and is Prandtl number.
For the fully developed laminar flow, . For the turbulent flow , 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 
for . Here the effective Reynolds number for each of the components, , is estimated by means of the average hydraulic diameters and .
For the forced convection, the heat transfer coefficient  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 in terms of pressure and enthalpy , we can rewrite the mass and energy conservation equations in terms of pressure and enthalpy or pressure and void fraction . In the case of two-phase control volume, the equations are (with ):
and in the case of single-phase control volume (e.g., vapor), we have
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 results in the following set of equations:
Here and 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 . The wall temperature is calculated for each region according to Eq. (2).
The mass fraction 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 , and the new density is defined by the void fraction in Eq. (4), while the relation between 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 .
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 , position of the front , void fraction , vapor () and liquid () enthalpies , and two wall temperatures before and after the front and . These variables are calculated at the center of control volume.
Next the mass and enthalpy fluxes are calculated at the boundaries of each control volume; see Figure 2, using the following equation:
The subindex for frictional factor and minor losses 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 .
For each control volume , the continuity condition is added as follows: for the serial connections of the , , and if the -th control volume consists of two parallel control volumes and .
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 , , , , and 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 () and non-equilibrium () 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 , 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].
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 . 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.
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 . The non-hyperbolicity was suppressed using so-called virtual mass term . 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.  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 .
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 . 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 of the mass transfer at the wall and at the interface:
where is the difference between saturated gas and liquid enthalpies for positive flux while it is for negative flux.
The heat transfer coefficients for the gas (subscript ) and liquid () at the wall () and at the interface () are defined by the following relations:
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.
are Prandtl and Grashof numbers, respectively, is the coefficient of thermal expansion, and is the hydraulic diameter. The maximum value of 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 .
where , , is the wall superheat, and is liquid subcooling temperature.
When the wall superheat exceeds , 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 .
For the pool boiling value of , we used Zuber  correlation
The temperature for the critical heat flux was estimated in this work using the approach proposed by Theler and Freis .
where is the isoentropic expansion factor that for ideal diatomic gases is 7/2 and is the specific gas constant. was further corrected for boiling flows using Iloeje-type correlations  similar to the one applied , see below.
When wall superheat exceeds , 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 .
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) , (ii) , (iii) , (iv) , and (v) .
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 discussed above with appropriately modified parameters. Transition to the single-phase heat transfer is initiated when dry-out transition is detected.
In the regime where the wall superheat is increasing from to , the heat flux was defined using the following interpolation :
Similarly, for the regime where the wall superheat is increasing from to , the transition boiling heat flux was interpolated as .
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 .
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 and with the friction factors for turbulent and laminar flow given by Churchill approximation .
where is given by
Friction factors are the same as in above, while coefficients and depend on the flow pattern .
The interface drag is given by , where interfacial drag coefficient depends on the flow pattern .
The algorithms discussed above were used to develop a number of applications related to the management of cryogenic flows including, e.g., fault detection , diagnostics , evaluation , 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)  and (ii) large-scale experimental transfer line build at KSC .
The first system  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 .
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 .
To address this problem, we developed inferential framework  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 . 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.
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  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  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 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.
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 .
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.
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  and large-scale experimental transfer line built at KSC . 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.
|D h||hydraulic pipe diameter|
|K n||minor losses|
|N u||Nusselt number|
|m ̇||mass flux|
|q ̇||heat flux|
|τ w||friction losses|
|c p||specific heat for constant pressure|
|h||heat transfer coefficient|
|h lg||latent heat of evaporation|
|β||liquid void fraction|
|Γ||mass flow rate|
|h||heat transfer coefficient|