Modeling and Analysis of the Flow Characteristics of Liquid Hydrogen in a Pipe Suffering from External Transient Impact

Abstract

Pipes can be subjected to external transient impacts such as accidental collision, which affects the safe operation of storage and transportation systems for liquid hydrogen. Fluid–structure coupling calculation for a pipe under external transient impact is performed, and the flow characteristics of liquid hydrogen in the pipe are analyzed. The pipe deforms and vibrates when suffering from external transient impact. Liquid hydrogen pressure in a cross-section plane increases along the pipe deformation direction. Additionally, external transient impact enhances the disturbance of liquid hydrogen near the pipe wall. The increased flow resistance and the energy induced by the deformed pipe both affect the flow of liquid hydrogen, and contribute to the fluctuated characteristics of liquid pressure drop. In addition, the phase state of liquid hydrogen remains unchanged in the pipe, indicating that little of the induced energy is transformed into the internal energy of liquid hydrogen. The work provides theoretical guidance for the safe operation of liquid hydrogen storage and transportation systems.

1. Introduction

Liquid hydrogen, which has been widely used as space cryogenic propellant [1], has the advantages of high gravimetric and volumetric energy densities [2], environmental friendliness and renewability. Moreover, hydrogen energy is considered as a future energy and exhibits promising application potentials in fuel cells and vehicles [3,4]. The flow characteristics of liquid hydrogen in pipes are vital for the safe and stable operation of storage and transportation systems. However, the pipe could be subjected to external transient impacts such as accidental collision, which affects the dynamics behavior of the pipe and the flow characteristics of liquid hydrogen in the pipe. In particular, when the external transient impact exceeds a certain strength, the pipe can even break and fail, resulting in the leakage of liquid hydrogen into the surrounding environment, and the safety of personnel and equipment is seriously threatened [5].

In past decades, great efforts have been dedicated to investigations regarding the leakage and dispersion of liquid hydrogen, caused by the rupture of pipes or storage vessels. The National Aeronautics and Space Administration (NASA) in the USA [6], the Federal Institute for Materials Research and Testing (BAM) in Germany [7] and the Health and Safety Laboratory (HSL) in the UK [8] have carried out liquid hydrogen spill experiments in open environments, between buildings and under pipe failure conditions, respectively. Verfondern et al. [9] conducted liquid hydrogen pool spread experiments on the surface of water and solid ground. Due to the huge cost and the safety issues related to the experiment, many researchers also investigated the leakage and dispersion of liquid hydrogen numerically [10,11,12,13,14,15,16,17,18,19,20,21], which revealed the underlying mechanisms and the effects of various conditions, such as spill source, season, wind flow, air humidity, etc.

However, there is little research focused on the flow characteristics of liquid hydrogen in a pipe subjected to an external transient impact. Instead, the flow and heat transfer of fluid in a pipe under vibration excitation has been reported in the literature. Liu et al. [22] experimentally studied the influence of sinusoidal vibration on the heat transfer characteristics of circular pipes and found that vibration could significantly affect the heat transfer process. Tian et al. [23] applied lateral vibration on a pipe and found that the wall heat transfer was significantly strengthened due to the rapid growth of the thermal boundary layer. For cryogenic fluid, Zheng et al. [24] numerically studied the influence of wall vibration on the boiling flow and heat transfer of liquid hydrogen in a horizontal pipe. The vibration of a circular pipe could enhance the heat transfer intensity between liquid hydrogen and the wall, and the enhancement is more obvious at a lower inlet velocity. Chen et al. [25] studied the influence of vibration on the filling process of liquid hydrogen in a horizontal circular pipe. The increase in vibration frequency and amplitude enhances the heat exchange between liquid hydrogen and the pipe, while the vibration destroys the stable boundary of two-phase flow and promotes the bubble separation, which affects the stability of the transportation system. The external transient impact and the vibration could both be classified as the external excitation. However, compared with vibration, the transient impact commonly has an extremely short duration and high intensity, which will cause the deformation of the pipe and energy input to the liquid flow.

