Heat transfer analysis of high pressure hydrogen tank fillings

An analytical model is developed to predict the tank shell temperature distribution. (cid:1) The analytical model is validated against experimental results. (cid:1) The inﬂuence of ﬁlling parameters on the tank wall peak temperature is assessed. (cid:1) The tank materials heat diffusivity is the critical parameter for fast ﬁllings.


Introduction
In the light of the EU objectives to reduce green house gases emissions, the focus on alternative fuels has gained interest over the recent years [1].In particular, highly pressurized hydrogen gas is viewed as a relevant candidate to match actual fossil fuels-powered vehicles ranges [2,3].Several authors have in addition identified mixed hydrogen-compressed natural gas (CNG) as potential fuel for internal combustion engines to reduce CO 2 and all related pollutants [4e7].
A major topic in the study of gas filling processes is the prediction of temperature evolution in tanks due to high compression.In order to ensure integrity of on-board storage tanks during filling, the current standards recommend a maximum temperature of 85 C.This led to the development of filling protocols of high pressure vessels that involve precooling of the hydrogen [8e10].
These filling protocols yet need to be extended to heavyduty gas tanks.Experimental investigations of hydrogen tank fillings are expensive and cannot be conducted to investigate a broad range of impacting parameters.While CFD simulations of such tank filling scenarios bring insights on heat transfer [11e14], the flow complexity renders them very expensive.
It has been assessed by several authors [15e23] that simpler thermodynamical models can accurately predict temperature evolutions of such systems despite crude assumptions.The basis of such models is the assumption that due to sufficient turbulent mixing, pressure, temperature and thermal properties of hydrogen vary only slightly within the tank during filling.The simplest models originally assumed that the tank shell temperature is homogeneous, which leads to a 0-d system [15,16].Such models combined energy balance for both hydrogen and the solid tank shell.However, as noted by Bourgeois et al. [18], heat diffusion through the tank shell plays a particularly important role, and it is therefore inaccurate to consider homogeneous temperature within the tank shell.
The main challenge of thermodynamic models is to correctly account for heat transfer from the hydrogen to the tank inner wall.Several authors have proposed improvements; some use constant heat transfer coefficient [15,24e27] and others Nusselt correlations [16,28].By linking the heat transfer coefficient to the strength of the incoming turbulent hydrogen jet, the latter allow for more precise descriptions of heat transfer and hydrogen temperature evolution.The different approaches are extensively described in the review by Bourgeois et al. [29].
In the present paper, first the physical model and the main assumptions are introduced, then a robust numerical scheme to couple the 0-d hydrogen description and the energy equation in cylindrical coordinates in the solid is introduced.The obtained conjugate heat transfer solver is validated against experiments.The aim was to assess the impact of filling parameters on the peak temperature within the tank shell, with a special focus on the effects of tank materials' thermal diffusivity.

Physical modeling
A conjugate heat transfer solver was developed based on thermodynamic relations and heat diffusion equation.This section presents the equations and assumptions used by the model.

Gas modeling
Homogeneous hydrogen temperature and pressure distributions are assumed within the tank.This assumption is well verified during the filling process, where mixing due to turbulence negates the buoyant stratification.The temperature evolution of the hydrogen is described using the 1st law of thermodynamics as applied to an instationary open system, combined with the fundamental thermodynamic equation The equation for the hydrogen temperature evolution then reads Hydrogen temperature and mass are also related by the real gas equation of state as pV ¼ m g RZT g : (4)

