Numerical Modeling for Rapid Charging of Hydrogen Gas Vessel in Fuel Cell Vehicle

: As a fuel for power generation, high-pressure hydrogen gas is widely used for transportation, and its efﬁcient storage promotes the development of fuel cell vehicles (FCVs). However, as the ﬁlling process takes such a short time, the maximum temperature in the storage tank usually undergoes a rapid increase, which has become a thorny problem and poses great technical challenges to the steady operation of hydrogen FCVs. For security reasons, SAE J2601/ISO 15869 regulates a maximum temperature limit of 85 ◦ C in the speciﬁcations for reﬁllable hydrogen tanks. In this paper, a two-dimensional axisymmetric and a three-dimensional numerical model for fast charging of Type III, 35 MPa, and 70 MPa hydrogen vehicle cylinders are proposed in order to effectively evaluate the temperature rise within vehicle tanks. A modiﬁed standard k - ε turbulence model is utilized to simulate hydrogen gas charging. The equation of state for hydrogen gas is adopted with the thermodynamic properties taken from the National Institute of Standards and Technology (NIST) database, taking into account the impact of hydrogen gas’ compressibility. To validate the numerical model, three groups of hydrogen rapid refueling experimental data are chosen. After a detailed comparison, it is found that the simulated results calculated by the developed numerical model are in good agreement with the experimental results, with average temperature differences at the end time of 2.56 K, 4.08 K, and 4.3 K. The present study provides a foundation for in-depth investigations on the structural mechanics analysis of hydrogen gas vessels during fast refueling and may supply some technical guidance on the design of charging experiments.


