Thermal simulations of a hydrogen storage tank during fast filling

A finite element analysis is performed on the heat transfer process across the tank walls to determine the temperature distributions of hydrogen storage tanks during fast filling. The accuracy of the numerical model is shown by comparison between the experimental measurements and the computed results. A sensitivity analysis of the tank wall thermal conductivity, specific heat capacity, density and heat transfer coefficient between the tank's external surface and the ambient air is carried out and the resulting effects are described. The properties of the tank's composite layer have a larger effect on the temperature history on the tank external surface than the properties of the plastic liner. The heat transfer coefficient between the tank's external surface and the environment has a negligible effect during the filling but a significant impact during the holding time. Increasing the liner thickness significantly decreases the temperature in the composite


Introduction
In recent years a large investigation has been undertaken on the filling and emptying of hydrogen tanks at the Institute for Energy and Transport of the Joint Research Centre (JRC). The investigation has been carried out with an extensive program of experiments in the JRC GasTeF Facility [1,2] supported by computational fluid dynamics (CFD) simulations [3e6]. Other research groups have been performing experiments and CFD simulations in the same field [7e24]. The interest in the topic comes from the fact that compressed hydrogen in tanks is the technological solution currently adopted by car manufacturers to address the issue of hydrogen storage in vehicles.
Hydrogen is compressed to a very high pressure (70 MPa) to store sufficient energy in the tank to achieve a driving range comparable to that of conventionally powered vehicles. The short filling time and the high final pressure in the tank result in a temperature increase that can potentially be harmful for the mechanical properties of the tank materials. According to international Standards and Regulations, the maximum temperature within the vehicle fuel system must be equal or less than 85 C [25e28].
Many CFD studies on the fast filling of hydrogen tanks have been performed in the past decade [3e6, 8,13,15e17,19e22]. The focus of the large majority of the investigations hitherto has been the gas temperature while in this work we identify the temperature of the tank walls and the heat transfer across the walls as the main areas of interest. In previous numerical studies of the same group of authors the analysis had been performed with CFD techniques. In the present work, a Finite Element Analysis code has been chosen since it is a more suitable tool to investigate physical processes in solids like the heat transfer across multiple layers of material. Initially, the simulation results have been compared with experimental measurements to assess the accuracy of the numerical model. In a second step, the numerical model has been used to perform an analysis of the heat transfer process in the tank walls. Finally, a sensitivity analysis of the effect of the relevant material properties on the heat transfer has been carried out.
Currently, two types of composite tanks are being considered as on-board storage solution for hydrogen powered cars. In the type III tank, the inner layer (liner) is made of aluminum alloy, while in the type IV it is made out of plastic (e.g. high density polyethylene). In both types of tanks, the outer wall is made of carbon fibre reinforced polymer. The main function of the liner is the leak-tightness, while the outer composite layer provides the structural strength. Since the thermal conductivity of the plastic liner is smaller than of the metal one, the heat transfer from the gas to the tank is slower in the type IV than in the type III. Consequently, given the same initial and boundary conditions type IV tanks have higher gas temperatures than type III tanks. Because of this, the type IV tank is selected for the current investigation.

The experimental data
Experimental testing of a type IV tank has been performed at the Institute for Energy and Transport of the Joint Research Centre in the compressed hydrogen Gas tanks Testing Facility (GasTeF) [1,2]. In the current work experimental data from a test cycle composed of a filling (330 s), holding (16,150 s), and an emptying stage (365 s) is used. Fig. 1 shows the measured gas temperature close to Point1 located on the inner surface of the plastic liner. Point1 position is indicated in Fig. 2. During the filling stage hydrogen is injected into the cylinder, producing an increase of the tank gas temperature. Maximum temperature is obtained at the end of the filling process. During the holding stage no hydrogen is injected or removed from the tank. After the holding stage the hydrogen is rapidly removed from the tank, resulting in a rapid temperature decrease.

