Numerical characterization of underexpanded cryogenic hydrogen gas jets

High-resolution direct numerical simulations are conducted for under-expanded cryogenic hydrogen gas jets to characterize the nearfield flow physics. The basic flow features and jet dynamics are analyzed in detail, revealing the existence of four stages during early jet development, namely, (a) initial penetration, (b) establishment of near-nozzle expansion, (c) formation of downstream compression, and (d) wave propagation. Complex acoustic waves are formed around the under-expanded jets. The jet expansion can also lead to conditions for local liquefaction from the pressurized cryogenic hydrogen gas release. A series of simulations are conducted with systematically varied nozzle pressure ratios and systematically changed exit diameters. The acoustic waves around the jets are found to waken with the decrease in the nozzle pressure ratio. The increase in the nozzle pressure ratio is found to accelerate hydrogen dispersion and widen the regions with hydrogen liquefaction potential. The increase in the nozzle exit diameter also widens the region with hydrogen liquefaction potential but slows down the evolution of the flow structures. s


I. INTRODUCTION
Hydrogen is being increasingly used as an alternative fuel for vehicles and aircraft propulsion. In order to facilitate large scale applications, hydrogen is likely to be delivered and stored in the liquid state (LH2) at cryogenic temperature under pressurized conditions. In the event of potential LH2 leaks, choked conditions will develop at the leaks, resulting in the formation of under-expanded hydrogen jets. The choked jets will typically experience flash evaporation, cryogenic boiling, and condensation/freezing of surrounding gases, accompanied by intense heat and mass transfer. The low temperature brings challenges to characterize cryogenic hydrogen jets due to the complex underlying physics.
Since 2011, various experimental and numerical studies have been conducted to understand the overall dispersion and concentration fields following the LH2 release. Experimental studies on underexpanded cryogenic hydrogen jets were first conducted by Veser et al. 1 They measured hydrogen concentration along the streamwise centerline for the inflow temperature from 35 K onward. It was found that hydrogen concentration in the 35 K case decayed slower than that of the warmer jets released at 80 K and 290 K. Friedrich et al. 2 conducted a series of experimental investigations with cryogenic hydrogen at temperatures of 34 K-65 K to study the behavior of cryogenic jets. They proposed a correlation for predicting the dependence of hydrogen concentration on the nozzle diameters and the cryogenic reservoir conditions. Hecht and Panda 3 measured the concentration and temperature fields of cryogenic hydrogen jets by using Raman scattering at laboratory scales. The effects of diameter, temperature, and pressure of nozzles were investigated. They found that the centerline mass fraction of hydrogen decayed at a similar rate with that under room temperatures, but the half-width of the hydrogen mass fraction decayed more slowly. Limited numerical investigations have also been carried out for cryogenic hydrogen jets. Venetsanos and Giannissi 4 numerically studied the pressurized hydrogen jets at cryogenic temperatures tested by Friedrich et al. 2 using the Reynolds-averaged Navier-Stokes (RANS) approach with the k-ε turbulence model and the ideal gas assumption for hydrogen. The notional nozzle approach was used to bypass the complex flow structures near the nozzle for computational efficiency. The predicted axial hydrogen concentrations were in reasonably good agreement with the measurements.
Considering the lack of studies on cryogenic jets, a brief review is included here for the previous studies of the under-expanded ARTICLE scitation.org/journal/adv hydrogen jets released at environmental temperatures. The nearfield zone of the under-expanded jet exhibits complex flow features including the mixing layer, turbulence transition, shock waves, shock-vortex interaction, etc., and has been the subject of numerous experimental, numerical, and theoretical investigations. [5][6][7][8][9][10][11] The flow characteristics from a nozzle are highly dependent on the ratio of the upstream total pressure, P 0 , to the ambient static pressure Pa, specifically the nozzle pressure ratio (NPR), NPR = P 0 /Pa. Depending on the nozzle pressure ratio, 12 the jets are generally characterized as follows: (1) 13 numerically investigated the effects of ambient pressure on the instantaneous evolution of the nearfield shock wave structures and tip vortices of hydrogen and methane jets. The increasing release pressure was found to result in faster initiation of the Mach disc and promote the radial penetration as well as the mixing of chemical species. They also studied the influence of the ambient temperature 14 and found that the increase in temperature results in the increase in the tip penetration and volumetric growth of jets. The formation of the preliminary vortex rings is attributed to the flow separation at the edge of the embedded shock. Tang et al. 15,16 focused on the starting evolution of a highly under-expanded jet and the Mack disk stabilization process. They observed the strong variation of pressure waves during the extremely high-pressure hydrogen jets. Hamzehloo and Aleiferis 17 studied the effect of the nozzle exit geometry on the under-expanded jet behavior. Their predictions revealed that the decrease in the nozzle length-to-diameter ratio would lead to a longer jet and the reduction in the length of the straight nozzle results in a stronger expansion. The diverging conical nozzle contributes to a higher penetration length than the straight nozzle with the same NPR.
In the present study, the detailed nearfield flow structures and transient physics of the under-expanded cryogenic hydrogen jets are analyzed for the first time using high-resolution direct numerical simulations (DNSs). The numerical setup mimics the experimental conditions of Hecht and Panda. 3 The influences of the nozzle pressure ratio and exit diameter will be analyzed. The predicted local liquefaction due to the complex wave interactions during jet evolution will also be discussed. The present study is focused on the near-field under-expanded region to gain insight on the complex underlying physics. As the measurements of Hecht and Panda 3 only covered the far field, no direct comparison will be made with the experimental data.