Introduction
The importance of the global climate and environmental protection challenges is expanding due to the increasing use of fossil fuels. Due to worsening environmental degradation and concerns about energy scarcity, it is vital to identify alternate energy sources [1].
Hydrogen is viewed as the most viable alternative energy form in the 21st century, since it is a clean, efficient, safe, and sustainable secondary energy [2]. Hydrogen is preferred over other fuels due to its high calorific value during burning, pure combustion products (H 2 O), and accessibility to sources [3]. The extremely low density of hydrogen gas (Standard Condition: 0.0899 kg/m 3 ), along with its storage and transportation challenges, prevents the widespread use of hydrogen energy [4].
The use of hydrogen as an automotive power source is anticipated to be crucial in lowering CO 2 emissions, which are now the biggest environmental problem [5]. According to the DOE (US Department of Energy) [6,7], the filling time with 5 kg of hydrogen should not exceed 200 s to be competitive with regular gas stations and appealing to potential customers. The biggest obstacle to the widespread deployment of FCVs is the noticeable temperature spike that occurs when refilling is sped up. In general, hydrogen is kept in a variety of forms, such as metal hydride, cryogenic liquid hydrogen, and compressed gas hydrogen [8]. Due to the established preparation technology, cheap cost, high energy efficiency, and practical hydrogenation, high-pressure gas hydrogen storage is the main utilization type in automobiles. To ensure the normal operation of FCVs, the American Society of Automotive Engineers SAE J2601 and International Standards ISO 15869 stipulate that the gas temperature should not exceed 85 • C during rapid filling [9][10][11][12]. This is mainly due to the limitations in the thermal properties of the materials currently available for the manufacturing of high-pressure hydrogen tanks.
Numerical modeling becomes a powerful tool to predict the fluid flow and thermodynamic behavior during hydrogen charging because the experimental test is relatively complex, has a long measure term, is expensive to construct, and the comprehensive analysis cannot be shown the temperature distribution inside the storage tank. Numerous researchers have carried out various numerical analyses on hydrogen rapid charging for FCVs. Dicken and Merida [13] developed a simplified two-dimensional axisymmetric model to predict the rise in gas temperature and pressure that occur during the filling of a hydrogen cylinder. Compared with an in situ measurement of the average temperature rise and temperature distribution, the results demonstrated that the proposed numerical model was able to predict the variation in temperature during refueling, and could be useful to determine the best locations for the onboard temperature sensor so that the local measurement best represents the mean gas temperature. Heitsh et al. [14] also constructed a three-dimensional numerical simulation using the ANSYS CFX software to depict the rapid filling process based on the experimental data recorded by Dicken and Merida [15]. Melideo et al. [16] developed a three-dimensional numerical model to study the hydrogen gas refueling. The effects of precooling on the temperature rise of the storage tank were comprehensively analyzed. Afterward, the effects of key parameters such as maximum temperature rise, state of charge, and energy cooling on the filling process [17] were numerically researched based on the proposed numerical model. With four different fuel-delivery temperatures, several filling rates, and initial pressure all taken into account, Miguel et al. [18] numerically studied the rapid filling scenarios in two different types of tanks. It was discovered that while the temperature rise and the final state of charge (SOC) decreased linearly with an increase in initial tank temperature, the final gas temperature increased linearly with both. Melideo et al. [19] quantitatively analyzed the hydrogen filling and emptying process to minimize the harm caused by the rise in temperature to the structural integrity of the hydrogen storage system. The findings revealed a 3 • C temperature differential at the conclusion of the refueling between measurements and calculations. The distribution of thermal stress was greatly impacted by the temperature stratification that was seen as a result of the prolonged defueling period. The heat transfer performance in a hydrogen gas storage-vessel under charging conditions was numerically investigated by Oh et al. [20]. They discovered that the end charging temperature lowers as the inflow temperature rises, and the final temperature rises as the average pressure ramp climbs. Comparative numerical investigations and analyses were carried out by Suryan et al. [12] to assess the efficacy of different turbulence models. The results demonstrated that when accuracy, convergence, and computational cost were taken into account, the realizable k-ε model was the best appropriate turbulence model for the hydrogen tank-filling issue. A three-dimensional numerical model was modified by Sapre et al. [21] to examine the effects of refueling factors on the storage density of compressed hydrogen in vehicle gas cylinders. The findings supported the notion that increasing the storage density of the tank is primarily influenced by pressure, filling time, supply temperature, and to a lesser extent by the temperature and filling rate. Meanwhile, Sapre et al. [22] investigated the stress and strain distribution in a mechanical performance study on hydrogen storage tanks. They investigated the structural integrity of storage tanks during filling and examined the location of the greatest stress and strain using the final parameters of various filling situations as the foundation for their mechanical study. In order to study the fast charging of a gas hydrogen storage vessel with a volume of 150 L and a pressure of 35 MPa, Zhao et al.
[8] suggested a two-dimensional axisymmetric numerical model. The proposed model and the measured experimental data had excellent agreement. Moreover, the thermodynamic performance in a hydrogen gas storage vessel during emptying was researched as well [23], based on the proposed numerical model. Wang et al. [24] numerically researched the fatigue test of a carbon fiber resin composite high-pressure hydrogen storage tank. By using hydrogen as the working medium, the fatigue property, failure behavior, and the influence of the safe hydrogen charging/discharging mode on the structure of the vehicle-mounted hydrogen storage tank were analyzed. Liu et al. [25] built a two-dimensional axisymmetric model to examine the fast charging and 10 min holding processes of a 150 L Type III and IV hydrogen storage vessel. The results displayed that for various tank types, the maximum temperature can emerge in various places (head/tail). The Type IV vessel's maximum temperature rise was higher than that of the Type III vessel. Thereafter, Liu et al. [26] numerically investigated the effect of the cylinder's inner and outer lining thickness on the temperature distribution.
As stated before, despite the fact that a lot of research has been completed on compressed hydrogen storage systems, most of it has focused on the rapid filling procedure, and the mechanical examination of the hydrogen vessel during filling or temperature rise control methods were rarely involved. The structural stability of the composite tank was obviously impacted by the high temperature and pressure created during refueling. As the execution of the hydrogen-related technologies depends on effective hydrogen storage techniques and mechanical analysis, an in-depth investigation of the modeling of hydrogen rapid charging should be conducted based on previous research. In this work, a thorough numerical model is developed to predict the interior temperature distribution and temperature rise of the storage tank while taking into account various heat-transfer modes. The thermodynamic performance during the fast refueling in FCV hydrogen storage tanks is predicted using the real gas equation of state, so that the compression effect of hydrogen gas can be taken into account. The thermal physical characteristics of the gas hydrogen come from NIST. The heat exchange between the tank wall and the ambient surroundings is considered as well. For the purpose of validating the developed numerical model, three series of filling experiment data are selected. After a detailed comparison, the numerical model is found to be suitable and can be used for in-depth investigations of hydrogen fast refueling. The current work may offer some technical guidance on the design of the filling experiment and the mechanical analysis of the tank structure, in addition to being important for the establishment of a compressed hydrogen refueling model.