Consequently, it is essential to carry out research on the flow characteristics of liquid hydrogen in a pipe suffering from an external transient impact. In the present work, a three-dimensional computational fluid dynamics (CFD) model predicting the flow of liquid hydrogen in a pipe is first designed and validated. Then, the fluid–structure coupling calculation for the pipe under external transient impact is performed in ANSYS Workbench 14.5 (ANSYS, Inc., Canonsburg, PA, USA). Finally, the flow characteristics of liquid hydrogen in the pipe under external transient impact are analyzed. The work provides theoretical guidance for the safe operation of liquid hydrogen storage and transportation systems, which could help to promote the smooth progress of space launch missions and the applications of hydrogen in fuel cells and vehicles, as well as ensuring the safety of personnel and equipment.

2. CFD Model Predicting the Flow of Liquid Hydrogen in a Pipe

2.1. Mathematical Model

The mixture multiphase model is adopted to predict the flow of liquid hydrogen in a pipe, and the conservation equations for mixture mass, mixture momentum and mixture enthalpy are:

$$\frac{\partial {\rho }_{\mathrm{m}}}{\partial t}+\nabla \cdot ({\rho }_{\mathrm{m}}{\overrightarrow{v}}_{\mathrm{m}})=0,$$

(1)

$$\frac{\partial }{\partial t}({\rho }_{\mathrm{m}}\overrightarrow{v_{\mathrm{m}}})+\nabla \cdot ({\rho }_{\mathrm{m}}\overrightarrow{v_{\mathrm{m}}}\overrightarrow{v_{\mathrm{m}}})=-\nabla p+\nabla \cdot \left[{\mu }_{\mathrm{m}}\left(\nabla \overrightarrow{v_{\mathrm{m}}}+\nabla {\overrightarrow{v_{\mathrm{m}}}}^{\mathrm{T}}\right)\right]+{\rho }_{\mathrm{m}}g+\nabla \cdot ({\displaystyle \sum _{i=1}^2{\alpha }_i{\rho }_i}{\overrightarrow{v}}_{\mathrm{dr},i}{\overrightarrow{v}}_{\mathrm{dr},i}),$$

(2)

$$\frac{\partial }{\partial t}{\displaystyle \sum _{i=1}^n({\alpha }_i{\rho }_i{h}_i)}+\nabla \cdot {\displaystyle \sum _{i=1}^n\left({\alpha }_i{\overrightarrow{v}}_{{}_i}\left({\rho }_i{h}_i+p\right)\right)}=\nabla \cdot \left(k_{\mathrm{eff}}\nabla T\right),$$

(3)

where ρ, v and μ are density, velocity and kinetic viscosity, respectively. k_eff is the effective conductivity coefficient, v_dr,i is the drift velocity for phase i and α is the volume fraction. The subscripts m, l and g represent the mixture phase, the liquid phase and the vapor phase, respectively. The physical properties of hydrogen and liquid hydrogen are calculated in terms of temperature in the present work, and the corresponding data are acquired from NIST REFPROP (National Institute of Standards and Technology, Gaithersburg, MD, USA) [26]. The mixture properties are calculated by the following equations:

$${\rho }_{\mathrm{m}}={\alpha }_{\mathrm{l}}{\rho }_{\mathrm{l}}+{\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}},$$

(4)

$$\overrightarrow{v_{\mathrm{m}}}=({\alpha }_{\mathrm{l}}{\rho }_{\mathrm{l}}\overrightarrow{v_{\mathrm{l}}}+{\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}\overrightarrow{v_{\mathrm{g}}})/{\rho }_{\mathrm{m}},$$

(5)

$${\mu }_{\mathrm{m}}={\alpha }_{\mathrm{l}}{\mu }_{\mathrm{l}}+{\alpha }_{\mathrm{g}}{\mu }_{\mathrm{g}}.$$

(6)

The phases move with different velocities, and the slip velocity is expressed as [27]:

$${\overrightarrow{v}}_{\lg }=\frac{{\tau }_{\mathrm{l}}}{f_{drag}}\frac{{\rho }_{\mathrm{l}}-{\rho }_{\mathrm{m}}}{{\rho }_{\mathrm{l}}}\overrightarrow{a},$$