A. Governing equations and numerical methods
Two-dimensional multicomponent compressible Navier-Stokes equations are directly solved to model the under-expanded jet dynamics, where ρ, ui, P, and T are the density, velocity, pressure, and temperature of the gas mixture, respectively. Y k is the mass fraction of the species k. The state equation for the ideal gas with multispecies is used to close the gas-phase equations. The assumption of the Fourier heat conduction and Fickian mass diffusion is utilized to consider the molecular contributions in the viscous terms. The Soret and Dufour effects are neglected. δij is the Kronecker delta function, and τij is the Newtonian viscous stress tensor, Sij is the strain tensor. μ is the shear viscosity. V k is the diffusion velocity of the species k and is approximated by A correction velocity V c j is added in the species transport equations to ensure the global mass conservation, where NS is the total number of species. D k and W k are the molecule weight and the mass diffusion coefficient of the kth species, respectively. et is the total energy, i.e., kinetic energy and internal (containing chemical) energy, which is defined as where c p,k is the specific heat capacity at a constant pressure and h 0 f ,k is the specific chemical formation enthalpy at the reference temperature, T ref . The heat flux qj is where λ is the thermal conductivity of the gas mixture. The transport properties including the viscosity, μ k , the heat conductivity, λ k , and the binary diffusion coefficient, D k , of each chemical species are obtained based on the kinetic theory. 18 In particular, the heat conductivity is calculated by using the modified Eucken model. The dynamic viscosity and the binary diffusion coefficient are computed according to the Chapman-Enskog theory, 19 and the semiempirical expressions proposed by Wake and Wassiljewa 20 are used to calculate the dynamic viscosity and heat conductivity of the mixture. The governing equations of the multi-phase jet flow are solved by utilizing our in-house code, which have previously been applied ARTICLE scitation.org/journal/adv to study a variety of compressible flows. 21,22 The conservative finite difference method is used to discretize the spatial terms of the governing equations on the Cartesian grid. The adaptive centralupwind sixth-order weighted essentially non-oscillatory (WENO-CU6) scheme 23 is used for the convection terms to facilitate the simulations of the main flow with low dissipation and achieve a proper resolution of the flow properties around the shock waves. A sixth-order symmetric compact difference scheme is applied for the viscous diffusion terms. Time-integration is realized by the explicit third order Runge-Kutta method.