Model configuration
A finite element model is used to analyse the transient heat transfer problem. The models consist of three sections: two stainless steel bosses at the ends of the cylinder, the plastic liner in contact with the hydrogen, and a carbon fibre reinforced polymer (CFRP-composite) on the outside. 2D and 3D models have been generated as shown in Fig. 2. The 3D model is more representative of the reality of the problem but the axisymmetric 2D model requires a significant shorter computational time. Only a 90 section is modelled in a 3D model since both the geometry and the applied loads are symmetrical.
8-node, second order, axisymmetric finite elements are used for 2D axisymmetric model while 20-node, second order, finite elements are used for the 3D model [29]. A mesh density study has been performed to verify the mesh-independence of the computed results. Table 1 lists basic finite element model properties.

Thermal loads and boundary conditions
A temperature of 17.8 C, equal to the initial temperature of the hydrogen storage tank assembly, is used as initial  condition. In the experiment, the cylinder is placed in the horizontal position. For the condition on the inner surface two distinct approaches are explored: (1) application of convection boundary condition by using heat transfer coefficients (HTC), calculated with a CFD model [5,30] and (2) direct temperature specification.
In the first approach, the values of the HTC at the interface between the gas and the liner surface have been extracted from CFD calculations [5] where a good agreement between the experimental data and the simulation results for the gas temperature history was demonstrated. The heat transfer coefficient at the hydrogen-tank interface is a function of space and time and it is calculated by the CFD code when the conjugate heat transfer is included in the computation. In Fig. 3, the computed HTC history at three locations at the interface between the gas and the liner is shown. The computed HTC value (200e1800 W/(m 2 K)) has a similar range as found by Immel et al. [23] for a 129L type IV tank where the calculated maximum, minimum and average values are about 2800, 200e300 and about 1400 W/(m 2 K), respectively.
The computational time for the CFD simulation of the holding and emptying stages in 3D tank geometry is extremely long due to the number of computational nodes and long holding stage. It is therefore not feasible to apply the HTC approach beyond the filling stage.
In the second approach, the hydrogen temperature, measured in the gas just below Point1 of the plastic liner [2], Fig. 1, is applied as a time-dependent temperature boundary condition. The gas temperature is higher than the 85 C threshold between 181 and 462 s. During the filling stage of 330 s the hydrogen flow from the injector produces strong recirculation in the tank, creating a uniform gas temperature field in the cylinder as shown in the CFD calculations and in the experimental measurements [5]. The temperature boundary condition is therefore prescribed uniformly on the inner surface of the plastic liner and the steel sections of the bosses that are in contact with the hydrogen. At the end of the filling, the flow velocities decay rapidly to zero, and the effect of buoyancy becomes relatively more important. The gas in contact with the internal surfaces of the tank is cooled down while the warmer gas tends to accumulate on the top of the tank due to buoyancy. The gas temperature stratification starts to develop, resulting in the highest hydrogen temperature at the top part of the storage tank (Point1). Applying the measured temperature at Point1 as a boundary condition for the whole inner surface of the liner is a conservative approach for the whole vessel. On the outer surface of the cylinder natural convection heat transfer boundary conditions are applied with a wall HTC of 6 W/(m 2 K) [5] and a gas temperature of 17.8 C.

Material parameters
Thermal conductivity, density and specific heat capacity (SHC) of the plastic liner and the composite are taken from literature [31], while standard values for stainless steel, Table 2.  i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 2 5 6 0 e1 2 5 7 1

Results and discussion
Model validation