(7)

where the drag force is

$$f_{drag}=\{\begin{array}{l}1+0.15R{e}^{0.687}{,}^{}{{}^{}}^{}Re\le 1000\\ 0.0183Re{,}^{}{{}^{}}^{}Re>1000\end{array}.$$

(8)

The Lee model [28] is used to simulate the liquid–vapor phase change and mass transfer process for hydrogen:

$$\frac{\partial }{\partial t}({\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}})+\nabla \cdot ({\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}\overrightarrow{v_{\mathrm{g}}})={\dot{m}}_{\lg }-{\dot{m}}_{\mathrm{gl}},$$

(9)

$$\{\begin{array}{l}{\dot{m}}_{\lg }={\gamma }_{\mathrm{ev}}\cdot {\alpha }_{\mathrm{l}}{\rho }_{\mathrm{l}}(T_{\mathrm{l}}-T_{\mathrm{sat}})/T_{\mathrm{sat}}{,}^{}i{f}^{}{T}_{\mathrm{l}}>T_{\mathrm{sat}}\\ {\dot{m}}_{\mathrm{gl}}={\gamma }_{\mathrm{con}}\cdot {\alpha }_{\mathrm{g}}{\rho }_{\mathrm{g}}(T_{\mathrm{g}}-T_{\mathrm{sat}})/T_{\mathrm{sat}}{,}^{}i{f}^{}{T}_{\mathrm{g}}<T_{\mathrm{sat}}\end{array},$$

(10)

where ${\dot{m}}_{\lg }$ and ${\dot{m}}_{\mathrm{gl}}$ are the mass transfer rates due to evaporation and condensation, respectively. The coefficients γ_ev and γ_con are specified as 0.25 s⁻¹ for the phase change of hydrogen [18]. The liquid phase and the vapor phase are in thermal equilibrium state with the same temperature.

The realizable k-ε model [29] is adopted for turbulent closure:

$$\frac{\partial }{\partial t}(\rho k)+\frac{\partial }{\partial x_j}(\rho k{u}_j)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\cdot \frac{\partial k}{\partial x_j}\right]+G_k+G_b-\rho \epsilon ,$$

(11)

$$\frac{\partial }{\partial t}(\rho \epsilon )+\frac{\partial }{\partial x_j}(\rho \epsilon u_j)=\frac{\partial }{\partial x_j}\left[\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\cdot \frac{\partial \epsilon }{\partial x_j}\right]+\rho C_{1}E\epsilon -\rho C_2\frac{{\epsilon }^2}{k+\sqrt{\nu \epsilon }}+C_{1\epsilon }\frac{\epsilon }{k}{C}_{3\epsilon }{G}_b,$$

(12)

where $C_1=\max \left[0.43,\frac{\eta }{\eta +5}\right]$, $\eta =\sqrt{2S_{ij}\cdot S_{ji}}\cdot \frac{k}{\epsilon }$, $S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$. G_k and G_b are the generation of turbulence kinetic energy due to the mean velocity gradients and the buoyancy, respectively. In addition, standard wall functions are adopted for the near-wall treatment.

2.2. Model Validation

The experimental data from Tatsumoto et al. [30] for liquid hydrogen flow in a pipe under low wall superheat conditions are used to validate the reliability of the developed model. Figure 1 illustrates the calculation domain for the flow of liquid hydrogen in a horizontal pipe schematically. The inner diameter of the pipe is 3 mm. The lengths of walls labeled with A and B are 75 mm and 100 mm, respectively. The wall labeled A is adiabatic. The wall labeled B is heated by a heater, and a temperature difference exists between the pipe wall and liquid hydrogen. Liquid hydrogen in saturated state (with the pressure and temperature of 0.7 MPa and 29 K) enters the pipe from the left plane, while the outlet (the right plane) is specified as out-flow boundary condition. The flow is composed of liquid hydrogen only at the entrance. The velocities of liquid hydrogen are 4.75 m/s, 16.7 m/s and 31.5 m/s, respectively.