B. Computational setup
The recent numerical evaluations of Su et al. 10 revealed that there were only minor differences in the predicted 2D and 3D flow structures, such as the scale of the under-expanded jet, the shock wave location, and the flow field (temperature, density, and velocity). The present study is, therefore, limited to 2D simulations for computational efficiency. As shown in Fig. 1, the computational domain spans a 2D domain with width Lx and length Ly. The subscripts x and y refer to the transverse and streamwise directions. In the present study, Lx equals 20 mm and Ly is 20 mm. The numerical setup mimics the experimental configurations of Hecht and Panda 3 with hydrogen (inflow) and air (coflow) injected from the bottom being represented by the inflow boundary condition. The non-reflected boundary conditions are applied on the left as well as the right boundaries following the study of Rudy and Strikwerda. 24 For the outlet conditions, the parameters are interpolated by assuming first-order derivatives. The thermodynamic properties of hydrogen, such as heat capacity and enthalpy, are calculated according to the Aly-Lee equation, 25 which is valid for a wide range of temperatures.
The inflow conditions are specified from the test conditions of Hecht and Panda 3 and summarized in Table I. A slow coflow air with a streamwise velocity of 0.3 m/s is imposed at the inlet boundary in line with the experimental setup. 3 The static pressure, Pa, and static temperature, Ta, of ambient air are 1 bar and 297 K, respectively. The coflow air is considered as the mixture of nitrogen and oxygen with the mass fraction of 0.77 and 0.23, respectively. The nozzle pressure increased from 3.0 bars to 5.0 bars as in the experiments. 3 As shown in Table I, the respective cases are named as LP, MP, and HP to investigate the effect of nozzle pressure on the jet behavior. Further details of the experimental conditions which were simulated in the numerical simulations can be found in Ref. 3. In addition, for the under-expanded jet of case HP, the nozzle exit diameters are varied from 0.75 mm to 1.5 mm for cases HPD1, HPD2, and HPD3. An under-expanded jet is expected when the nozzle pressure is above 2 bars. Considering the present pressure and temperature of the nozzle, the stagnation conditions are all located in the vapor/gas phase. 3 As the nozzle conditions are away from the saturation curve, the current inlet of the hydrogen jet is purely low temperature gas.

C. Code verification and grid sensitivity study
The calculation of Sod's shock tube problems 26 is simulated first to verify the capability of the code to capture the complex compression and expansion waves during jet under-expansion. The initial condition and the associated distributions of the parameters are As shown in Fig. 2, the predictions are in excellent agreement with the analytical result with the predictions almost completely overlapping with the exact solution. Besides, there is no oscillation of pressure and density around the shock discontinuity. The physical properties of hydrogen are obtained by solving the Navier-Stokes equations, associated with the state equation for an ideal gas with multi-species. In order to verify the treatment of the physical properties of cryogenic hydrogen adopted in the code, the computed density values from the temperatures and pressures at the nozzle exit in the experiments are shown in Table II. There is no difference between the density value in the experiment and that presently calculated for case LP. With the increase in the nozzle pressure ratio, the difference increases to 1.2% for case HP. The effect of such minor differences on the predicted jet flow features is neglectable. The grid independence for the simulations of case HP, which is the under-expanded jet with a nozzle pressure ratio of NPR = 5.0, is tested by three different meshes with uniform grid sizes of Δ = 10 μm, 20 μm, and 30 μm, respectively. Figure 3 shows the instantaneous distributions of the density gradients at t = 30 μs with different grid sizes. The streamwise profiles of pressure and the Mach number along the transverse centerline are shown in Fig. 4. The predictions obtained by using the grid size of Δ = 20 μm are almost the same as that with the grid size of Δ = 10 μm, indicating that further refinement of the grid resolution from Δ = 20 μm has a minimum effect on the predictions. Therefore, uniform Cartesian grids with the grid size of Δ = 20 μm are applied in the subsequent numerical simulations.

