Numerical Investigation on Influence of Gas and Turbulence Model for Type III Hydrogen Tank under Discharge Condition

Abstract

The high-pressure gaseous hydrogen (HPGH2) storage method is widely used owing to the low density of hydrogen gas at ambient temperature and atmospheric pressure. Therefore, rigorous safety analysis of the filling and discharging of compressed gas in a hydrogen tank is required to achieve reliable operational solutions for the safe storage of hydrogen. In this study, the behavior of compressed hydrogen gas in a hydrogen tank was investigated for its discharge. Numerical models for the adaptation of gas and turbulence models were examined. Gas model effects were examined to account for hydrogen gas behavior at the discharge temperature and pressure conditions. Turbulence model effects were analyzed to consider the accuracy of each model: the assessment of the turbulence models was compared in terms of the turbulence intensity. From the study of gas model effect, the Redlich–Kwong equation was found to be one of the realistic gas models of the discharging gas flow. Among the turbulence models, the shear stress transport model and Reynolds stress model predicted the compressed gas behavior more accurately, showing a lower turbulence intensity than those of the realizable and renormalization group models.

1. Introduction

Hydrogen is one of the best alternatives for the energy sources in the future, owing to its calorific value and environment-friendly properties [1,2], being a clean and sustainable energy source, solving air pollution and global warming problems. Because production of hydrogen is from carbon-free resources, it can affect greenhouse gas reduction [3]. The increasing energy demand and the depletion of available fossil energy sources have also necessitated the development of a new technique for the generation of hydrogen energy. Because hydrogen has a low density at ambient temperature and atmospheric pressure, hydrogen gas storage and transport technologies serve as technical obstacles, limiting the widespread use of hydrogen energy [4,5].

Hydrogen gas must be pressurized at ambient temperature to ensure the safety of gas storage and transportation [6]. Hydrogen storage methods can be categorized into pressurized gas, cryogenic liquid, and solid fuel storage [4]. Due to technical simplicity and high filling–releasing rate, high-pressure gaseous hydrogen (HPGH2) storage has become the well-known method [7].

Composite pressure vessels and tanks are suitably utilized to store compressed gaseous hydrogen because the low density of hydrogen requires it to be stored at a high pressure, particularly for adaptation in vehicular applications. The storage tank should be lightweight and of sufficiently high strength to withstand the applied high pressure [8,9,10,11,12]. Hydrogen tanks are composed of carbon fiber-reinforced plastic and metal or polymer liners to ensure a low weight, where the liner provides a gas-tight seal and the outer laminate provides the structural integrity necessary to the tank.

The prediction of compressed hydrogen gas (CHG) behavior is important because rigorous safety analysis of its filling and discharging in a hydrogen tank is required to achieve reliable operational solutions for the safe storage of hydrogen. Furthermore, accurate prediction of the gas state in a high-pressurized tank requires the knowledge of the process conditions of the tank in operation. The performance and characteristics of CHG behavior in a tank are determined under the conditions of pressure and temperature of the gas. These characteristics are typically determined experimentally. However, owing to experimental limitations, such as costs, and prohibitions due to safety reasons concerning high pressure, numerical calculations can be conducted to predict compressed gas behavior in a tank. Such numerical calculations should be based on the understanding of compressed gas behavior in a high-pressurized gas state under different conditions of pressure, temperature of the gas, and configuration characteristics of the tank. To address these issues, several researchers have conducted numerical analysis for predicting the gas behavior in a hydrogen tank. Xiao et al. [13] analytically obtained the hydrogen temperature in a tank using a simplified thermodynamic model of a compressed hydrogen storage tank as a function that fits experimental temperatures. Kuroki et al. [14] proposed a model to obtain the transient temperature rise of hydrogen when pre-cooled hydrogen is heated by the filling equipment at a refueling station. The abovementioned numerical analyses are important for controlling storage tanks; however, most numerical studies on compressed gas behavior in a pressurized tank have been conducted for filling and refueling cases.

A few studies on discharge cases have been performed by some researchers. Guo et al. [15] conducted a gas discharging study in which they conducted experiments and presented a model for simulating and describing the thermal behavior during a cycling test. Dhaou et al. [2] conducted an experimental and numerical study of a closed metal–hydrogen reactor, and proposed a theoretical model for describing the dynamic behavior of the reactor.