Figure 2 shows the comparison between the simulation results and the experiment, where ΔT is the temperature difference between the pipe wall and the liquid hydrogen. It can be seen that the simulation result and the experimental data are basically consistent, with the maximum relative error less than 15%. The complex two-phase flow and phase change behaviors of liquid hydrogen, and the negligible sensor errors in the experiment, contribute to the deviation between the simulation and experiment. It is considered that the developed model in the paper could be used to simulate the flow and heat transfer of liquid hydrogen in the pipe under low wall superheat conditions.

3. Calculation Strategy for the Fluid–Structure Coupling of the Pipe under External Transient Impact

3.1. Physical Model

Pipes for transportation of cryogenic fluid commonly have a double-layer structure, and there is a thermal insulation layer outside the pipe, so as to reduce the heat exchange with the surrounding environment. The outer wall of the pipe is first subjected to the external transient impact and the impact strength is reduced. However, if the external transient impact is strong enough, the double-layer structure of the pipe can be ruptured and fail. In the present work, the pipe is simplified as a single-layer structure under the transient impact which has been attenuated by the outer layer. The inner diameter, wall thickness and length of the pipe are 100 mm, 2.5 mm and 3 m, respectively. The pipe material is austenitic stainless steel. The density, Young’s modulus, Poisson ratio, yield strength and tensile strength of the selected material are, respectively, 7860 kg/m³, 210 GPa, 0.279, 1615 MPa and 1950 MPa [31].

A physical model for the fluid–structure coupling is illustrated in Figure 3. For the solid domain, SOLID187 is used for the finite element modeling. SOLID187 is defined by 10 nodes, and has three degrees of freedom at each node. It has plasticity, hyper elasticity, creep, stress stiffening, large deflection and large strain capabilities [32]. Both the entrance (left plane) and the exit (right plane) of the pipe are set as fixed support, which are prevented from moving and deforming. The surface of the pipe is subjected to atmospheric pressure (101,325 Pa), and the center of the pipe is subjected to the external transient impact which is perpendicular to the surface. The commonly used external impact curve includes triangular wave, semi-sinusoidal wave, etc. [33], with the duration in milliseconds. In fact, the semi-sinusoidal wave and the triangular wave are approaching, especially when the duration is short and the intensity is high. In the present work, the external transient impact is characterized by a triangular wave, with a duration of 20 ms and a peak of 10,000 N, as shown in Figure 4. It resembles the impact by a falling rock with a mass of 10 kg and a velocity of 7 m/s (the amplification factor = 2, and normal recovery coefficient = 0.42) [34]. Liquid hydrogen enters the pipe with a velocity of 8 m/s and a temperature of 20.35 K, while the outlet pressure is set as 101,325 Pa. The internal pressure of the pipe is above 101,325 Pa which varies with location due to the flow resistance, and liquid hydrogen is slightly subcooled in the pipe. The viscosity of liquid hydrogen at the outlet is 1.33×10⁻⁵ Pa·s [26], and the Reynolds number is 60,150, which indicates turbulent flow. As the heat loss of a well-insulated cryogenic pipe is commonly less than 1 W/m [35], the leaked heat energy through the pipe is relatively low, especially when compared with the energy induced by the external transient impact, and the wall is considered adiabatic and heat transfer with the environment is neglected. In addition, the wall is non-slip.

3.2. Calculation Method

The fluid–structure coupling calculation is performed in ANSYS Workbench 14.5. The fluid flow module is solved by ANSYS Fluent, while the transient force module is solved by Mechanical APDL (ANSYS, Inc., Canonsburg, PA, USA). The data are transmitted between the two modules through the fluid–structure coupling surface (i.e., pipe wall adjacent to the fluid) in the system coupling module. The flow of liquid hydrogen is calculated firstly, and the fluid pressure is transmitted to the transient analysis module. Then, the transient force calculation is performed and the displacement data are transformed to the flow calculation module. The convergence target for the data transfer control is 10⁻³. The step size is 0.001 s. The external transient impact is applied at 0.7 s when the flow is already steady, and the total simulation time is 1.3 s. The time in the following analysis is set as 0 when the external transient impact is applied. Specifically, for fluid flow calculation, the PISO algorithm is chosen for pressure–velocity coupling, and the standard scheme is used for pressure. A second-order upwind scheme is adopted to solve the momentum, volume fraction, turbulent kinetic energy, turbulent dissipation rate and energy equations. The standard initialization method is chosen. Time step for the flow calculation is also 0.001 s, and the convergence criterion is 10⁻⁴.