III. UNSTEADY JET FEATURES
As mentioned earlier, the present study is focused on the underlying physics in the near-field under-expanded region, which was not reported in the work of Hecht and Panda. 3 Indeed, it is challenging to resolve the under-expanded region in laboratory tests due to the existence of multiple shock waves and steep gradients of the field variables. DNS is, hence, a viable predictive tool to provide such insight. In this section, the DNS results of the unsteady flow structures of case HP with NPR = 5.0 are analyzed in detail. The jet tip penetration, H tip , is one of the key properties for under-expanded gaseous jets, which characterize the spread features of hydrogen. Figure 5 shows the evolution of H tip , which is shown by the scatters. The jet tip penetration is found to have a quasi-linear relationship with time, as depicted by the dashed lines. There appear to be two stages during the initial development of the jet. From the jet start to the time t = 40 μs, the slope of the H tip to time quasi-linear curve is about k = 0.25. From t = 40 μs, the slope is found to decrease with a value of k = 0.21, which indicates that the jet spreads more slowly.
The early stages (from t = 10 μs to t = 40 μs) for the development of the near-nozzle flow structures are shown in Fig. 6 for case HP, which presents the instantaneous distributions of the density gradient, hydrogen mass fraction, and pressure, and Fig. 7 shows the streamwise variations of pressure along the transverse centerline at different times, corresponding to the time regime in Fig. 6. The early jet evolution can be divided into four stages, namely, (a) initial penetration (t = 0 μs-10 μs), (b) establishment of near-nozzle expansion (t = 10 μs-20 μs), (c) formation of downstream compression (t = 20 μs-30 μs), and (d) wave propagation (t = 30 μs-40 μs). The flow physics in each stage is further analyzed as follows: (a) Initial penetration (t = 0 μs-10 μs): A leading shock (LS) is formed from the pressurized hydrogen jets spouting out from the nozzle exit. It spreads downstream, resulting in the sudden pressure increase due to the compression effects, as shown by the pressure distribution in Fig. 6. At the beginning of the jet, the strong LS gradually weakens as it propagates, as shown by the decrease in its pressure in Fig. 7(a). This is attributed to the spatial weakening effects from the quick rise of its volume during the continuous spreading, as shown by the decreasing wave pressure from t = 10 μs to t = 40 μs. Due to the density difference between the hydrogen with cryogenic temperature and the ambient air, a jet interphase with a high-density gradient occurs around the jet boundary. In the upstream regime from the interphase, a strong expansion wave (EW) is generated in the inner region of the jet due to the large density  The jet pressure decreases continuously until reaching the minimum pressure around the streamwise location y = 1.5 mm from the nozzle exit. It is found that the location of the minimum pressure in Fig. 7(b), which is the starting point of the expansion, stays at the same location as the jet advances. The expansion waves reflect from the jet boundary, where the pressure equals the ambient pressure and results in compression waves, which converge toward the inner of the jet flow and coalesce to form the oblique shock waves. These are also referred to as the intercepting shock wave, as indicated by the pressure increase from y = 1.6 mm in Fig. 7(b). A second intercepting shock is generated around y = 2.5 mm at t = 20 μs, resulting in a further pressure increase. On the other hand, the compression waves reflect from the jet boundary, resulting in the intercepting expansion waves facing to the outer jet region, and the pressure decreases accordingly. propagates downstream and becomes stronger, as shown by the increases in the pressure peaks in Fig. 7(c), and a curved shock wave is formed finally, as shown by the density gradient at t = 30 μs. (d) Wave propagation (t = 30 μs-40 μs): The curved shock stabilizes around y = 2.8 mm, and its strength increases continuously until the pressure, P/Pa, reaches a maximum value of 2.6. The LS and EW continue to propagate downstream, and their strengths gradually decrease.
In the above evolution processes, both the flow structures and jet shapes change. From the initial mushroom shape, the jet boundary evolves into a quasi-rectangle shape, as shown by the density gradient in Fig. 6.
The complex pressure wave formed in the cryogenic jet results in some regions where the partial pressure of hydrogen (P H2 ) is higher than its saturated vapor pressure (Pvap), creating the conditions for localized hydrogen liquefaction. The hydrogen liquefaction potentiality (HLP) is calculated to analyze the possibility of local liquefaction as follows: If the HLP is higher than zero, the occurrence of hydrogen liquefaction is expected to occur theoretically. The red dashed lines in Fig. 6 denote the regions where HLP > 0, indicating that localized liquefaction can occur in the under-expanded jets of cold hydrogen gas. The liquefaction tends to occur after the EW with relatively lower local pressures and temperatures, and the region has an unsteady variation due to the spatial evolution of wave structures. In particular, the liquefaction regions expand with the decrease in the post-EW pressure from time t = 30 μs-40 μs.
The downstream development of the cryogenic hydrogen jets is shown in Fig. 8 by the density gradient. The jet tip propagates more slowly during this period. The strong shear results in the rolling and shedding of vortices, which is associated with a large amount of hydrogen entrained by the coflow air. The velocity of the jet tip decreases due to momentum transfer from the hydrogen jet to the coflow air, resulting in the increase in the local velocity. Complex acoustic waves are formed around the jet flow. The potential liquefaction region decreases continuously as the expansion wave weakens.