Geometry
A two-dimensional axisymmetric and a three-dimensional numerical model were simultaneously taken into consideration and constructed for various experimental requirements to precisely forecast the variation in the internal temperature of the tank during the filling process. Three groups of filling experiments carried out by Dicken and Merida [13] and Liu [27] were adopted for model validation, and labeled as experiments No. 1, No. 2, and No. 3. The geometric model of the selected Type III hydrogen vessel is depicted in Figure 1, and the pertinent dimensions are listed in Table 1.

Materials
Due to the high-pressures during refueling of hydrogen storage tanks, it is necessary to take into account the compressibility impact of hydrogen gas. The thermodynamic properties of real hydrogen gas are taken from the NIST database using Fluent 2020R2. The thermal physical properties of solid materials are given in Table 2. Table 2. Thermal physical properties of solid materials [13,27].

Materials
Density

Boundary Conditions and Initial Settings
The cylinder inlet, inner wall, and outside wall boundary conditions are evaluated in detail during the model establishment. In order to match the recorded pressure at the cylinder's intake, the experiment's settings-based pressure inlet boundary conditions are adopted at the inlet of the gas vessel. A total temperature boundary condition is likewise applied at the inlet, and the flow direction is parallel to the inlet centerline. The conservation of mass, momentum, turbulent kinetic energy, and dissipation rate are solved at the inner wall of the cylinder using the non-slip boundary condition. In order to simplify the numerical model, the constant heat transfer coefficient is adopted to reflect the heat exchange between the outer wall of the cylinder and the environment, with values of 5 W/(m 2 ·K), 5 W/(m 2 ·K), and 10 W/(m 2 ·K) for experiments No. 1-3, respectively. The ambient temperature, which is 303 K, 285.3 K, and 293.4 K for experiments No. 1-3, respectively, is assumed to remain constant during refueling. The convergence residuals of continuity equation, momentum equation, k turbulent kinetic energy equation, and ε turbulent kinetic energy dissipation rate equation all have the same value of 10 −3 , while the convergence residual of energy conservation equation has a value of 10 −6 .
The initial temperature and pressure of the gas inside the cylinder serve as the initial conditions of the numerical model. The cylinder walls are assumed to have the same temperature as the gas. The gas pressure and temperature inside the cylinder are treated as being uniform at this point in the filling process. Figure 2 displays the detailed variations in the experimental parameters. It seems that while the gas temperature appears to fluctuate in diverse ways over time, the gas pressure appears to climb practically linearly. The reasonableness of the calculation results can be ensured by accurately arranging the circumstances for numerical simulation. The experiment results are fitted, and the fitted equations are inserted into the numerical model via the user-defined function (UDF) as boundary conditions, improving the accuracy of the simulated results. It is obvious that the fitted profiles closely match the experimental data with the maximum deviation of <10%. The boundary conditions and initial settings are therefore appropriate and acceptable.
temperature as the gas. The gas pressure and temperature inside the cylinder are treated as being uniform at this point in the filling process. Figure 2 displays the detailed variations in the experimental parameters. It seems that while the gas temperature appears to fluctuate in diverse ways over time, the gas pressure appears to climb practically linearly. The reasonableness of the calculation results can be ensured by accurately arranging the circumstances for numerical simulation. The experiment results are fitted, and the fitted equations are inserted into the numerical model via the user-defined function (UDF) as boundary conditions, improving the accuracy of the simulated results. It is obvious that the fitted profiles closely match the experimental data with the maximum deviation of <10%. The boundary conditions and initial settings are therefore appropriate and acceptable.