Heat transfer modeling
The heat flux _ Q g/w from the gas to the tank inner wall is written as according to Newton's law of cooling, where the heat transfer coefficient k g/w between gas and tank inner wall is determined using experimental Nusselt correlations.As both forced and natural convection contribute to the overall heat transfer during filling, the effective Nusselt number is computed as following the approach of [30], where the parameter n usually takes a value between 1 and 4. The forced convection Nusselt number reads Nu for ¼ a*Re 0:67 inl ; as initially suggested by [31], where the Reynolds number Re inl is calculated for inlet nozzle conditions as In the present study the value of the parameter a is determined by fitting experimental data.It is common to use the Prandtl number in such forced convection correlations.For hydrogen the value of the latter is almost constant over the whole range of pressure and temperature of interest, close to 0.7 [16], and therefore is embodied in the constant parameter a.The natural convection Nusselt number is computed according to the empirical correlation Nu nat ¼ 0:104*Ra 0:352 g ; taken from [32], where the Rayleigh number Ra g is calculated for hydrogen gas conditions as Following [31], the parameter n is set to 1, such that the Nusselt correlation finally reads The heat flux from the tank outer wall to the ambient air is used as a boundary condition for the cylindrical heat diffusion solver.A tank on-board a vehicle is typically enclosed, and thus only subject to natural convection during filling.Accordingly, the heat transfer coefficient is set to a constant value of k w/a ¼ 8 [W/m 2 /K] [33].

Tank shell modeling
Although the assumption of homogeneous temperature holds for the gas phase during filling, is it not valid for the tank shell.The Biot number is a formal indicator of the non-uniformity of the temperature inside the solid phase.While it never is explicitly calculated in [18,24], the values of Bi g/w based on their data ranges from 100 ~500.These high values show that it is important to account for thermal diffusion through the tank shell.
Hydrogen Type IV tanks are composed of two layers, i.e., a thin thermoplastic liner on the inside, wrapped by a composite carbon fiber shell (CFRP).The tank is approximated by a cylinder with the same internal exchange area as well as the same internal and external diameters as the real tank.The temperature profile across the tank shell is obtained by solving the cylindrical heat diffusion equation Typically, two aluminium bosses are also enclosed at the top and bottom of the tank.In the present work, the results from [18], where it is estimated that 25% of the heat generated during filling is transported to the atmosphere through those bosses, are imposed rather than modeling the latter.

Conjugate heat transfer solver
The conjugate heat transfer problem is composed of two solvers, one for the hydrogen (fluid) and the other for the tank shell (solid).The fluid solver computes the hydrogen temperature and mass by solving Eqs. ( 3) and ( 4), and the solid solver computes the temperature profile through the tank shell by solving Eq. ( 14).As the inner wall temperature T w,i is needed to compute the heat flux _ Q g/w in the fluid solver, and the latter is needed as boundary condition for the solid solver, coupling is required.For this purpose, a Flux Forward Temperature Back (FFTB) algorithm is implemented in Matlab.The flow chart of the FFTB algorithm [34] is shown in Fig. 1.The convergence criterion is based on evaluating the temperature residuals of the hydrogen and of the first cell of the tank shell.
Usually, no more than 3 iterations are required to lower the residuals below the tolerance, chosen here as 10 À6 .The energy equation is solved in cylindrical coordinates using the Finite Volume approach.According to Eq. ( 14), the semi-discrete energy balance of Finite Volume cell U reads The temperature gradients are computed at the inner and outer interfaces of U.
For the results to be grid-independent, and in order to keep the ratio of spatial stepping from one cell to the next close to 1 across the tank shell, the liner and CFRP wrapping are discretized by 12 and 50 cells respectively.The gradients are discretized using a 2nd order scheme, and implicit Euler time integration is employed.
Thermodynamic and transport properties of hydrogen, namely c p , b, h, Z, m and k, are taken from the NIST database [35]; as functions of temperature and pressure.

Validation experiments
Hydrogen filling experiments of a type IV tank were carried out at the Hydrogen Refueling Station (HRS) of Empa (Du ¨bendorf, Switzerland).Table 2 presents the characteristics of the experimental tank.
Eight type K thermocouples were mounted on a specifically designed "tree", as depicted in Fig. 2, to measure the hydrogen temperature in the tank.Another thermocouple and a pressure sensor are placed near the inlet.It was observed that if the inlet velocity is lower than 5 m/s, temperature stratification occurs towards the end of the filling phase.This is consistent with the observations made by [18,36].An injector with a nozzle diameter of 3 mm was chosen to ensure high inlet velocities, such that turbulent mixing nullify the buoyant stratification.
Data from four hydrogen tank filling experiments were used to validate the proposed model.Thereby, varying initial pressure and ambient temperature values were considered.Table 1 describes the different parameters for each experiment.The Average Pressure Rampe Rate (APRR) was determined based on the CEP-A70 fueling protocol (ver.2010), and the tank was filled up to 70 MPa in every experiment.According to this protocol, hydrogen has to be cooled down to À30 C prior to the dispenser.This task was performed by the HRS built-in heat exchanger.

