Two-Phase Flow Modeling of Cryogenic Loading Operations

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.

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. $\alpha$ 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.

$$\begin{array}{l}{\rho }_{,t}+{\left(\mathit{\rho u}\right)}_{,x}=0\\ {\left(\mathit{\rho u}\right)}_{,t}+{\left(u\left(\mathit{\rho u}\right)\right)}_{,x}+p_x+\frac{1}{A}{\tau }_w{l}_w+\mathit{\rho g}\text{sin}\theta =0\\ {\left(\mathit{\rho e}\right)}_{,t}+{\left(\mathit{\rho uH}\right)}_{,x}+\frac{1}{A}{\dot{q}}_w{l}_w=0\end{array}$$ E1

Here $H$ is enthalpy, $\rho$ is density, and $u$ is velocity of the flow. ${\tau }_w$ is the friction losses, $l_w$ is the pipe length, $g$ is the gravity, $\theta$ is the angle of the pipe with horizontal plane, ${\dot{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{\rho }_w{A}_w{\left(T_w\right)}_{,t}=\mathit{h\pi D}\left(T-T_w\right)+h_{\mathit{amb}}\mathit{\pi D}\left(T_{\mathit{amd}}-T_w\right)+\kappa A_w{\left(T\right)}_{,\mathit{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_{\mathit{amb}}$ and with the fluid of temperature $T$ flowing in the pipe. $c_w$ is the specific heat, ${\rho }_w$ is density, and $\kappa$ is thermal conductivity of the wall material. $h$ is the heat transfer coefficient at the inner wall surface, and $h_{\mathit{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 $\alpha$ and $\left(1-\alpha \right)$ for the liquid. Each volume can contain pure liquid with $\alpha =0$ , pure gas with $\alpha =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

$$\begin{array}{l}\alpha =\frac{A_v}{A}=\frac{A_v}{A_v+A_l}\text{and}x=\frac{M_v}{M_v+M_l}={\left(1+\frac{\left(1-\alpha \right){\rho }_v}{\alpha /{\rho }_{l}S}\right)}^{-1}\end{array}$$ 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]:

$$\begin{array}{l}\rho ={\mathit{\alpha \rho }}_v+\left(1-\alpha \right){\rho }_l;H={\mathit{xH}}_v+\left(1-x\right)H_l\end{array}$$ 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:

$$\begin{array}{l}{\left(\mathit{A\rho }{u}^2\right)}_{,x}+\frac{1}{A}{\left({\tau }_w{l}_w\right)}_{\mathit{tp}}=p_{,x}-\mathit{\rho gsin}\left(\theta \right)\end{array}$$ E5

The model equations are closed by adding equation of state ${\rho }_{g,l}$ $=$ ${\rho }_{g,l}\left(p,H_{v,l}\right)$ 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]

$$\begin{array}{l}\alpha \left(x,p\right)={\left(1+a\left(p\right)\left(\frac{1}{x}-1\right)\right)}^{-1}\end{array}$$ E6

where the factor $a$ is

$$\begin{array}{l}a={\left(\frac{{\rho }_v}{{\rho }_l}\right)}^{n1}{\left(\frac{{\mu }_l}{{\mu }_v}\right)}^{n2}\end{array}$$ E7

and $\mu$ is the viscosity of the fluid. Coefficients $n_1$ and $n_2$ can have different values, e.g., according to Thome [35] $n1=0.89,n2=0.18$ , while it is suggested to use $n1=0.67,n2=0$ in [34].

The estimated average void fraction $\alpha$ can be used to evaluate the heat transfer coefficient from the tube wall to the TP mixture in the TP region $\left(n=2\right)$ as $h_2=\alpha h_v+\left(1-\alpha \right)h_l$ , where $h_{v,l}={\left(\frac{\kappa }{D_h}\mathit{Nu}\right)}_{v,l}$ , $\mathit{Nu}$ is Nusselt number, and $\mathit{Pr}$ is Prandtl number.

For the fully developed laminar flow, $\mathit{Nu}=3.66$ . For the turbulent flow $\left(\mathit{Re}>2300\right)$ , 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]

$$\mathit{Nu}=0.0214\left({\mathit{Re}}^{4/5}-100\right){\mathit{Pr}}^{2/5}\left(1+{\left(D_h/L_2\right)}^{2/3}\right)$$ E8

for $0.5<\mathit{Pr}<1.5$ and

$$\mathit{Nu}=0.012\left({\mathit{Re}}^{0.87}-280\right){\mathit{Pr}}^{2/5}\left(1+{\left(D_h/L_2\right)}^{2/3}\right)$$ E9

for $1.5<\mathit{Pr}<500$ . Here the effective Reynolds number for each of the components, ${\mathit{Re}}_{v,l}=\left(\frac{\dot{m}{D}_h}{\mathit{\mu A}}\right)$ , is estimated by means of the average hydraulic diameters $D_{h,v}=\mathit{\alpha D}$ and $D_{h,l}=\left(1-\alpha \right)D$ .

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

$$h_{\mathit{amb}}=\frac{2\kappa }{\mathit{\pi D}}{\mathit{Nu}}_{\mathit{amb}},{\mathit{Nu}}_{\mathit{amb}}={\mathit{Nu}}_{\mathit{lam}}+{\mathit{Nu}}_{\mathit{\text{turb}}},$$ E10

where laminar and turbulent contributions to the Nusselt number

$${\mathit{Nu}}_{\mathit{lam}}=0.0664{\mathit{Re}}^{1/2}{\mathit{Pr}}^{1/3}\mathit{\text{and}}{\mathit{Nu}}_{\mathit{\text{turb}}}=\frac{0.037{\mathit{Re}}^{0.8}\mathit{Pr}}{1+2.443{\mathit{Re}}^{-0.1}\left({\mathit{Pr}}^{2/3}-1\right)}$$ 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 $\rho$ 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 $\alpha$ . In the case of two-phase control volume, the equations are (with $\beta =1-\alpha$ ):

$$\begin{array}{l}\left(\alpha \frac{\partial {\rho }_v}{\partial p}+\beta \frac{\partial {\rho }_l}{\partial p}\right)\frac{\mathit{dp}}{\mathit{dt}}+\left({\rho }_v-{\rho }_l\right)\frac{\mathit{d\alpha }}{\mathit{dt}}=\frac{{\dot{m}}_L-{\dot{m}}_R}{\mathit{LA}}\\ \left(\alpha \frac{\partial {\left(\mathit{\rho H}\right)}_v}{\partial p}+\beta \frac{\partial {\left(\mathit{\rho H}\right)}_l}{\partial p}-1\right)\frac{\mathit{dp}}{\mathit{dt}}+\left({\rho }_v{H}_v-{\rho }_l{H}_l\right)\frac{\mathit{d\alpha }}{\mathit{dt}}=\frac{{\dot{m}}_L{H}_L-{\dot{m}}_R{H}_R}{\mathit{LA}}+\frac{4h\left(T_w-T\right)}{D}\end{array}$$ E12

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

$$\begin{array}{l}\frac{\partial \rho }{\partial p}\frac{\mathit{dp}}{\mathit{dt}}+\frac{\partial \rho }{\partial H}\frac{\mathit{dH}}{\mathit{dt}}=\frac{{\dot{m}}_L-{\dot{m}}_R}{\mathit{LA}}\\ \left(\frac{\partial \left(\mathit{\rho H}\right)}{\partial p}-1\right)\frac{\mathit{dp}}{\mathit{dt}}+\left(\frac{\partial \mathit{\rho H}}{\partial H}\right)\frac{\mathit{dH}}{\mathit{dt}}=\frac{{\dot{m}}_L{H}_L-{\dot{m}}_R{H}_R}{\mathit{LA}}+\frac{4h\left(T_w-T\right)}{D}\end{array}$$ E13

Here subindexes $R$ and $L$ refer to the right and left boundaries; ${\rho }_v$ and ${\rho }_l$ are the saturated density for vapor and liquid in Eq. (12), while in Eq. (13) for the single-phase, $\rho$ 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:

$$\begin{array}{l}\left(\alpha \frac{\partial {\rho }_v}{\partial P}+\beta \frac{\partial {\rho }_l}{\partial p}\right)\frac{\mathit{dp}}{\mathit{dt}}+\left({\rho }_v-{\rho }_l\right)\left(\frac{\mathit{d\alpha }}{\mathit{dt}}-\beta \frac{\mathit{dln}\left(L_2\right)}{\mathit{dt}}\right)=\frac{{\dot{m}}_L-{\dot{m}}_R}{L_{2}A}\\ \left(\alpha \frac{\partial {\rho }_v{H}_v}{\partial p}+\beta \frac{\partial {\rho }_l{H}_l}{\partial p}-1\right)\frac{\mathit{dp}}{\mathit{dt}}+\left({\rho }_v{H}_v-{\rho }_l{H}_l\right)\left(\frac{\mathit{d\alpha }}{\mathit{dt}}-\beta \frac{\mathit{dln}\left(L_2\right)}{\mathit{dt}}\right)\\ =\frac{{\dot{m}}_L{H}_L-{\dot{m}}_R{H}_v}{L_{2}A}+\frac{4h}{D}\left(T_w-T\right)\\ \frac{\partial \rho }{\partial p}\frac{\mathit{dp}}{\mathit{dt}}+\frac{\partial \rho }{\partial H}\frac{\mathit{dH}}{\mathit{dt}}+\left(\rho -{\rho }_v\right)\frac{\mathit{dln}\left(L_3\right)}{\mathit{dt}})=\frac{{\dot{m}}_L-{\dot{m}}_R}{\left(L_3\right)A}\\ \left(\frac{\partial \mathit{\rho H}}{\partial p}-1\right)\frac{\mathit{dp}}{\mathit{dt}}+\frac{\partial \mathit{\rho H}}{\partial H}\frac{\mathit{dH}}{\mathit{dt}}+\left(\mathit{\rho H}-{\rho }_v{H}_v\right)\frac{\mathit{dln}\left(L_3\right)}{\mathit{dt}}=\frac{{\dot{m}}_L{H}_l-{\dot{m}}_R{H}_R}{\left(L_3\right)A}+\frac{4h}{D}\left(T_w-T\right)\\ \end{array}$$ 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 $\alpha$ in Eq. (4), while the relation between $x$ and $\alpha$ 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 $\alpha$ , vapor ( $v$ ) and liquid ( $l$ ) enthalpies $H_{v\left(l\right)}$ , and two wall temperatures before and after the front $T_{w1}$ and $T_{w2}$ . These variables are calculated at the center of control volume.

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

$$\begin{array}{l}\frac{{\dot{m}}_j^2}{A_n^2}\left(\frac{1}{{\rho }_n}\left(f_n\frac{L_n}{D_n}+\sum K_n\right)+{\left(\frac{x^2}{{\mathit{\gamma \rho }}_g}+\frac{{\left(1-x\right)}^2}{\left(1-\gamma \right){\rho }_l}\right)}_n\right)=-\mathrm{\Delta }{p}_j-\mathrm{\Delta }{z}_j{\rho }_j\mathit{gsin\theta }\end{array}$$ 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: ${\dot{m}}_{k,\mathit{\text{in}}}={\dot{m}}_{k-1,\mathit{\text{out}}}$ for the serial connections of the $k-1$ , $k$ , and ${\dot{m}}_{k-1,\mathit{\text{out}}}={\dot{m}}_{k1,\mathit{\text{in}}}+{\dot{m}}_{k2,\mathit{\text{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$ , $\alpha$ , 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\ne u_l$ ) and non-equilibrium ( $T_g\ne 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]

$$\begin{array}{l}{\left(A{\mathit{\alpha \rho }}_g\right)}_{,t}+{\left(A{\mathit{\alpha \rho }}_g{u}_g\right)}_{,x}=A{\mathrm{\Gamma }}_g\\ {\left(A{\mathit{\alpha \rho }}_g{u}_g\right)}_{,t}+{\left(A{\mathit{\alpha \rho }}_g{u}_g^2\right)}_{,x}+\mathit{A\alpha }{p}_{,x}=A{\mathrm{\Gamma }}_g{u}_{\mathit{ig}}-\mathit{A\alpha g}{\rho }_g{z}_{,x}-{\tau }_{\mathit{gw}}{l}_{\mathit{wg}}-{\tau }_{\mathit{gi}}{l}_i\\ {\left(A{\mathit{\alpha \rho }}_g{E}_g\right)}_{,t}+{\left(A{\mathit{\alpha \rho }}_g{E}_g{u}_g\right)}_{,x}=A{\mathrm{\Gamma }}_g{H}_g-\mathit{Ap}{\alpha }_{,t}-{\left(\mathit{pA\alpha }{u}_g\right)}_{,x}+{\dot{q}}_{\mathit{gw}}{l}_{\mathit{wg}}+{\dot{q}}_{\mathit{gi}}{l}_i\end{array}$$ E16

and liquid phases

$$\begin{array}{l}{\left(A{\mathit{\beta \rho }}_l\right)}_{,t}+{\left(A{\mathit{\beta \rho }}_l{u}_l\right)}_{,x}=-A{\mathrm{\Gamma }}_g\\ {\left(A{\mathit{\beta \rho }}_l{u}_l\right)}_{,t}+{\left(A{\mathit{\beta \rho }}_l{u}_l^2\right)}_{,x}+\mathit{A\beta }{p}_{,x}=-\mathit{A\beta g}{\rho }_l{z}_{,x}-{\tau }_{\mathit{lw}}{l}_{\mathit{wl}}-{\tau }_{\mathit{li}}{l}_i-A{\mathrm{\Gamma }}_g{u}_{\mathit{il}}\\ {\left(\mathit{A\beta }{E}_l{\rho }_l\right)}_{,t}+{\left(\mathit{A\beta }{E}_l{\rho }_l{u}_l\right)}_{,x}=-\mathit{Ap}{\beta }_{,t}-{\left(\mathit{pA\beta }{u}_l\right)}_{,x}+{\dot{q}}_{\mathit{lw}}{l}_{\mathit{wl}}+{\dot{q}}_{\mathit{li}}{l}_i-A{\mathrm{\Gamma }}_g{H}_l.\\ \end{array}$$ 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$ :

$$\begin{array}{l}{\rho }_w{c}_w{d}_w\frac{\partial T_w}{\partial t}=h_{\mathit{wg}}\left(T_g-T_w\right)+h_{\mathit{wl}}\left(T_l-T_w\right)+h_{\mathit{amb}}\left(T_{\mathit{amb}}-T_w\right).\end{array}$$ 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 ${\mathrm{\Gamma }}_g={\mathrm{\Gamma }}_{\mathit{wg}}+{\mathrm{\Gamma }}_{\mathit{ig}}$ of the mass transfer at the wall and at the interface:

$${\mathrm{\Gamma }}_{\mathit{wg}}=\frac{{\dot{q}}_{\mathit{wl}}}{H_g^{\ast }-H_l^{\ast }};{\mathrm{\Gamma }}_{\mathit{ig}}=\frac{{\dot{q}}_{\mathit{li}}+{\dot{q}}_{\mathit{gi}}}{H_g^{\ast }-H_l^{\ast }};$$ E19

where $H_g^{\ast }-H_l^{\ast }$ is the difference between saturated gas and liquid enthalpies $H_{g,s}-H_l$ for positive flux $\mathrm{\Gamma }>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:

$$\begin{array}{l}{\dot{q}}_{\mathit{wg}}=h_{\mathit{wg}}\left(T_w-T_g\right);{\dot{q}}_{\mathit{ig}}=h_{\mathit{ig}}\left(T_{l,s}-T_g\right);\\ {\dot{q}}_{\mathit{wl}}=h_{\mathit{wl}}\left(T_w-T_l\right);{\dot{q}}_{\mathit{il}}=h_{\mathit{il}}\left(T_{l,s}-T_l\right).\end{array}$$ 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_{\mathit{cb}}$ was calculated for four flow regimes [6, 60]. The value of $\left(D_h{h}_{\mathit{cb}}/\kappa \right)$ was taken to be $\sim$ 4.36 in the forced laminar regime, while in forced turbulent regime, it was calculated as $0.023\cdot {\mathit{Re}}^{0.8}{\mathit{Pr}}^{0.4}$ [53]. The coefficient $\left(D_h{h}_{\mathit{cb}}/\kappa \right)$ for natural convection was calculated as $0.1{\left(\mathit{Gr}\cdot \mathit{Pr}\right)}^{1/3}$ for laminar and as $0.59{\left(\mathit{Gr}\cdot \mathit{Pr}\right)}^{1/4}$ for turbulent flows [61]. Here

$$\mathit{Pr}=\frac{\mu C_p}{\kappa }\text{and}\mathit{Gr}=\frac{{\rho }^{2}g{\beta }_T\left(T_w-T_{l\left(g\right)}\right)D^3}{{\mu }^2}$$ E21

are Prandtl and Grashof numbers, respectively, ${\beta }_T$ is the coefficient of thermal expansion, and $D_h$ is the hydraulic diameter. The maximum value of $h_{\mathit{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_{\mathit{onb}}$ .

The onset of nucleate boiling correlations is based on the analysis of the balance between mechanical and thermodynamical equilibrium [3]. The corresponding $T_{\mathit{onb}}$ and heat flux ${\dot{q}}_{\mathit{onb}}$ are [62, 63, 64].

$$T_{\mathit{onb}}=T_s+F\left(1+\sqrt{1+\frac{2\mathrm{\Delta }{T}_{\mathit{sub}}}{F}}\right),{\dot{q}}_{\mathit{onb}}=\frac{B}{{\mathit{Pr}}^2}\mathrm{\Delta }{T}_{\mathit{sat}}^2=h_{\mathit{cb}}\left(T_{\mathit{onb}}-T_l\right)$$ E22

where $B=\frac{{\rho }_g{h}_{\mathit{lg}}{\kappa }_l}{8\sigma T_s}$ , $F=\frac{h_{\mathit{cb}}{\mathit{Pr}}_l^2}{2B}$ , $\mathrm{\Delta }{T}_{\mathit{sat}}=T_{\mathit{onb}}-T_s$ is the wall superheat, and $\mathrm{\Delta }{T}_{\mathit{sub}}=T_s-T_l$ is liquid subcooling temperature.

When the wall superheat exceeds $\mathrm{\Delta }{T}_{\mathit{sat}}=T_{\mathit{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 ${\dot{q}}_{\mathit{chf}}$ .

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

For the pool boiling value of ${\dot{q}}_{\mathit{chf},0}$ , we used Zuber [58] correlation

$${\dot{q}}_{\mathit{chf},0}=\frac{\pi }{24}{h}_{\mathit{lg}}{\rho }_g{\left(\frac{\mathit{\sigma g}\left({\rho }_l-{\rho }_g\right)}{{\rho }_g^2}\right)}^{1/4}{\left(\frac{{\rho }_l}{{\rho }_l+{\rho }_g}\right)}^{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].

$${\dot{q}}_{\mathit{chf}}={\dot{q}}_{\mathit{chf},0}\left({\alpha }_{\mathit{cr}}-\alpha \right)\left(1+a_1\left(\frac{{\rho }_l{c}_l\mathrm{\Delta }{T}_{\mathit{sub}}}{{\rho }_g{h}_{\mathit{lg}}}\right)\right.+\left.a_2{\mathit{Re}}_l+a_3\sqrt{\frac{{\mathit{Re}}_l{\rho }_l{c}_l\mathrm{\Delta }{T}_{\mathit{sub}}}{{\rho }_g{h}_{\mathit{lg}}}}\right),$$ E24

where ${\alpha }_{\mathit{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\times {10}^{-3}$ [57], and ${\alpha }_{\mathit{cr}}=0.96$ [55] for water.

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

The temperature $T_{\mathit{chf}}$ for the critical heat flux was estimated in this work using the approach proposed by Theler and Freis [67].

$$T_{\mathit{chf}}=\frac{T_s}{1-\frac{T_s{R}_g}{h_{\mathit{lg}}}\mathit{log}\left(2k_g+1\right)},$$ 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_{\mathit{chf}}$ was further corrected for boiling flows using Iloeje-type correlations [66] similar to the one applied $\mathrm{\Delta }{T}_{\mathit{mfb}}$ , see below.

When wall superheat exceeds $\mathrm{\Delta }{T}_{\mathit{chf}}=T_{\mathit{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_{\mathit{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 $\mathrm{\Delta }{T}_{\mathit{mfb}}=T_{\mathit{mfb}}-T_s$ was estimated by Berenson as [68, 69].

$$\begin{array}{l}\mathrm{\Delta }{T}_{\mathit{mfb},0}=0.127\frac{{\rho }_g{h}_{\mathit{lg}}}{{\kappa }_g}\times {\left[\frac{g\left({\rho }_l-{\rho }_g\right)}{{\rho }_l+{\rho }_g}\right]}^{2/3}{\left[\frac{\sigma }{g\left({\rho }_l-{\rho }_g\right)}\right]}^{1/2}{\left[\frac{{\mu }_g}{\left({\rho }_l-{\rho }_g\right)}\right]}^{1/3}\end{array}$$ E26

Iloeje [55, 66] has corrected Berenson equation to take into account the dependence of the $\mathrm{\Delta }{T}_{\mathit{mfb}}$ on the quality and mass flux of the boiling flows in the form

$$\mathrm{\Delta }{T}_{\mathit{mfb}}=c_1\mathrm{\Delta }{T}_{\mathit{mfb},0}\left(1-c_2{\chi }_e^{c_3}\right)\left(1+c_4{G}^{c_5}\right),$$ E27

where ${\chi }_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\times {10}^{-3}$ [57], and ${\alpha }_{\mathit{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_{\mathit{f}}=c_1\cdot {\left[\frac{g{\rho }_g{\kappa }_g^2\left({\rho }_l-{\rho }_g\right){\tilde{h}}_{\mathit{lg}}{c}_{\mathit{pg}}}{D\left(T_w-T_{\mathit{spt}}\right){\mathit{Pr}}_g}\right]}^{0.25}\left(1-c_2{\chi }_e^{c_3}\right)\left(1+c_4{G}^{c_5}\right)$$ 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

$${\dot{q}}_{\mathit{mfb}}=h_{\mathit{f}}\mathrm{\Delta }{T}_{\mathit{mf}}.$$ E29

In the region of single-phase gas flow, the heat transfer coefficient is given by $h_{\mathit{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 $\mathrm{\Delta }{T}_{\mathit{onb}}$ to $\mathrm{\Delta }{T}_{\mathit{chf}}$ , the heat flux was defined using the following interpolation [6]:

$${\dot{q}}_{\mathit{nb}}=y^n{\dot{q}}_{\mathit{onb}}+\left(1-y^n\right){\dot{q}}_{\mathit{chf}},$$ E30

where $n$ is constant, $y$ is defined as $\left(T_w-T_{\mathit{onb}}\right)/\left(T_{\mathit{chf}}-T_{\mathit{onb}}\right)$ , while ${\dot{q}}_{\mathit{chf}}$ and ${\dot{q}}_{\mathit{onb}}$ are given by Eqs. (24) and (22), respectively.

Similarly, for the regime where the wall superheat is increasing from $\mathrm{\Delta }{T}_{\mathit{chf}}$ to $\mathrm{\Delta }{T}_{\mathit{mfb}}$ , the transition boiling heat flux was interpolated as [6].

$${\dot{q}}_{\mathit{tb}}=f_{\mathit{tb}}\cdot {\dot{q}}_{\mathit{chf}}+\left(1-f_{\mathit{tb}}\right){\dot{q}}_{\mathit{mfb}},$$ E31

where $f_{\mathit{tb}}={\left(\frac{T_w-T_{\mathit{mfb}}}{T_{\mathit{chf}}-T_{\mathit{mfb}}}\right)}^2$ , while $T_{\mathit{chf}}$ , $T_{\mathit{mfb}}$ , and ${\dot{q}}_{\mathit{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 ${\tau }_{\mathit{wl}}=0.5f_{\mathit{wl}}{\rho }_l{u}_l^2$ and ${\tau }_{\mathit{wg}}=0.5f_{\mathit{wg}}{\rho }_g{u}_g^2$ with the friction factors $f_{\mathit{wl}\left(g\right)}$ for turbulent and laminar flow given by Churchill approximation [70].

The two-phase pressure drop ${\left(\frac{\mathit{dp}}{\mathit{dz}}\right)}_{2\varphi }$ is given by the Lockhart-Martinelli correlations [71]. The pressure losses are partitioned between the phases as follows [7]:

$$\begin{array}{l}{\tau }_{\mathit{wg}}{l}_{\mathit{wg}}={\alpha }_g{\left(\frac{\mathit{dp}}{\mathit{dz}}\right)}_{2\varphi }\left(\frac{1}{{\alpha }_g+{\alpha }_l{Z}^2}\right),{\tau }_{\mathit{wl}}{l}_{\mathit{wl}}={\alpha }_l{\left(\frac{\mathit{dp}}{\mathit{dz}}\right)}_{2\varphi }\left(\frac{Z^2}{{\alpha }_g+{\alpha }_l{Z}^2}\right)\end{array}$$ E32

where $Z^2$ is given by

$$Z^2=\left(f_{\mathit{wl}}{\mathit{Re}}_l{\rho }_l{u}_l^2\frac{{\alpha }_{\mathit{wl}}}{{\alpha }_l}\right)/\left(f_{\mathit{wg}}{\mathit{Re}}_g{\rho }_g{u}_g^2\frac{{\alpha }_{\mathit{wg}}}{{\alpha }_g}\right),$$ E33

Friction factors $f_{\mathit{wg}\left(l\right)}$ are the same as in above, while coefficients ${\alpha }_{\mathit{wl}}$ and ${\alpha }_{\mathit{wg}}$ depend on the flow pattern [7].

The interface drag is given by ${\tau }_{\mathit{ig}}=-{\tau }_{\mathit{il}}=\frac{1}{2}{C}_D{\rho }_g\left|u_g-u_l\right|\left(u_g-u_l\right)$ , 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 $\sim$ 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 $\sim$ 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 $\sim$ 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.

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