Grid Generation
The numerical calculation area includes three domains, i.e., the hydrogen gas region, aluminum liner, and CFRP layer. Both structured and unstructured meshes are used for computation. The solid region is generated with the structured grid, and the gas region is

Grid Generation
The numerical calculation area includes three domains, i.e., the hydrogen gas region, aluminum liner, and CFRP layer. Both structured and unstructured meshes are used for computation. The solid region is generated with the structured grid, and the gas region is generated with the unstructured grid. To precisely model the flow and heat transfer close to the wall, the boundary layer grids are incorporated. To lower the simulated inaccuracy, a high-resolution mesh is mapped in the entrance region where the largest gas velocity is formed. Figure 3 displays the details of the computational mesh. As the grid independence verification was performed in the original literature [13,27], it is not the focus of this study and thus not introduced here. to the wall, the boundary layer grids are incorporated. To lower the simulated inaccuracy, a high-resolution mesh is mapped in the entrance region where the largest gas velocity is formed. Figure 3 displays the details of the computational mesh. As the grid independence verification was performed in the original literature [13,27], it is not the focus of this study and thus not introduced here. The number of divided cells and nodes for experiments No. 1, 2, and 3 are 14,349/10,750, 147,258/637,973, and 17,885/12,622, respectively.

Governing equations
Once high-pressure hydrogen gas is charged into the storage tank, considerable heat exchange occurs between the injected gas and the tank wall, due to the high velocity gas flow. That is to say, obviously unsteady thermodynamic variation occurs in the hydrogen gas storage tank. Here, the Reynolds-averaged Navier-Stokes governing equations are utilized to simulate the hydrogen gas refueling process.
The mass conservation equation is given as follows: where ρ is gas density; t is the time; and v is the Favre average velocity.
The equation of the momentum conservation yields the following equation:

Governing Equations
Once high-pressure hydrogen gas is charged into the storage tank, considerable heat exchange occurs between the injected gas and the tank wall, due to the high velocity gas flow. That is to say, obviously unsteady thermodynamic variation occurs in the hydrogen gas storage tank. Here, the Reynolds-averaged Navier-Stokes governing equations are utilized to simulate the hydrogen gas refueling process.
The mass conservation equation is given as follows: ∂ρ ∂t where ρ is gas density; t is the time; and v is the Favre average velocity. The equation of the momentum conservation yields the following equation: where p is the pressure; = τ is the stress tensor; F y is the gravity source term; and u and I represent the dynamic viscosity and unit tensor, respectively. For the two-dimensional axisymmetric model, the influence of gravity can be neglected, so F y = 0.
The structure of the gas injection at the inlet primarily controls the flow field that forms in the storage tank. As a turbulence model, a modified standard k-ε turbulence model [13] 7 of 12 is used. This model is a semi-empirical model [28] that can resolve issues of a similar nature. The turbulence kinetic energy and turbulent dissipation rate are provided below: ∂ρk ∂t ∂ρε ∂t where G k and G b represent the turbulence kinetic energy generated by the average velocity gradient and the modified turbulence kinetic energy buoyancy term; β is the thermal diffusion coefficient of hydrogen; Y M represents the contribution of wave expansion to the overall energy dissipation; and M t is the turbulent Mach number, expressed as follows: where a represents the speed of sound and u t refers to the turbulent viscosity. In the above equations, C 1ε , C 2ε , C u , σ k and σ ε are model constants, with values of C 1ε = 1.52, C 2ε = 1.92, C u = 0.09, σ k = 1.0 and σ ε = 1.3.
In the fluid region, the energy equations are given as follows: k ge f f = k g + c p u t Pr t (14) where h is the enthalpy of gas; k geff is the effective thermal diffusion coefficient; and k g is the thermal diffusion coefficient.