Results and discussion
Fig. 3a presents the results of a simulation compared against thermocouple measurements for Experiment 1.The pressure at the inlet nozzle, p inl , is shown in the background.The simulation is in good agreement with the measured temperatures.The 2-Norm of the error over the whole filling process is 1.3 K, which is of the same order as the precision of 1.1 K of the type K thermocouples.The model, however, is not capable of reproducing the quick temperature rise at the end of the filling phase, as it is a 3-dimensional effect involving the inlet nozzle.The temperature is slightly overestimated after the end of the filling, as the inlet mass flow rate is assumed to be zero.According to Eq. ( 7), the forced convection is then zero,    while in reality turbulent motions persist for some time and increase heat transfer to the inner wall.
The value of the parameter a in the forced Nusselt correlation, Eq (7), is determined by fitting the gas temperature curves to the thermocouple data, using the least square approach.The obtained correlation reads Nu for ¼ 0:32 |ffl{zffl} a Re 0:67 inl ; where the value of a is similar to results from [16,37].The evolution of the incoming hydrogen temperature T inl is split into two phases.During the first phase of Experiment 1, the temperature decreases nearly linearly during approximately 20 s.This can be attributed to the progressive mixing of the leftover gas in the pipes at ambient temperature and the cooled gas coming out of the heat exchanger.Afterwards during the second phase, the inlet temperature remains close to 245 K. Fig. 3b shows the temperature profile through the tank shell at different times during filling.The interface of the two material layers is clearly visible, with the interface radius being r int ¼ 135.4 mm.Heat is transported faster through the carbon wrapping layer, as indicated by the ratio a lin /a CFRP- ¼ 0.33 of heat diffusivities, see Table 3 for the thermal properties of the reference tank materials.
Interestingly, during such fast fillings (under 3 min), the temperature in the outer layer does not rise.The ambient temperature barely influences what happens inside the tank, the interesting parameter is then rather the starting temperature T init .However, it is common to consider hydrogen and tank shell initially in equilibrium with the atmosphere, i.e., T g ¼ T w ¼ T a , and it is in that sense that ambient conditions influence the filling processing.Fig. 4 highlights the importance of the different terms in the hydrogen energy balance, Eq. ( 3), during filling.Compression is responsible for the overall heat increase in the tank, which is counterbalanced by heat transfer to the inner wall and inflowing pre-cooled hydrogen.While heat increase due to compression is very fast, the cooling effect due to the other two terms require more time to be effective, leading to the two-phase temperature rise described previously.
Fig. 5 shows the ratio of forced, Eq. ( 7), to natural, Eq. ( 9), convection Nusselt number.It illustrates that forced convection is the dominant flow driver for the heat transfer.Natural convection can therefore be neglected during filling, as also assumed in [16,18].
The verified model was further used to investigate the influence of several parameters of the filling process on the temperature distribution in the tank shell.The critical variable is the temperature of the tank inner wall, which should not exceed 85 C (the highest temperature in the tank shell is always reached at the inner wall, i.e., at r ¼ 129.4 mm).The    impact of the initial temperature T init , the inlet temperature T inl , the inlet nozzle radius r inl , and the filling time t full on this critical variable is investigated.For this purpose, a reference case is considered, whose filling parameters are defined in Table 4.The previously described first phase of the inlet temperature evolution is modeled as linear decrease from the starting temperature of 290 K, to T inl after 30 s.Note that in order to modify the filling time, the APRR is not set according to the CEP-A70 protocol, but to an arbitrary choice for demonstrative purposes.Unless specified otherwise, from now on this set of parameters is considered for every filling scenario.Fig. 6 shows inner wall temperature evolutions during filling for different inlet temperatures; 6a for winter and 6b for summer conditions, with T init ¼ 270 K and T init ¼ 310 K respectively.
As expected, the lower the inlet temperature, the lower the maximum temperature reached at the inner wall.It is clear that reducing the cooling of the hydrogen ahead of the dispenser leads to excessive temperatures in the liner, as is the case with T inl ¼ 275 K, under summer conditions.Fig. 7a shows the inner wall temperature for different nozzle radii.In the 3 different cases presented, the mass flow rate and thermodynamic state of the inflowing hydrogen are the same, and therefore the inlet velocity is lower for the bigger inlet radii.Fig. 7b depicts the inlet velocity evolution during filling.Note that for the smallest inlet radius r inl- ¼ 1.5 mm, the Mach number reached Ma ¼ 0.13.While this value is relatively high, it is sufficiently small such that compressibility effects due to high speeds can be neglected.
As a result of those higher velocities, the overall heat transfer coefficient k g/w from the gas to the inner wall is much higher in the case of small inlet radii, as shown in Fig. 8.Note however that the heat flux ratio is not as high as the heat transfer coefficient ratio.According Eq. ( 5), the temperature difference DT g/w between the gas and the inner wall is smaller.Indeed, as k g/w is higher, the temperature of the inner wall quickly rises to the hydrogen temperature.In order to maximize heat transfer, and thus to keep DT g/w sufficiently high in the case of highly turbulent inflow, heat must be removed from the inner tank shell layer at a higher rate.The heat conductivity of the tank materials is then expected to have a significant impact on the peak temperature at the inner wall.
Fig. 9 presents the inner wall temperature evolution of the reference case, where the typical Biot number is Bi g/w ~600, and cases with lower and higher heat conductivity, thus higher and lower Biot numbers, respectively.The cases are modified by multiplying the thermal conductivity, and thus the thermal diffusivity, by different factors, resulting in different Biot numbers.While the specific heat capacity c p of the tank materials also has an impact on the temperature evolution of the system, the heat conductivity is chosen as the variable parameter to stay consistent with the Biot number definition, Eq. ( 13).
The inner wall temperature in the case of higher conductivity follows the reference case during the first phase of the filling, then rises slower until the end.In the opposite case of lower heat conductivity, the inner wall temperature rises slowly during the first phase, then its slope is similar to that of the reference case.In both cases, the temperature reached at the end of filling is lower.These two behaviors are made better visible in Fig. 10a and b, and are further explained by considering the evolution of heat fluxes.
Fig. 11 depicts the ratio of heat fluxes of cases with different Biot numbers.Note that in the four presented cases, the value of the heat transfer coefficient only changes up to a few percent, meaning that the differences of heat flux are primarily due to differences of DT g/w .
In the case of higher thermal diffusivity, T w,i follows the increase of T g , resulting in smaller heat flux at the beginning of the filling.However, due to the high diffusivity of the tank materials, heat is quickly transferred to outer layers of the tank shell.The temperature difference DT g/w then grows with time.
In the case of lower thermal diffusivity, more time is required for the inner wall temperature to increase.As the hydrogen temperature rises quickly in the first phase, so does DT g/w , and in turn _ Q g/w .The inner wall layer collects a lot of heat from the hydrogen due to this high heat flux, but the heat diffusion through the tank shell is much slower, and thus heat is stored in the liner layer.This stored heat then creates a bottleneck by accumulating, and so DT g/w and in turn _ Q g/w decrease.
In the four presented cases, whichever is the behavior of the heat flux, the ratio against the reference case is superior to 1, resulting in every case in a smaller peak temperature.However, it needs to be mentioned that the type IV tank considered here has a limited capacity, that is about 1.3 kg of hydrogen at 70 MPa.It is expected that for light-(3e5 kg) or heavy-(7e10 kg) duty tanks, as the amount of heat produced by compression is greater, a low thermal diffusivity in the tank shell would result in a bottleneck, therefore in a higher inner wall temperature.
Fig. 12 presents the inner wall temperature evolution for different filling times.The shorter the filling time, the higher the inlet velocity and in turn the higher the overall heat transfer coefficient.As previously discussed, for fast filling of light-duty tanks, the outer wall to ambient heat transfer has no influence.Indeed, due to slow heat diffusion, the outer wall temperature T w,o still is at its initial value.However, for slower   fillings the resistance to heat transfer comes from the limited heat transfer coefficient k w/a , as it is only subject to natural convection.The heat generated by compression is then trapped within the tank shell, taking up to several hours to cool down under natural ventilation conditions.For heavyduty tanks, the higher amount of heat that needs to be removed may lead to excessive temperatures.