3.3. Computational Mesh and Mesh Independency Test

Computational meshes for the numerical analysis are shown in Figure 5 and Figure 6. For the solid domain, mesh is locally refined near the impact region. For the fluid domain, unstructured mesh is generated, and is also locally refined near the impact region, with a minimum mesh size of 0.0004 m, a maximum size of 0.006 m and a growth rate of 1.20. The moving-grid model is used for flow calculation due to the deformation of the calculation domain. Smoothing and remeshing methods are adopted for the dynamic mesh with implicit updating. The entrance and outlet of the pipe are stationary, the interior-fluid zone is deforming, and the pipe wall is system coupling.

Mesh independency tests are conducted to ensure that the numerical results are independent of mesh size, and the results are illustrated in Figure 7. The vertical deformation of the impact point along the -z direction at 0.01 s and the pressure drop between the entrance and the exit of the pipe at 0.01 s are the variables selected for comparison. It can be seen that the relative deviation of the pipe vertical deformation calculated at 286,562 cells from that at 372,865 cells (for the solid domain) is 1.31%, and the relative deviation of the pressure drop calculated at 683,121 cells from that at 872,563 cells (for the fluid domain) is −1.05%. Considering the mesh independency and the computation efficiency, the mesh with 286,562 cells for the solid domain and the mesh with 683,121 cells for the fluid domain are chosen.

4. Flow Characteristics of Liquid Hydrogen in the Pipe under External Transient Impact

Dynamics behavior of the pipe which influences the flow of liquid hydrogen is firstly observed. The total deformation of the pipe at 0.01 s with liquid hydrogen inside is shown in Figure 8. The vertical deformations of the impact point with and without liquid hydrogen in the pipe are shown in Figure 9. The time “0” corresponds to the time when the impact is applied and not to the actual simulation time. It can be seen that the pipe will deform and vibrate when suffering from external transient impact. Compared with the condition of no fluid in the pipe, the maximum amplitude decreases slightly from 9.78 mm to 9.70 mm when there is fluid. Additionally, with fluid in the pipe, the vibration amplitude of the pipe attenuates rapidly, and the positive amplitude attenuates to 6.9% of the maximum amplitude 0.2 s after the transient impact. The pipe vibration becomes unobvious 0.3 s after the transient impact. On the contrary, when no fluid exists in the pipe, the amplitude of pipe vibration decreases slowly with time and, even 0.8 s after the impact, it has not attenuated. It can be concluded that energy induced by the external transient impact is dissipated by both the pipe and the fluid. Moreover, the energy dissipation by the fluid is higher than that by the pipe.

The contours of the pressure distribution of liquid hydrogen on the vertical plane (y = 0) are shown in Figure 10, and the pressure distribution contours on the cross-section plane (x = 1.5 m) are shown in Figure 11. Before the pipe is subjected to the external transient impact (i.e., t = 0 s), the pressure of liquid hydrogen decreases uniformly along the flow direction, with the pressure in the cross-sectional direction being nearly constant. After suffering from the external impact, the pipe is deformed and consequently works on the liquid hydrogen. From Figure 10 and Figure 11, it can be seen that partial energy induced by the external transient impact is converted into the pressure energy of liquid hydrogen. The pressure distribution of liquid hydrogen in the cross-section direction is no longer constant. Instead, liquid pressure increases along the deformation direction. For example, the pipe deforms in the -z direction at 0.01 s, and the liquid pressure increases along the same direction. As shown in Figure 9a, the deformation amplitude of the pipe decreases with time, and the pressure distribution characteristics of liquid hydrogen gradually return to the state before the external transient impact is applied. Though the external impact is transient, the pressure distribution is influenced for approximately 0.3 s.