IV. EFFECTS OF THE NOZZLE PRESSURE RATIO
Further simulations are run to investigate the effect of the nozzle pressure ratio on the nearfield physics. The NPR was systematically varied from 3.0 to 5.0. Figure 9 illustrates the evolution of jet tip penetration as time advances. Generally, all three cases have a quasi-linear increase in jet tip penetration. With nearly the same H tip at time t = 5 μs, the increase in the NPR results in the faster increase in H tip . With the increase in the inject pressure and accordingly higher density, the momentum of the hydrogen jet increases, resulting in higher jet penetration speed.
The early stages of the near-nozzle flow structures for cases LP, MP, and HP from time t = 10 μs-40 μs are shown in Fig. 10 by the density gradients. It is found that the length of the jet root increases with the increase in nozzle pressure ratios. The shape of the jet keeps almost unchanged with a round head during the penetration for the low NPR cases, and there is no hydrogen liquefaction for cases LP and MP. For case LP, the expansion and compression during the jet evolution are weak and the wave structures inside the jet are unobvious, as shown in Fig. 10. The resulting partial pressure of hydrogen is found to be always lower than the saturated vapor pressure of hydrogen, which renders liquefaction impossible.
The further downstream development of the jets with a low NPR is shown in Fig. 11. For case LP with NPR = 3.0, the boundary of the jet root is straight and only weak intercepting shocks are formed from the reflected compression waves. As the NPR increases to 4.0, the boundary of the jet root becomes a convex-to-concave shape with a stronger curved shock generated inside the jet. As the jet flow develops, the smooth round boundary of the jet head is disrupted due to the shearing with the coflow air and the shedding of vortices starts at a more upstream position as the NPR decreases. In addition, the acoustic waves around the jet gradually weaken and tend to disappear at time t = 80 μs for case LP.