The deviatoric stress tensor
= τ e f f can be denoted as: where the stress tensor u is replaced by u e f f = u + u t .

Heat Transfer on the Vessel Wall
In the solid region, the heat conduction through the tank wall and external natural convection equations are expressed as follows: where ρ w , T w , and k w represent the density, temperature and thermal conductivity of the tank wall, respectively. Additionally, h out is the external natural convection coefficient, and T amb is the ambient temperature.

Real Gas Model
Different from the Redlich-Kwong and modified Benedict-Webb-Rubin (MBWR) real gas equation of state adopted by Dicken and Merida [13] and Liu [27], the Helmholtz equation of state for normal hydrogen of Leachman et al. [29] supplied by the National Institute of Standards and Technology (NIST) is adopted and utilized to integrate the compressibility effects of high-pressure hydrogen.

Model Validation
A thorough study is conducted in the same environment as the boundary conditions to confirm the viability and accuracy of the produced numerical model. It is clear to see in experiments No. 1 and No. 2 the internal average temperature rises with the increase in charging time. Moreover, the simulated results have the same temperature rise trend and match well with the experimental data. Within the first few seconds, the average temperature rises quickly. Thereafter, the rate of ascent slows. In general, the numerical model can forecast the temperature distribution's variation trend within the storage tank. As Figure 4a shows, the numerical analysis errors reported by Liu [27] are within ±3%, while the inaccuracy of this model is within ±2%. According to this model, the final mean temperature differences are 2.56 K and 4.08 K, respectively. Figure 4b only displays a portion of the findings of the numerical analysis due to the high computational cost of the three-dimensional simulation, but it still provides useful reference information.
Both Figure 4 and 5 have been revised. Please check for details.    To increase the accuracy and improve the application scope of the proposed numerical model, the experimental data obtained by Dicken and Merida [13] are also chosen and adopted to conduct the model validation. Figure 5 illustrates the detailed parameter comparisons between the numerical simulation and experiment test. It is clear to see that the mean temperature, inflow mass flow rate, and internal pressure calculated by this proposed numerical model are all consistent with the experimental results and the other investigators' findings.
The comparison of the mean temperature of hydrogen gas between the numerical simulation and experimental data is shown in Figure 5a. It is easy to find that the simulated gas temperature and experimental data are in good agreement. The final discrepancy between the experimental results and the numerical simulation is only 4.3 K, with an error of <2%. The average anticipated temperature is, in the same way as other numerical models, consistently higher than the experimental data, suggesting that the numerical analysis may be more conservative. Qualitatively, the model accurately depicts the change in the cylinder's interior temperature, which rapidly increases at the start of filling, when the gas's rate of compression is at its highest. The interaction between the geometric model, inlet boundary conditions, and actual gas properties is the primary cause of the prediction deviations.      , and correlation between the mean temperature and inlet mass flow rate (d) [12,13,20].
The comparison of the internal pressures is displayed in Figure 5b. Since the curve fitting has large deviations in the experimental results, the inlet pressure boundary condition piecewise section experiences the largest inaccuracy. With the numerical and measured pressure of 12.32 MPa and 11.43 MPa, respectively, the highest inaccuracy of the numerical model is 7.24%. Figure 5c shows a comparison of the inlet mass flow rate between the numerical model and experimental results. It is evident that the inlet mass flow rate rapidly rises to the maximum at the beginning, then gradually decreases due to the increase in hydrogen density inside the storage tank. Generally, both the tank pressure and hydrogen mass within the tank increase with the hydrogen charging time. With the incremental increase in tank pressure, the difference between the filling pressure and tank pressure reduces, which makes it difficult to refuel the storage tank; thus the inlet mass flow rate decreases. The maximum mass flow rates are 0.072 kg/m 3 , 0.087 kg/m 3 , 0.090 kg/m 3 , and 0.092 kg/m 3 . As the filling comes to an end, the mass flow rate drops to 0.02 kg/m 3 .
The heat transfer between the fluid and the tank walls is largely influenced by the inlet mass flow rate; consequently, the gas temperature in the storage vessel also experiences considerable changes. The variations in the mean temperature of hydrogen gas within vessels under different inlet mass flow rates are shown in Figure 5d. The results demonstrate that the average temperature increases with the change of the whole inlet mass flow. During the refueling of the first 5 s, the temperature rises of different models in experiment No. 3 are 24.95 K, 33.89 K, 30.86 K, and 31.08 K. This tendency is comparable to Figure 5a, namely great temperature increases during the first few seconds of filling, followed by a gradual drop in the temperature rise rate. However, the relationship between the average temperature rise and inlet mass flow rate is affected by many other conditions, including the heat transfer of the tank wall, the velocity distribution inside the tank, and others. Therefore, in order to explain the current connection, further study is required to thoroughly investigate the relationship between the mass flow rate and the rise in fluid temperature in the storage tank. Figure 6 shows the temperature distribution within the compressed hydrogen cylinder. It is easy to find that experiment No. 1 lasts for 180 s, experiment No. 2 for 96 s, and experiment No. 3 for 8 s, 20 s, and 37 s. As Figure 6 shows, the maximum temperature predicted by the numerical model occurs near the end of filling and close to the inlet. This numerical finding is supported by experimental and other numerical temperature fields [13,14,25]. The gas temperature in storage tanks simulated by the two-dimensional axisymmetric model appears to be more uniform (shown in Figure 6a) when compared to that predicted by the three-dimensional model (shown in Figure 6b). This is mainly because gravity is taken into account in the three-dimensional numerical simulation. Due to the relatively lower temperature of the incoming hydrogen gas, apparent temperature changes are observed near the cylinder inlet.
Processes 2023, 11, x FOR PEER REVIEW 11 of 13 cause gravity is taken into account in the three-dimensional numerical simulation. Due to the relatively lower temperature of the incoming hydrogen gas, apparent temperature changes are observed near the cylinder inlet. Meanwhile, the plume flow appears on the centerline of the storage tank as a result of the action of the charging gas jet. As shown in Figure 6c, a backflow region occurs at the back of the cylinder when the incoming gas hits the end of the storage tank, which contributes to heat exchange enhancement between the fluid and inner wall. Additionally, with the last of filling times, the effect of heat transfer becomes more obvious, which Meanwhile, the plume flow appears on the centerline of the storage tank as a result of the action of the charging gas jet. As shown in Figure 6c, a backflow region occurs at the back of the cylinder when the incoming gas hits the end of the storage tank, which contributes to heat exchange enhancement between the fluid and inner wall. Additionally, with the last of filling times, the effect of heat transfer becomes more obvious, which is clearly shown in Figure 6a.

Conclusions
In this study, a numerical model of the refueling of a compressed gas cylinder in an FCV with both a two-dimensional axisymmetric and a three-dimensional model has been built and compared to the three groups of experimental data. The numerical model is able to predict the final mean temperature with differences in final average temperatures being 2.56 K, 4.08 K, and 4.3 K. The numerical results demonstrate excellent agreement with experimental data, and the calculation deviation is <2%.
To a certain extent, the proposed numerical model can wholly capture the tendency towards a temperature rise, and displays the internal temperature distribution of the storage cylinder; hence, where there is a rapid temperature growth within the first 5 s, followed by a steady decrease. This is different from the three-dimensional simulation, which considers the gravitational effect, the temperature in the bulk region of gas which surrounds the inlet region is more uniform in the two-dimensional axisymmetric simulation, and there coexists a large temperature gradient at the inlet of the cylinder. In addition, the backflow formed at the end of the cylinder greatly enhances the heat transfer between the tank wall and fluid. Consequently, a significant temperature gradient that develops and forms in the tank wall needs attention. As the obvious temperature distribution usually leads to serious thermal stress and strain, and further causes the failure of the structural integrity of the hydrogen gas vessel, the associated mechanical analysis of the gas vessel during refueling will be conducted in our future research.