Most numerical studies on compressed gas behavior in a hydrogen tank have been conducted from the perspective of filling and refueling, with the estimation of the filling time and temperature increase, without considering the effects of turbulence and gas models. Hence, the validity of the numerical results of compressed gas behavior in hydrogen tanks remains uncertain. There are limited studies on numerical approaches for examining compressed gas behavior in a hydrogen tank during discharging. However, the need for such studies is increasing with the growing demand for various storage tank system operations and their control [1,2].

In this study, compressed gas behavior in a hydrogen tank during discharging was investigated. Numerical models were validated for the adaptation of gas and turbulence models. Gas model effects were investigated to account for the deviation with experimental results in the hydrogen gas behavior at the discharge temperature and pressure. Turbulence models were analyzed to assess the turbulence model accuracy in terms of the turbulence intensity.

2. Background and Theory

2.1. Background

Hydrogen storage technology can be categorized into solid, gas, and liquid storage systems based on the storage form of hydrogen. One of the common storage methods is the CHG storage method, in which hydrogen is stored under high pressure, e.g., 15–70 MPa. The low-density characteristics of hydrogen require a high-pressurized storage system; therefore, the storage tank should be lightweight and rigid to withstand high pressure. A CHG storage tank is composed of a liner in the inner part of the tank and a laminate on the tank body, where the liner and the laminate provide a gas-tight seal and structural integrity to the tank, respectively.

CHG storage tanks can be classified as types I–V tanks based on the material and arrangement, as listed in

Table 1. Type of hydrogen storage tank.

Type	Materials	Hydrogen Storage Pressure
Type I	All metal	17.5~20 MPa
Type II	Metal liner with hoop wrapping	26.3~30 MPa
Type III	Metal liner with full composite wrapping [16]	35 MPa~70 MPa
Type IV	Plastic liner with full composite wrapping [16]	70 MPa

. Type I tanks are composed of metals, e.g., steel and aluminum tanks. Type II tanks are metal tanks (aluminum) with filament windings, such as glass fiber/aramid or carbon fiber, around them. Aluminum/glass and steel/carbon or aramid tanks can withstand maximum pressures of 29.9 MPa, respectively. Type III tanks are made from composite materials, such as fiber glass/aramid or carbon fiber, with a metal liner (aluminum or steel) to withstand maximum pressures of up to 70 MPa. Type IV tanks are made from composite tanks, such as carbon fiber, with a polymer liner (thermoplastic). They can withstand an approximate maximum pressure of 70 MPa.

2.2. Theory

2.2.1. Real Gas Model

Various gas models have been adopted for gas flow analysis to accurately simulate gaseous nitrogen, natural gas, and hydrogen flows. For a hydrogen gas flow, four real gas models are typically implemented to describe its behavior [17]: Redlich–Kwong equation [18], Soave’s modified Redlich–Kwong equation [19], Aungier’s modified Redlich–Kwong equation [20], and Peng–Robinson equation [21].

For the Redlich–Kwong equation, the state relation of a gas is given as follows:

$$p=\frac{RT}{V_m-b}-\frac{a}{\sqrt{T{V}_m}\left(V_m+b\right)}$$

(1)

where a and b are constants that correct for the attractive potential of molecules and volume, respectively, and are expressed as

$$a=\frac{\left(0.42748R^2{T}_c{}^{2.5}\right)}{p_c},b=\frac{\left(0.08662R{T}_c\right)}{p_c}$$

(2)

where $V_m$, $p_c$, $T_c$ are molar volume, pressure at the critical point, and temperature at the critical point, respectively.

This equation was modified by Soave [19] by the Soave function, $\alpha (T,\omega )$, which is associated with temperature $(T)$ and the acentric factor $(\omega )$. The modified Redlich–Kwong equation is expressed as follows:

$$p=\frac{RT}{V_m-b}-\frac{a\alpha }{V_m\left(V_m+b\right)}$$

(3)

where a and b are expressed as

$$a=\frac{\left(0.42747R^2{T}_c{}^{2.5}\right)}{p_c},b=\frac{\left(0.08664R{T}_c\right)}{p_c}$$