Conclusions
A conjugate heat transfer model combining thermodynamic evolution of the fluid and energy equation in cylindrical coordinates in the solid tank shell has been developed.Experiments were conducted using a Type IV tank, and the model was assessed by comparing predicted temperature evolutions to experimental data, showing good agreement.
The verified model was further used for heat transfer analysis.It is well known that the critical parameter is the peak temperature in the tank shell, which should not exceed 85 C.
The temperature evolution of the gas is driven by the compression of hydrogen in the tank, while heat transfer from hydrogen to the tank inner wall is determined by forced convection.Apart from the temperature of the inflowing precooled hydrogen, the initial hydrogen temperature at the start of the filling is of major importance; for fast fillings more than ambient temperature per se.
Similarly, a smaller nozzle inlet radius induces higher heat transfer and thus leads to lower peak temperatures.However, this conclusion only is valid if the thermal diffusivity of the tank materials is high enough, such that heat gets quickly enough transported to the outer layers of the tank shell and from there to the atmosphere.
The study of the impact of thermal diffusivity shows that, for small-scale tanks, both very large and very low thermal diffusivities result in lower peak temperatures.Nonetheless, for heavy-duty tanks, as the amount of heat produced by compression of hydrogen is larger, the bottleneck created by reduced tank material diffusivity results in greater peak temperatures, eventually reaching the threshold value.
The extension of the proposed conjugate heat transfer model for larger tanks faces some issues [29].Indeed, as the aspect ratio L/D of the tank gets increased, the turbulent jet created by the nozzle is not reaching the end of the tank [38].The assumption of homogeneous hydrogen temperature breaks down, as