Filling stage
The comparison of the temperature history between the experiment and the simulations at Point3 (external surface) of the composite during the filling stage of 330 s is depicted in Fig. 4. Labels HTC refer to the first approach (convection boundary condition). Labels T refer to the second approach (temperature boundary condition). With the first approach, the computed temperatures are very close to the measured temperatures on the external surface for the first 150 s. Later, they diverge and the maximum temperature at the end of the process is underestimated by less than 2 C. With the second approach, the computed temperatures are initially slightly overestimated but at the end of the filling stage the values are closer to the measured ones compared to the first approach. The temperature boundary condition results in a better overall agreement of the computed external surface temperatures with the measured ones. Computational times of both approaches can be significantly decreased by increasing the time increment from 1 s to 10 s. This only slightly affects the computed temperatures, see Fig. 4 (top graph). The differences in results of the 2D and 3D calculations are practically negligible, Fig. 4, therefore only 2D models are used in the following sections. Labels T-1, T-2 and T-5 in Fig. 4 (bottom graph) represent sub-cases where the measured gas temperature has been decreased by 1, 2 and 5 C and then applied as a temperature boundary condition on the inner surface. The reason   for the sensitivity study on the temperature boundary condition is the uncertainty in the sensor positions. These are attached to the ribs of an umbrella-like structure, see Fig. 5. When the structure is inserted through a small opening into the tank, the ribs with the thermocouples are completely aligned with the tank central axis. Once the ribs are opened inside the vessel, it is not possible to verify the exact position of the sensors. The sensor was placed as close as possible to the tank wall, but it was not possible to determine the exact distance between the sensor and the wall. Because of the heat transfer and the temperature gradient in the gas in the region close to the wall, the actual inner surface temperature could be a few degrees lower than the measured temperature during the filling stage. This effect has been investigated numerically by decreasing the inner surface temperature. As shown in Fig. 4 (bottom graph), this effect has only a minor influence on the heat transfer during the filling stage.

The test cycle
Due to prohibitively high computational times for obtaining wall HTCs using CFD analysis, only the second approach is applied to the complete test cycle. Consistent with the previous comparison in the filling stage and as depicted in Fig. 4 (top graph), no significant differences can be observed in the results between 1 s and 10 s time increments and between the 2D and 3D cases in the whole cycle as shown in Fig. 6. Application of the unreduced internal surface temperatures results in a slight overestimation of the external surface temperatures. The maximum temperature is reached slightly later than in the experiment. Reducing the inner surface temperature by 1 C results in the best match with the measured values. Decreasing the inner surface temperature by 2 or 5 C causes an under-estimation of the maximum temperature, making those assumptions non-conservative. In order to be on the conservative side of the predictive model capabilities, the unreduced internal surface temperatures are applied to the following analyses.