Variations in the pressure drop between the entrance and exit of the pipe with time are shown in Figure 12. When the pipe is subjected to the external transient impact, the pressure drop of liquid hydrogen in the pipe presents fluctuating characteristics. Figure 13 illustrates the turbulent kinetic energy contours near the pipe wall at the time of 0 and 0.01 s, respectively. It is apparent that the external transient impact enhances the disturbance of liquid hydrogen near the wall of the pipe. In addition, the changes in the flow cross-sectional area and flow direction increase the resistance of the liquid flow. On the other hand, the deformed pipe wall works on the fluid and increases the pressure energy of fluid. The increased flow resistance and the energy induced by the deformed pipe both affect the flow of liquid hydrogen, and contribute to the fluctuating characteristics of liquid pressure drop. In addition, the variation amplitude of the pressure drop gradually decreases with time and approaches the initial state.

The phase state of liquid hydrogen in the pipe is also investigated. The volume fraction of vapor phase in the pipe remains 0 even after the external transient impact, as shown in Figure 14. It should be noted that liquid hydrogen is saturated at the exit and subcooled inside the pipe, and the wall is adiabatic. The result indicates that little energy of the external transient impact is transformed into the internal energy of liquid hydrogen. Instead, the energy is dissipated by the fixed support constraints at both ends of the pipe, and transformed into the pressure energy and kinetic energy of the fluid. However, if there is heat loss through the pipe wall, the strengthened liquid turbulence will enhance the heat transfer between liquid hydrogen and the wall [24], which could possibly induce two-phase flow in the pipe, threatening the safe operation of the liquid transportation system. Specifically, the vapor phase in the pipe could cause pressure oscillation and pipe vibration. A geyser may occur and lead to the failure of the pipe [36]. In addition, the normal workings of equipment in the transportation system (e.g., flowmeter and pump) are influenced.

In the present work, the pipe is not ruptured during the impact. It is apparent that the action region of the external transient impact suffers from the strongest stress intensity. Establishing a linearization path at the shortest distance along the wall thickness direction, the membrane strength is calculated to be 98.5 MPa and the membrane plus bending stress is calculated to be 433.2 MPa at the time of 0.01 s, which are within the designed stress intensity. The study predicts and reveals the flow characteristics of liquid hydrogen in a pipe suffering from external transient impact, which could provide theoretical guidance for formulating the emergency plan, identifying the potential accident by monitoring the fluid flow parameters, etc. and thus the safe operation of liquid hydrogen storage and transportation systems. In view of theoretical guidance for the safe handling of liquid hydrogen, mechanical analysis of the pipe under external transient impact, the flow of liquid hydrogen in a pipe under pipe rupture conditions, etc. need to be investigated further.

5. Conclusions

A three-dimensional CFD model predicting the flow of liquid hydrogen in a pipe is designed and validated, and the fluid–structure coupling calculation for the pipe under external transient impact is performed. The flow characteristics of liquid hydrogen in the pipe under external transient impact are then numerically analyzed. The main conclusions are as follows:

The mixture multiphase model and realizable k-ε model could be used to simulate the flow and heat transfer of liquid hydrogen in a pipe.

The pipe will deform and vibrate when suffering from external transient impact. The pipe vibration attenuates rapidly, and the positive amplitude attenuates to 6.9% of the maximum amplitude 0.2 s after the transient impact, while the amplitude of pipe vibration decreases slowly with time when no fluid exists in the pipe.

After suffering from the external impact, the pressure of liquid hydrogen in the cross-section plane increases along the pipe deformation direction, and the pressure distribution is influenced for approximately 0.3 s though the external impact is transient. External transient impact enhances the disturbance of liquid hydrogen near the wall of the pipe. The increased flow resistance and the energy induced by the deformed pipe both affect the flow of liquid hydrogen, and contribute to the fluctuating characteristics of the liquid pressure drop. In addition, the phase state of liquid hydrogen remains unchanged in the pipe, i.e., the volume fraction of vapor phase in the pipe remains 0, indicating that little energy of the external transient impact is transformed into the internal energy of liquid hydrogen.