Fig. 1 e
Fig. 1 e Flow chart of FFTB algorithm.

Fig. 2 e
Fig. 2 e Scheme of the tank with the thermocouple tree.Thermocouple positions are indicated in red.

Fig. 3 e
Fig. 3 e Simulation of Experiment 1 filling conditions.Comparison between simulated hydrogen temperature and thermocouples (T1 -T8) data (a) and inner wall temperature profile evolution (b).The vertical dashed line (b) shows the boundary between the two tank materials.

Fig. 4 e
Fig. 4 e Impact of the different terms in the energy balance equation, Eq. (3), for Experiment 1. Term 1 refers to compression, _ Q g/w to heat exchange and Term 3 to cooling due to pre-cooled inflowing hydrogen.

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 7 ( 2 0 2 2 ) 2 3 0 6 0 e2 3 0 6 9

Fig. 5 eFig. 6 e
Fig. 5 e Ratio of forced to natural convection Nusselt number during filling for Experiment 1.

Fig. 7 e
Fig. 7 e Inner wall temperature evolution (a) and inlet velocity (b) for different inlet nozzle radii at T init ¼ 290 K.

Fig. 8 e
Fig. 8 e Heat transfer coefficient and heat flux ratios, computed for r inl ¼ 1.5 mm against r inl ¼ 6 mm.

Fig. 9 e
Fig. 9 e Inner wall temperature evolutions of the reference case with higher and lower Biot numbers, Bi g/w,ref ~600.

Fig. 10 e
Fig. 10 e Temperature profile across the tank shell over time for different Biot numbers, Bi g/w ¼ 3000 (a) and Bi g/w ¼ 120 (b).Black curves represent temperature isolines separated by 5 K.The vertical dashed red line indicates the materials boundary.

Fig. 11 e
Fig. 11 e Heat flux ratio, computed for different Biot numbers against reference case.

Fig. 12 e
Fig. 12 e Inner wall temperature evolution for different filling time, at T init ¼ 290 K.

Table 1 e
Filling parameters of the different experiments.

Table 2 e
Tank characteristics.

Table 3 e
Tank material properties.