Thermal analysis
In Fig. 7 the temperature and the heat flux distributions through the thickness are shown at selected times. At 0 s, the temperature and the heat flux profile are completely flat. As soon as the filling begins, the gas temperature and the heat flux increase. Consequently the temperature of the tank wall increases, initially in the liner and then in the composite. The gas temperature increases quickly in the initial stages of the filling whereby large temperature gradients between the internal and the external surface of the liner develop and generate large heat fluxes as shown by the blue curve in Fig. 7.
As the filling proceeds, the tank wall temperature increases, including on the external surface, as described by the red and green curves at 189.3 s and 329.3 s respectively. The temperature profile has a different slope at the interface between the plastic and the composite (Point2) due to the different thermal properties of the two materials.
The temperature at Point1 of the inner surface of the plastic liner reaches its highest value of 93.95 C at the end of the filling stage. The maximum temperature of the composite at this time is 65.68 C i.e. about 30 C smaller than the inner surface temperature. Hence the plastic liner acts as a thermal shield for the composite. After 330 s, the filling is stopped and the gas temperature starts to decrease but still remains above the threshold until 469 s. The temperature inside the composite layer continues to increase temporarily while the heat fluxes in the liner and partially in the composite decrease, as shown by the cyan curve at 469.3 s. Although at the end of the filling the gas temperature exceeds the 85 C threshold, the temperature in the composite layer remains well below the limit during the entire test cycle. Since the heat transfer is driven by the temperature gradient between adjacent layers, the heat transfer (flux) becomes slower as the temperature differences between layers decrease and the temperature and heat flux profile flatten as depicted at 7999.3 s. With a flat temperature profile, the heat transfer is very slow and it takes a very long time before the whole system returns to the initial temperature of 17.8 C. After more than 4 h, the temperature on the external surface is still by about 10 C higher than the initial 17.8 C.
At the end of the test cycle a fast emptying of the tank is initiated. The fast depressurization of the vessel generates a quick drop of the gas temperature and the heat transfer direction reverses. The temperature and heat flux profiles change completely as shown at 16,845 s.
The temperature distribution in the tank wall is shown in Fig. 8 at three different times: at the end of the filling when the maximum temperature (93.95 C) on the inner surface of Fig. 5 e Position of thermocouples [2].  the plastic liner is reached; at the time (1609 s) when the maximum temperature (53.84 C) of the tank external surface is attained; and at the end of the fast emptying (16,845 s). For a large portion of the tank cylindrical section, the temperature distribution is uniform. The thickness of the composite material is smaller in the shoulder region of higher curvature at the boundary between the dome and the cylindrical section. During the filling and holding stages the temperature of the shoulder region is therefore higher than in the straight sections. The situation reverses at the end of the emptying. At the two ends of the tank, the temperature distribution is quite different compared to the other parts of the tanks due to the larger thermal diffusivity of the metal bosses compared to the plastic liner and the composite material. The external temperatures at the ends are therefore higher than in the other parts of the tanks during the filling and holding time and lower at the end of the emptying. The small asymmetry in the temperature distribution between the two ends at 329 s is due to the constant temperature of the incoming gas.

Material sensitivity analysis
In literature different material property values can be found as shown in Table 3. To analyze the sensitivity of the applied material properties, conductivity, SHC and density are modified by multiplying values in Table 2 by 0.5, 0.8, 1.2, and 1.5. The external surface temperature at Point3 is computed for different values of the material properties. Different material properties of the tank affect the heat transfer and therefore the gas temperature history inside the vessel. In this investigation, the gas temperature history is fixed and the focus of the investigation is the heat transfer and the temperature gradient inside the tank wall materials and the final result on the temperature on the outer surface of the tank. Fig. 9 shows temperature histories at Point3 for different liner properties. As expected, increased thermal conductivity  i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 2 5 6 0 e1 2 5 7 1 produces a higher heat transfer and therefore a higher temperature on the external surface, causing a higher and earlier occurrence of the maximum temperature. After that time, the difference between the computed temperatures decrease. After about 4000 s there is practically no impact of the different conductivity since the temperature differences across the material are very small as shown in Fig. 7 for the time of 7999.3 s. The thermal conductivity has a larger effect than the SHC. During the filling stage, the temperature at Point3 only slightly Fig. 9 e Temperature history at Point3 for different plastic liner material properties. Left: filling stage. Right: full testing cycle. i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 2 5 6 0 e1 2 5 7 1 depends on the liner's SHC. Decreased SHC results in higher temperatures. Beyond the end of the filling stage the liner's SHC has a negligible effect. The density has the same effect on the temperature history as the liner's SHC. That is consistent with the expression of the thermal diffusivity (a) that measures the capability of a material to conduct thermal energy (k), relative to its capability to store thermal energy (r c p ), Eq. (1).