(4)

and the Soave function, $\alpha$, and temperature are expressed as follows:

$$\alpha ={\left(1+\left(0.48508+1.55171\omega -0.17613{\omega }^2\right)\left(1-T_r^{0.5}\right)\right)}^2,T_r=T/T_c$$

(5)

where $\omega$ is the acentric factor for the species, and the Soave function, $\alpha$, is expressed for hydrogen as follows: $\alpha =1.202\exp \left(-0.30288T_r\right)$.

Aungier’s modified Redlich–Kwong equation is

$$p=\frac{RT}{V_m-b+c}-\frac{a}{V_m\left(V_m+b\right)T_r^n}$$

(6)

where a, b, and n are expressed as

$$\begin{array}{c}a=\frac{\left(0.42747R^2{T}_c{}^2\right)}{p_c},b=\frac{\left(0.08664R{T}_c\right)}{p_c},\\ n=0.4986+1.1735\omega +0.4754{\omega }^2\end{array}$$

(7)

The Peng–Robinson equation is expressed as follows:

$$p=\frac{RT}{V_m-b}-\frac{a\alpha }{V_m{}^2+2b{V}_m-b^2}$$

(8)

where a, b, and $\alpha$ are expressed as

$$\begin{array}{c}a=\frac{\left(0.45724R^2{T}_c{}^2\right)}{p_c},b=\frac{\left(0.07780R{T}_c\right)}{p_c},\\ \alpha ={\left(1+\left(0.37464+1.54226\omega -0.26992{\omega }^2\right)\left(1-T_r{}^{0.5}\right)\right)}^2,T_r=T/T_c\end{array}$$

(9)

2.2.2. Turbulence Models

(a) Realizable $k-\epsilon$ Model

Transport equations for turbulent kinetic energy ($k$) and dissipation rate ($\epsilon$) are expressed as

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

(10)

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

(11)

where $C_{3\epsilon }$ is turbulence constant for the Realizable $k-\epsilon$ model other constants are expressed: $C_1=\max \left[0.43,\frac{\eta }{\eta +5}\right]$, $\eta =S\frac{k}{\epsilon }$, and S means strain rate magnitude expressed by $S=\sqrt{S_{ij}{S}_{ij}}$.

$G_k$ and $G_b$ represent the turbulence kinetic energy generated by the mean velocity gradients and buoyancy, respectively, and $Y_M$ is the contribution of the fluctuating dilatation in the compressible turbulence to the overall dissipation rate. ${\sigma }_k$ and ${\sigma }_{\epsilon }$ are the turbulent Prandtl numbers for $k$ and $\epsilon$, respectively. $C_2$ and $C_{1\epsilon }$ are turbulence constants.

(b) Renormalization Group (RNG) $k-\epsilon$ Model

The RNG $k-\epsilon$ model is the modified $k-\epsilon$ model is based on the original Lauder and Spalding model [22]. The transport equations are summarized below.

$$\rho \frac{Dk}{Dt}=P-\rho \epsilon +\frac{\partial }{\partial x_i}\left(\mu {\alpha }_k\frac{\partial k}{\partial x_i}\right)$$

(12)

$$\rho \frac{D\epsilon }{Dt}=\frac{k}{\epsilon }\left(C_{\epsilon 1}P-C_{\epsilon 2}\rho \epsilon \right)-\rho R+C_{\epsilon 3}\rho \epsilon \nabla ·U+\frac{\partial }{\partial x_i}\left(\mu {\alpha }_{\epsilon }\frac{\partial \epsilon }{\partial x_i}\right)$$

(13)

where the production, P, and R are defined as follows:

$$P=2C_{\mu }\rho \frac{k^2}{\epsilon }\left[S_{ij}{S}_{ij}-\frac{1}{3}{\left(\nabla ·U\right)}^2\right]-\frac{2}{3}\rho k\nabla ·U,R=\frac{C_{\mu }{\eta }^3\left(1-\eta /{\eta }_0\right){\epsilon }^2}{1+\beta {\eta }^3}·\frac{{\epsilon }^2}{k}$$

(14)

The constants used in the RNG model are summarized below.