V. EFFECTS OF THE NOZZLE EXIT DIAMETER
In order to investigate the effects of the nozzle exit diameter on the near field characteristics and the liquefaction potentiality, DNS was conducted with different exit diameters from 0.75 mm to 1.50 mm. The evolutions of jet tip penetration are shown in Fig. 12. The increase in the nozzle exit diameter results in the acceleration of hydrogen dispersion as jets spouting from larger nozzles have higher initial mass flow rates of hydrogen, and the higher initial momentum leads to the faster penetration. The difference of jet tip penetration between cases HP and HPD1 reduces to very small after 75 μs.
The early stages of the jet development for cases HP, HPD1, HPD2, and HPD3 are shown in Fig. 13. For case HPD1 with a small nozzle exit, the jet shape and its evolution process are almost the same with those of case HP, but the evolution process tends to be quicker. For example, the jet head with a quasi-rectangle shape at time t = 30 μs for case HPD1 is found to be similar to the shape for case HP at t = 40 μs, as denoted by the black lines of density gradients in Fig. 13. This indicates that the flow structures evolve more rapidly with a smaller nozzle exit. However, the penetration speed of the jet tip fails to accelerate at this small nozzle diameter. As the nozzle exit diameter increases to 1.25 mm, the dispersion regime expands in the transverse direction and the flow structures evolve more slowly. The curved shock wave inside the jet root is generated at time t = 40 μs in case HPD2. With a further increase in the nozzle exit diameter to 1.5 mm, the jet shape evolves even more slowly and the mushroom shape of the jet is found to be almost unchanged during the early-stage evolutions of the jet. It can be observed that the HLP regime has an inhomogeneous distribution, as indicated by the red lines at time t = 30 μs. The conditions for liquefaction are found in narrow regions surrounded by the intercepting shock and the jet boundary in the upstream of the EW. The HLP quickly becomes less than zero from t = 40 μs. The unsteady variation of the liquefaction is mainly due to the compression-expansion complex during the under-expanded jet.
The further developments of the jets with different nozzle exit diameters from t = 50 μs-80 μs are shown in Fig. 14. The jet with a small nozzle exit tends to be very unsteady with the shearing between hydrogen jets, and coflow air induces large-scale vortices, indicating that the flow structures for case HPD1 develop more quickly. For case HPD2, the jet tip is smooth with a semicircular boundary at time t = 80 μs. The quasi-rectangle shape of the jet head occurs at time t = 50 μs for case HPD2 and t = 60 μs for case HPD3. The area of liquefaction potentiality increases with the increase in nozzle exit diameter, associated with the increase in initial mass flow rates of hydrogen.
For different nozzle exit diameters, the streamwise pressure variations along the transverse centerline during the jet evolution are shown in Fig. 15. The expansion and compression waves are formed quickly near the nozzle for case HPD1. The expansion of the highpressure hydrogen from the nozzle is achieved at time t = 10 μs. The establishment of the curved shock wave with the pressure peak, P/P 0 , around 2.8 is formed at time t = 30 μs. Subsequently, this shock wave gradually weakens from time t = 30 μs to t = 40 μs as the jet develops. A second compression wave with lower pressures is generated in the downstream region. The expansion wave which induces the local liquefaction gradually weakens with the decrease in the pressure drop, and the HLP regime gradually narrows, as shown in Fig. 13. For case HPD2 with a larger nozzle exit, the expansion and compression waves near the nozzle need a longer time to stabilize. The compression wave reached its pressure peak at time t = 50 μs. The streamwise length for the formation of these waves increases with the increase in the nozzle exit diameter. For the curved shock wave, the streamwise location increases from y = 2.8 mm for case HPD1 to y = 3.9 mm for case HPD2. As the nozzle exit further enlarges for case HPD3, the stabilization of the nearfield flow structures progressed more slowly and the formation of the curved shock is finally achieved around time t = 60 μs. In addition, the variation of the nozzle exit diameter is found to have a minor influence on the strength of the expansion and compression waves for the under-expanded jet.

VI. CONCLUSIONS
High-resolution multi-component direct numerical simulations have been conducted for under-expanded cryogenic hydrogen jets to characterize the near-field flow physics. The DNS results also shed light on the influences of the nozzle pressure ratio and exit diameter on the nearfield jet dynamics. The potential of liquefaction in the cryogenic hydrogen gas jet is analyzed through the dynamics of pressure evolution during the complex expansion-compression processes.
The results indicate that the early evolution of the jet can be divided into four stages, which are (a) initial penetration, (b) establishment of near-nozzle expansion, (c) formation of downstream compression, and (d) wave propagation. The strong expansion formed in the jet head can lead to localized liquefaction due to the difference between the partial pressure of hydrogen and the saturated vapor pressure of hydrogen.
The penetration of the jet is approximately linearly dependent on time. The increase in the nozzle pressure ratio and nozzle exit diameter results in a faster rise of the jet tip penetration. The variation in the nozzle pressure ratio not only affects the hydrogen dispersion but also changes jet shapes. The jet head varies from a round shape for the low nozzle pressure ratio to a quasi-rectangle shape for the high nozzle pressure ratio. The increase in the nozzle exit diameter leads to a slow evolution of the flow structures and delays the stabilization of the waves further downstream. The region with positive hydrogen liquefaction potentiality expands with the increasing nozzle exit size.