Liner
Since both the density and the SHC appear in the denominator, they have the same effect on the thermal diffusivity and therefore on the temperature history.
The maximum temperature in the composite layer occurs in Point2, see Table 4. Consistent with the temperature history in Point1, the thermal conductivity has a relevant impact on the maximum temperature in Point2 while density and SHC have a negligible effect. Decreasing the thermal conductivity lowers heat transfer in the liner and decreases the maximum temperature in Point2. The lowest maximum temperature in Table 4 is reached with the lowest thermal conductivity of the plastic layer.
Composite Fig. 10 shows the sensitivity of Point3 temperature history on the composite layer material properties. These have a larger effect than differences in the material properties of the plastic liner. This is mainly due to the significantly larger thickness and density of the composite layer which results in a significantly larger mass (and therefore thermal capacity) of the composite material compared to the plastic liner.
The increased composite's thermal conductivity accelerates the heat transfer across the material. The maximum temperatures on the external surface of the cylinder are higher and occur earlier in time. Similarly to the plastic liner, the effect of a different conductivity of the composite layer is larger when larger temperature gradients occur in the material. As soon as the temperature gradients decrease, the differences in the temperature histories reduce, Fig. 10.
Increasing the SHC or the density increase the capacity to store thermal energy and dampens the temperature increase.
As shown in Table 4, decreasing the density, the SHC, and the conductivity increases the maximum temperature in Point2. Whereas a decrease of the thermal conductivity in the composite layer increases the temperature of Point2, it decreases the temperature of Point3, Fig. 10.

External HTC
The effect of external wall HTC on the external tank temperature is illustrated in Fig. 11. The best overall match with the experimental data is obtained using the value of 7 W/(m 2 K). However, the maximum temperature at Point3 is slightly underestimated, which makes the usage of such value nonconservative. The simulation results for the value of 6 W/ (m 2 K) are also very close to the experimental data and conservative. Therefore the value of 6 W/(m 2 K) is the most appropriate.
Applying the value of 3, 9 or 12 W/(m 2 K) results in negligible differences in the temperature history during the filling stage, Fig. 11, and significant differences in the holding stage. The time is the critical factor to explain the different impact of the external heat transfer coefficient between the filling and holding stage. During the filling process the temperature difference between the external surface of the vessel and the environment is still too small (less than 10 C) to produce a large heat transfer between the tank and the environment and therefore the effect of a different HTC is limited. During the holding time, there is enough time for a larger temperature difference to develop between the cylinder surface and the environment and therefore the heat transfer with the environment becomes larger and the impact of a different HTC becomes higher.
As described in Table 4, the value of the HTC has negligible effect on the maximum temperature in Point2.

Liner thickness
In Fig. 12 temperature histories at Point2 at the linercomposite interface and at Point3 on the external tank surface are shown for different liner thicknesses. Increasing the liner thickness decreases the maximum temperature at Point2 from 70.9 C to 55.6 C. Since the temperature at Point2 is the maximum temperature reached in the composite layer, decreasing it by about 15 C is a relevant effect, especially since the threshold value of 85 C was set mainly due to the behaviour of the composite material.

Conclusions
Simulations of the filling and holding stages for a hydrogen type IV tank were performed with a finite element code. The model accuracy was initially validated by comparing the simulation results with the experimental measurements. A good agreement between the experiment and the simulation was shown for the temperature history at the middle section of the external surface of the vessel. No significant difference could be found in the simulation results between the 2D and 3D models, suggesting that 2D model suffices for the problem description. Different time increments, 1 s and 10 s, were applied with no relevant effects on the calculation results.
A sensitivity study was carried out on relevant material properties: thermal conductivity, specific heat capacity (SHC),   The thermal conductivity, the SHC, and the density of the composite layer have a relevant effect on the external temperature history during the filling and the first part of the holding time. Their effect is strongly decreasing in the second part of the holding time.
In general changes in the properties of the composite layer have a larger impact on the external surface temperature compared to changes in the properties of the plastic liner, due to the significantly smaller thickness and density of the plastic liner compared to the composite layer.
Different values of the external heat transfer coefficient have a negligible effect during the filling stage while they have a significant impact during the whole holding time.
Increasing the thickness of the plastic liner has a relevant impact on the maximum temperature in the composite layer, increasing the role of thermal shield that the liner plays for the composite layer.  i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 2 5 6 0 e1 2 5 7 1