$$\begin{array}{c}{C}_{\epsilon 1}=1.42,C_{\epsilon 2}=1.68,C_{\epsilon 3}=\frac{-1+2C_{\epsilon 1}-3m_1\left(n_1-1\right)+{\left(-1\right)}^{\delta }\sqrt{6}·C_{\mu }{C}_{\mu }\eta }{3},\\ C_{\eta }=\frac{\eta \left(1-\eta /{\eta }_0\right)}{1+\beta {\eta }^3},C_{\mu }=0.0837\end{array}$$

(15)

The parameters used in the constants are as follows:

$\beta =0.012$, ${\eta }_0=4.38$, $\eta =\frac{Sk}{\epsilon }$, $S={\left(2S_{ij}{S}_{ij}\right)}^{1/2}$, $\delta =1$ when $\nabla ·U<0$ and $\delta =0$ when $\nabla ·U>0$, and ${\alpha }_k={\alpha }_{\epsilon }=1.39$.

(c) Shear Stress Transport (SST) Model

The SST model is a combination of the $k-\epsilon$ model of Launder and Sharma and the $k-\omega$ model of Wilcox. The eddy viscosity is defined as a function the turbulence kinetic energy and the specific dissipation rate as follows:

$${\mu }_t=\rho \frac{k/\omega }{\max \left(1;\mathsf{\Omega }{F}_2/a_l\omega \right)}$$

(16)

The maximum value of the eddy viscosity is bounded with the turbulent shear stress by the specific dissipation rate time, $a_l$. This effect is obtained using an auxiliary function, $F_2$, and the absolute value of the vorticity, $\mathsf{\Omega }$. The auxiliary function is defined as a function of the wall distance, $l$, as follows:

$$F_2=\tan \left\{{\left(\max \left[2\frac{\sqrt{k}}{0.09\omega l};\frac{500\mu }{\rho l^2\omega }\right]\right)}^2\right\}$$

(17)

In this model, there are two transport equations using a blending function, $F_1$. The transport equation for the turbulent kinetic energy and specific dissipation rate is expressed as follows:

$$\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_{j}k-\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right)={\tau }_{ij}{S}_{ij}-{\beta }^{\prime }\rho k\omega$$

(18)

$$\frac{\partial \left(\rho \omega \right)}{\partial t}+\frac{\partial }{\partial x_j}\left(\rho u_j\omega -\left(\mu +\frac{{\mu }_t}{{\sigma }_{\omega }}\right)\frac{\partial \omega }{\partial x_j}\right)=P_{\omega }-\beta \rho {\omega }^2+2(1-F_1)\frac{\rho {\omega }_{\omega 2}}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}$$

(19)

where $P_{\omega }$ is the production term of $\omega$, which is expressed using vorticity as follows:

$$P_{\omega }\equiv 2\gamma \rho \left(S_{ij}-\frac{\omega S_{nn}{\delta }_{ij}}{3}\right)S_{ij}=\gamma \rho {\mathsf{\Omega }}^2$$

(20)

The blending function, $F_1$, is defined as

$$F_1=\tanh \left\{{\left(\min \left[\max \left[\frac{\sqrt{k}}{0.09\omega l};\frac{500\mu }{\rho l^2\omega }\right];\frac{4\rho {\sigma }_{\omega 2}k}{C{D}_{k\omega }{l}^2}\right]\right)}^4\right\}$$

(21)

where $C{D}_{k\omega }$ is a cross-diffusion term as follows:

$$C{D}_{k\omega }=\max \left[\frac{2\rho {\sigma }_{\omega 2}}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j};{10}^{-20}\right]$$

(22)

The turbulence constants in the model are jointly denoted by the symbol, $\varphi$, defined by blending the coefficients in the $k-\omega$ and $k-\epsilon$ models, and denoted as ${\varphi }_1$ and ${\varphi }_2$, respectively [23], as follows:

$$\varphi =F_1{\varphi }_1+\left(1-F_1\right){\varphi }_2,\varphi =\left\{{\sigma }_k,{\sigma }_{\omega },\beta ,\gamma \right\}$$

(23)

(d) Reynolds Stress Model (RSM)

The RSM, a second-order closure model, can be obtained by solving the following transport equations [24]:

$$\frac{\partial }{\partial t}\left(\rho \langle u_i{}^{\prime }{u}_j{}^{\prime }\rangle \right)+C_{ij}=D_{\tau ,ij}+\frac{\partial }{\partial x_k}\left(\mu \frac{\partial }{\partial x_k}\langle u_i{}^{\prime }{u}_j{}^{\prime }\rangle \right)+P_{ij}+{\mathsf{\Phi }}_{ij}-{\epsilon }_{ij}$$

(24)

where

$$C_{ij}=\frac{\partial }{\partial x_k}\left(\rho \langle u_k\rangle \langle u_i{}^{\prime }{u}_j{}^{\prime }\rangle \right)$$

(25)

$$P_{ij}=-\rho \left(\langle u_i{}^{\prime }{u}_k{}^{\prime }\rangle \frac{\partial \langle u_j\rangle }{\partial x_k}+\langle u_i{}^{\prime }{u}_k{}^{\prime }\rangle \frac{\partial \langle u_i\rangle }{\partial x_k}\right)$$

(26)

and

$$D_{\tau ,ij}=-\frac{\partial }{\partial x_k}\left(\frac{{\mu }_t}{{\sigma }_k}\frac{\partial \langle u_i{}^{\prime }{u}_j{}^{\prime }\rangle }{\partial x_k}\right)$$

(27)

The pressure–strain term, ${\mathsf{\Phi }}_{ij}$, treats the transport and production processes of the Reynolds stress. Because this term mainly plays an important role in precisely predicting anisotropic flow characteristics, it determines the structure of turbulent flows.

3. Numerical Implementation

3.1. Governing Equations

For the governing equations of a gas fluid with a turbulent flow, the conservation of mass and momentum is considered with Reynolds density-averaged variables in the Cartesian tensor notation [25]. The conservation equations of mass and momentum are expressed as follows:

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_i}\left(\rho u_i\right)=0$$

(28)

$$\frac{\partial }{\partial t}\left(\rho u_i\right)+\frac{\partial }{\partial x_j}\left(\rho u_i{u}_i\right)=-\frac{\partial p}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\mu \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}-\frac{2}{3}{\delta }_{ij}\frac{\partial u_l}{\partial x_l}\right)\right]+\frac{\partial \left(-\rho \overline{u_i{}^{\prime }{u}_j{}^{\prime }}\right)}{\partial x_j}$$

(29)

where $\rho$, p, t, and ${\delta }_{ij}$ are the density, pressure, time, and Kronecker delta tensor, respectively. The velocity components denoted by $u_i$ (i = 1, 2, 3), where subscript i = 1, 2, 3 denotes the x, y, and z directions (subscripts j, l have the same definitions), respectively. Superscripts ${}^{\prime }$ and $\overline{}$ denote the turbulent fluctuating component and the Reynolds time-averaged component, respectively.

3.2. Calculation Assumptions and Conditions

Figure 1 shows the cross-section of a 74-L Type III hydrogen tank. The diameter and dimensions of the tank are provided in

Table 2. The diameter and dimensions of Type III hydrogen tank.

Item		Length (mm)
Cylinder type	Total length (mm)	1030
	Inner diameter (mm)	354
	Outer diameter (mm)	427
	Thickness (mm)	27–30

. The tank cylinder is consisted of a combination of a liner and a fiber laminate. Therefore, a fluid subdomain (tank interior filled with hydrogen) and different solid subdomains (internal liner, external composite laminate, and bosses at the tank end) were considered as the computational domain. The material properties and discharge conditions used in numerical calculation were considered in the works of Guo et al. [15] and Suryan et al. [26,27]

The numerical calculations are based on the following assumptions:

(1)

Natural convection coefficient for the heat exchange of the laminate is constant, and the laminate has isotropic heat transfer characteristics.

(2)

The outlet flowing temperature of hydrogen gas is a constant.

(3)

The thicknesses of the liner and laminate are constants.

(4)

The temperature of the entire tank is considered as ambient temperature because the tank is exposed to the ambient temperature for sufficient time.

Unsteady, compressible Naiver–Stokes equations were solved. Numerical calculations were performed for hydrogen discharge case at a constant flow rate. The working fluid was a hydrogen gas flow, which was simulated using different real gas models. A non-slip boundary condition was applied at every fluid–solid interface. For the treatment of mesh around the bosses, an inflated layer of hexahedral elements was generated. In overall, a structured and non-adaptive grid (hexahedral grid) was generated with 200,000 nodes. To achieve smoothness, convergence base of the residual for the velocity and for the mass, energy balance was considered below 10⁻⁴ and 10⁻⁶, respectively. A high-resolution advection scheme that was two orders of magnitude below the default convergence level of the solver was used [28,29]. Governing equations with the boundary conditions were solved using the finite volume method, the computation was carried out by using the commercial computational fluid dynamics (CFD) package, Fluent (ANSYS Fluent, version 18.0) [30].

4. Results

4.1. Gas Model Effect

Validation studies were performed to verify the adoptability of the computational methodology and numerical formulation by comparing the numerical results to the experimental results. Suryan et al. [26,27] conducted validation studies for the filling of a hydrogen tank. In that study, it was found that the Redliche–Kwong model provides the smallest deviation from the experimental results, among the real gas models. Moreover, Suryan et al. [26,27] examined the adoptability of turbulence models for the filling of a hydrogen tank. For the filling of hydrogen in the tank, it was observed that the realizable $k-\epsilon$ model and the RSM are suitable turbulence models considering the accuracy, convergence, and computational cost.

For hydrogen discharge case, the effects of the different gas models were examined by comparing the numerical results with the experimental results of Guo et al. [15] to determine the suitable real gas and turbulence models.

The experiment for the hydrogen discharge case was performed with gas cycling test which is consisted of fast refueling and slow discharging process. After fast refueling of hydrogen gas, the discharging experiment was conducted when the temperature of the hydrogen gas and the cylinder itself became generally uniform [15]. The detailed procedures are as follows: at the discharging process, with opening of the outlet valve, the hydrogen gas was released to atmosphere. During the discharge process, the gas pressure and temperature in the cylinder were measured by a pressure transducer and thermocouples, respectively. For the experiment of hydrogen discharge case, 331 K of initial gas temperature inside the cylinder was observed, while the ambient temperatures were 305 K. The gas pressure changed from 36.3 MPa to 2.5 MPa after discharge of hydrogen gas during 1140 s [15].

Regarding adoptability of the gas models, the following four gas models were considered: Redlich–Kwong equation [18], Soave’s modified Redlich–Kwong equation [19], Aungier’s modified Redlich–Kwong equation [20], and Peng–Robinson Equation [21].

The calculations for the hydrogen tank were conducted with the inclusion of compressibility effects. Therefore, the behavior of hydrogen gas was different from those obtained by the real gas models. Figure 2 and Figure 3 show the results of the averaged temperatures and pressure of the inner gas tank for the different gas models and experimental results, respectively. The inner temperature decreased with the discharge time. During the early stage of discharge at 300 s, the results of the real gas models deviated by 0.5% relative to the experimental results. This is due to the occurrence of abrupt pressure and temperature changes with the increase in the range of the discharging velocity of the gas in the compressible gas stage. After 400 s, the deviation of the results was reduced to 0.15–0.4%. In this state, the thermal mass of gas in the tank decreased with the reduction in the gas discharging velocity. The pressure and temperature changes also decreased; therefore, the gas state was more stable than that during the early stage.

For the gas models of Soave’s modified Redlich–Kwong, Aungier’s modified Redlich–Kwong, Peng–Robinson, and Redlich–Kwong equations, the averaged deviations from the experimental results were 0.17, 0.39, 0.23, and 0.15%, respectively. Among them, the Redlich–Kwong equation showed the smallest deviation from the experimental results; therefore, the gas model for the Redlich–Kwong equation was regarded as one of the realistic models for the gas flow. For the Redlich–Kwong equation, gas behavior was described effectively, reflecting gas phase properties to implement gas flow behavior: similar results were reported by Suryan et al. [26,27]. However, Soave’s modification introduces a Soave function, involving temperature and the acentric factor. Acentric factor is a conceptual number meaning of measure of non-sphericity of molecules. For the Aungier’s modified equation, it also introduces molar volume effect. These modified equations are effective on modeling especially for the real gas effects when hydrocarbons are involved [26,27]. Therefore, the gas model for the Redlich–Kwong equation was effectively described for the hydrogen gas flow.

4.2. Turbulence Model Effect

The turbulence models—realizable $k-\epsilon$, RNG, SST, and RSM models—were assessed for the CHG average temperature. The average gas temperature results of the different turbulence models were compared to the experimental results. Figure 4 shows the comparison of the average gas temperature during discharging obtained by the different turbulence models. The results showed that the inner temperature decreased with the discharge time. The maximum deviation from the experimental results was 0.5% at 300 s. The deviation from the experimental results decreased with the discharge time. The SST model showed the closest agreement to the experimental results, among the turbulence models, with 0.23% deviation. The realizable $k-\epsilon$ model predicted poorly among the turbulence models.

Between the two models, the realizable $k-\epsilon$ and RNG models, RNG models used an isotropic eddy viscosity formulation. Turbulent kinetic energy production terms was summed from the formulation, indicating more accurate predictions for RNG model than the realizable k − ε model [25].

Figure 5, Figure 6, Figure 7 and Figure 8 shows the turbulence intensities along the tank length obtained by the different turbulence models at 150 s, 250 s, 350 s, and 500 s. Typically, the intensity increased near the discharging area, where there is a large pressure decrease and velocity increase. The realizable $k-\epsilon$ model predicted a higher turbulence intensity, particularly near the discharging area, with a higher velocity increase. The RNG model also predicted a higher turbulence intensity; however, the realizable $k-\epsilon$ model yielded a higher turbulence intensity than the RNG model. This was because the excess turbulence kinetic energy was overestimated, particularly for the discharging area. In the RNG model, the turbulence energy dissipation equation was revised to modify the differential stress model and deal with flows with a high rate of deformation [25]. In addition, the rates of turbulence energy production and energy dissipation was influenced by modifying constant of $C_{\epsilon 2}$. Especially in the discharging area, reducing the value of $C_{\epsilon 2}$ eventually makes possible an increase the energy dissipation rate, reducing the turbulent kinetic energy [25]. Therefore, the turbulence intensity obtained by the RNG model was slightly lower than that by the realizable $k-\epsilon$ model, indicating that the RNG mode provides more accurate predictions than the realizable k − ε model.

In contrast, the SST model and RSM predicted lower turbulence intensities than the realizable $k-\epsilon$ and RNG models. The amount of turbulent kinetic energy is related to the amount of turbulent kinetic intensity [25]. Higher turbulent kinetic intensity was shown by the higher turbulent kinetic energy in the realizable $k-\epsilon$ model. On the other hand, lower kinetic energy was predicted for the SST model and RSM especially at the discharging area.

5. Concluding Remarks

In this study, compressed gas behavior in a hydrogen tank was investigated during discharging. Numerical models were investigated for the adaptation of gas and turbulence models. Gas model effects were examined to account for the hydrogen gas behavior at the discharge temperature and pressure. Turbulence model effects were analyzed to compare the accuracy of each model relative to experimental results, and the turbulence models were assessed in terms of the turbulence intensity. The main conclusions are as follows:

• Regarding the gas model adoptability, four gas models were considered, and the numerical and experimental results were compared for the discharge of CHG. Among the gas models, Redlich–Kwong equation was regarded as one of the realistic models for gas flow.
• The realizable $k-\epsilon$, RNG $k-\epsilon$, SST, and RSM turbulence models were assessed for the CHG average temperature. The SST model results showed the closest agreement with the experimental results among the results of all the turbulence models, whereas the realizable $k-\epsilon$ model yielded poor predictions.
• The turbulence intensity obtained by the different turbulence models increased near the discharging area, where there was a large pressure decrease and velocity increase.
• The turbulence intensity yielded by the RNG model was slightly lower than that by the realizable $k-\epsilon$ model, indicating that the RNG model provides more accurate predictions than the realizable $k-\epsilon$ model.
• For the SST model and the RSM, these models predicted a lower turbulence intensity than the realizable $k-\epsilon$ and RNG models.
• The present findings regarding suitable gas and turbulence models for predicting CHG flow behavior can be effectively applied in research area of hydrogen applications. The results will help in identifying a suitable numerical model for formulating an optimum tank filling strategy that is both efficient and safe.