ǫǰǡǟ
Vol. 2: )3125*Ą38ġ43

An Attempt of Simulating the Real Time Filling of H2 Cylinder at 70MPa

An Attempt of Simulating the Real Time Filling of H2 Cylinder at 70MPa
M. M. A. Alam1*, T. Setoguchi1, S. Matsuo2
1

Institute of Ocean Energy, Saga University, 1, Honjo-machi, Saga, 840-8502, Japan.
2
Department of Advanced Technology Fusion, Saga University, Japan.

This report presents the computational fluid dynamics (CFD) studies on real time filling of H2 gas in cylindricaltank at 70MPa. A standard type IV, 39.1L vessel composed of plastic and Carbon Fibered Reinforced Plastic (CFRP) for
liner and laminate is used in the present study. A density based algorithm of ANSYS Fluent 13.0 is applied to solve thetime dependent 3D RANS equations with SST k-Z turbulence model. User defined real gas model is used for density
computation to model compressibility effects. The energy equation for the working gas is coupled to the energy
equation for the solid walls in order to estimate the convective heat transfer from the compressed gas to tank wall.
Duration of the tank filling of 163 s which is almost the same to real time filling at refueling station, is used in thepresent work.
Key Words : compressible flow, convection, hydrogen tank, numerical simulation, real gas

1. Introduction
As a sustainable and green energy carrier, hydrogen gas can meet the current needs of energy without sacrificing thecultural and economic growth, the environment (global worming, climate change, environmental degradation, etc.) or theability of future generations to live even better than we do today [1]. Hydrogen powered vehicles highlight the need of storinghydrogen (H2) gas for its clean and efficient conversion to energy. At present, pressurised gas cylinders represent the simplest
and most mature technology for H2 storage whose feasibility has already been demonstrated at working pressure up to 70MPa[2]. Pressurised cylinders are also the preferred near-term option for direct-hydrogen and hybrid fuel cell vehicles. Refuelling
vehicle cylinders involves transferring high-pressure hydrogen gas from the fuelling station tanks through a dispenser into thevehicle cylinder. The dispenser controls the rate of filling, ensures safe operation, and delivers the rated mass of gas [2]. The
total time required to refuel a cylinder is important for consumer acceptance and it can have significant financial
consequences in commercial applications. However, rapid filling of hydrogen increases the gas temperature. Current
regulation imposes limits on the maximum average gas temperature allowed within the cylinders for vehicular applications;the average H2 gas temperature cannot exceed 85oC [4-6]. Therefore, monitoring and feedback of temperature data inside thevessel and on the tank wall are essential during the real time filling.
Some experimental and numerical works may be found in literature. To evaluate temperature distribution within thecylinder wall, a simple approach consists in coupling a single gas temperature evolution prediction and a one-dimensional
conduction calculation within the wall [7]. Such a method was also applied in Air Liquide R&D [8] and allowed tounderstand the influence of filling parameters. Experimental studies using measurement devices placed inside the vessel havealso been carried out, allowing evaluating the temperature in few points [8-10] or in a more extended region using a large
matrix of sensors [3]. CFD simulations have also been performed to 35MPa filling for a duration inferior to one minute[11-13]. In the present work, CFD simulations are conducted on the real time filling of H2 gas in cylindrical tank at 70MPa.
Real gas effects are incorporated by using user defined real gas model to account for the deviation of hydrogen gas from ideal
gas behaviour. Filling time of 163 s that mostly needed at practical refuelling station is used as the duration of filling in thepresent work.
2. Numerical model and computational methodology*

Date of Manuscript Acceptance 2014.7.31
E-mail: dralam@me.saga-u.ac.jp

38

39

M. M. A. Alam, T. Setoguchi, S. Matsuo

2.1 Governing equations
The simulation of real time filling of hydrogen gas in pressurized cylindrical vessel requires the solution of a system ofgoverning equations in the unsteady form. The governing equations used in the present work are time dependent
compressible Reynolds and Favre averaged Navier-Stokes (RANS) equations. For the closure of the non-linear term from theconvective acceleration and Reynolds stress, the well known SST k-Z turbulence model [14-16] is used.
The conservation of mass and momentum do not apply at the tank walls. The energy equation for the walls is given by;w
Uhwt  k s T(1)where the sensible enthalpy is given by;T³ c ps dT

h

(2)

Tref

The energy equation for the working gas is coupled to the energy equation for the solid walls.
In real time filling of hydrogen tank involving high pressure and temperature compressibility effects are significant, and the
behavior of hydrogen gas is expected to deviate from ideal gas. Ideal gas simulation may not be sufficient, and therefore,
Sakoda et al. [17] equation of state (EOS) is used in the present work for real gas modeling. The functional form of the virial
EOS proposed by Sakoda et al. [16] can be expressed by,
Z 1 B T U  C T U 2  D T U 3
b1  b2 exp TrBT

c1  c2Tr

CT
DT

d1Tr

 k3

1

k

(3)

 b3Tr 1  b4Tr

 c3Tr

 k2

(4)

 k4

(5)

 k5

(6)

In the above equations, B, C and D are the second, third, and forth virial coefficients, respectively, and Tr =T/Tc, where Tc isthe critical temperature (33.145K for H2). The coefficients k1 – k5, b1 – b4, c1 – c3, and d1 weredetermined using non-linear
least-squares fitting to minimize the sum of squares,
SS

I

J

K

i 1

j 1

k 1

2
2
2
¦ wi2 Z Exp, i  Z Cal , i  ¦ w 2j BExp, j  BCal , j  ¦ wk2 C Exp, k  C Cal , k(7)

where w is the weight, I, J, and K are the number of density data, the second virial coefficients, and the third virial coefficients,
respectively. The determined coefficients are listed in Table 1. The specific heat, cpg, thermal conductivity, kg, speed of sound,a, surface tension,V and viscosity, P are computed employing equations proposed by Leachmann et al. [18].
Table 1 Coefficients in the virial EOS for Eqs. 3-6
b1
b2
b3
b4
c1
c2
c3

0.161 777×10-1-0.102 798×100
0.129 501×100
0.628 764×10-1
0.288 234×10-3
0.203 867×10-2-0.779 969×10-3

d1
k1
k2
k3
k4
k5

0.716 809×10-4
0.0612
2.71
2.53
4.85
1.09

2.2 Numerical methods
The density based solver of ANSYS Fluent 13.0 [19] is used to solve coupled three dimensional system of equations.
Hydrogen gas is selected as working fluid from ANSYS Fluent 13.0 database. Spatial derivatives are discretized by 3rd order
MUSCL scheme that is a blend of central differencing and second-order upwind schemes. Second-order central difference
scheme is used for viscous terms. Second-order implicit time integration scheme is used to discretize the temporal derivatives.
Convective heat transfer from the gas to tank wall is estimated by wall discritization for local heat transfer coefficient.
Radiation heat transfer between the gas and wall is assumed to be negligible since temperature difference between the gas andwall is very small.

3:

An Attempt of Simulating the Real Time Filling of H2 Cylinder at 70MPa

2.2 Computational conditions
Schematics of the computational geometry and grids are presented in Fig. 1. Type IV vessel having the volume of 39.1L,
inner radius of 0.28mm, liner and laminate thicknesses are of 3 and 22 mm, respectively, is used in the present computational
work. Model geometry is divided into two domains: fluid domain filled with hydrogen gas and solid domain involving liner
and laminate regions on the tank wall and inlet tube. Liner is made of plastic and laminate is constructed from carbon fiber
reinforced plastics (CFRP). Properties of liner and laminate materials used in the computation are shown in Table 2. Inlet tube
is protruding into the tank up to a distance of 70 mm. Structured clustered grid is used consisted of 222,752 hexahedral
elements in fluid domain and 109,156 hexahedral elements in solid domain. Grids are densely clustered in near wall regions
to capture flow features in boundary layers. Enough nodes are located in the structured grid in the solid region to computewall heat transfer accurately.
Pressure inlet boundary condition is specified at tube inlet. Inlet pressure and temperature are varied with time, as shown in
Fig.2, to replicate the experimental conditions. Curve fitted experimental data of pressures and temperatures are read into thenumerical model as a UDF (user defined function). No slip boundary conditions are specified at inner walls. Energy equation
for gas is coupled to energy equation for walls. A constant heat transfer coefficient of 10 W/m2 K is specified at the outer wall.
As the filling completed in a couple of minutes, a constant ambient temperature of 283.65 K is specified. Pressure andtemperature of the gas is assumed to be uniform within the tank before filling. The flow field is initialized with initialconditions of hydrogen gas in the vessel i.e. pressure 1.32MPa and temperature 279.83 K. Initial temperature of the tankwalls is assumed to be same as the initial temperature of gas.
For optimum time step size, an initial guess 'T=10-4s, is made based on the initial gas velocity and minimum cell size.
Table 2 Properties of liner and laminate materials
Material

Thermal conductivity (W/m.K)

Density (kg/m3)

Linerplastic

1578

1.17

1286

Laminate

CRFP

1075

1.14

1375

Fig. 1 Schematics of tank geometry and computational grids

Inlet pressure
Inlet temperature

60

20

40
0
20

-20

0

-20

Tempreature, oC

40

80

Pressure, MPa

Type IV tankSpecific heat (J/kg.K)

0

40

80

120

Time, sFig. 2 Inlet pressure and temperature

-40
160

41

M. M. A. Alam, T. Setoguchi, S. Matsuo

Subsequently, the solution is reiterated using reduced time step sizes. For grid convergence investigation, initially hexahedral
coarse mesh is used and then a finer mesh with double the original mesh density is created. The grid is also refined in theshear layer between the incoming jet and gas already present in the tank. The problem is reiterated using the time step size of'T=10-4s and refined mesh. However, shorter time steps and finer meshes turned out to be computationally expensive.
Therefore, the selected time step and mesh size are used in all subsequent simulations.
3. Results and discussion
Velocity vectors along the tank mid plane (xy-plane) at 8 s and 12 s from the commencement of filling are presented in Figs.
3(a) and 3(b), respectively. Gas velocity and velocity gradients are found to be highest in the vicinity of tank axis. It isobserved that the incoming jet strikes the bottom of the tank, and flow turns along the rear surface towards the inlet forming asecondary flow region on either sides of inlet jet. These secondary flow regions shift positions and vary in size with respect totime during the fill. Flow velocity gradually decreases, as shown in Fig. 3(b), with the increase in pressure as density of gas
increases as the tank gets filled up. Flow at the inlet is turbulent throughout. Development of flow field within the tank isdominated by the structure of turbulent jet flowing from the inlet. More turbulence near the inlet results in enhanced heattransfer to the wall in this region.
Contours showing in Fig. 4 are illustrated the variations of density along the mid plane of the flowfield. From the contour
sequences, it is observed that density of gas increases with the increase in tank pressure. A nearly conical inlet jet plume attransonic speed can be observed and initially, it extends from the inlet to about tank mid-section in the longitudinal direction;
while the jet plume spreads about the mid-section in the radial direction. The plume increases initially and decreases again insize as filling proceeds due to compression of gas and ‘growth’ of the bulk region. However, the jet plume shows symmetriccharacteristics during filling of tank.
Temperature variation during filling is presented in Fig. 5. As per the results, tank may be divided into two regions: thenear inlet region with a high temperature gradient and a bulk flow region with no significant change in temperature. A nearly
conical inlet gas plume can be observed and initially, it extends from the inlet to about tank mid section in the longitudinal
direction. Local temperatures in the plume region are significantly lower than mean gas temperature. High temperature
gradient in the vicinity of inlet is due to the influence of the cooler incoming gas. Also high level of turbulence due to(a) t =8 s(b) t =12 sFig. 3 Velocity vectors on xy-planet =8 sFig. 4 Density contours on xy-planet =12 s42

An Attempt of Simulating the Real Time Filling of H2 Cylinder at 70MPa

t =4 st =12 sFig. 5 Temperature contours on xy-plane

350

Temperature, K/LQQHU

Temperature, K/DPLQDWH/LQQHU

350

325

300

275

0

0.2 0.4 0.6 0.8 1 1.2
Normalized distance, [m/m]/DPLQDWH

4s
8s

325

12 s

300

275
-0.2

0 0.2 0.4 0.6 0.8 1
Normalized distance, [m/m](a) Front end(b) Rear end

Fig. 6 Temperature distribution on vessel wall

recirculating gas mixing with high velocity inlet jet facilitates a high rate of heat transfer to the walls. Temperature is more orless same outside the inlet plume region.
Figures 6(a) and 6(b) show the temperature-time histories along wall at the front end and rear end of tank, respectively.
High temperature gradient is experienced in the liner, while the temperature gradient gradually diminishes along the width,
and reaches to zero after one-third of the laminate section. Due to the low thermal conductivity of CFRP, temperature of outersurface of the laminate has not recorded any appreciable change at any location during filling. However, comparing with theother locations a significant increase in temperature is recorded at the rear end of the tank.
7. Conclusion
A 3D numerical model for simulating the real time filling of hydrogen gas into cylindrical tank at 70MPa has beendeveloped. Thermo-fluid dynamic characteristics of hydrogen gas during filling have been investigated. Time dependent
compressible RANS equations were solved along with SST k-Ȧ turbulence model. Real gas equation of state was consideredto accurately compute the density and temperature of the gas flowfield. The numerical model predicted the temperature fieldwith a good accuracy. A conical jet plume was observed near tube exit region initially, and found extending and spreading toabout tank mid-section in the longitudinal and radial direction, respectively. Size of the jet plume was decreasing in size as
filling proceeds, and found increasing of bulk region. Temperature distribution along tank wall was obtained.
References(1) Scott, D., “Template for Sustainability”, International Journal of Hydrogen Energy, Vol. 29, No. 1 (2004), pp. 17-22.
(2) Colom, S., Weber, M., and Barbier, F., “Storhy: A European development of composite vessels for 70MPa Hydrogen storage”,
World Hydrogen Energy Conference, Brisbane (2008).

43

M. M. A. Alam, T. Setoguchi, S. Matsuo(3) Dicken, C.J.B., and Merida, W., “Measured effects of filling time and initial mass on the temperature distribution within ahydrogen cylinder during refuelling”, Journal of Power Sources, Vol. 165 (2007), pp. 324–336.
(4) ISO. Gaseous hydrogen and hydrogen blends: land vehicle fuel tanks. ISO/TS 15869. Geneva, Switzerland: International
Organization for Standardization (2009).
(5) SAE J2579, Vehicular hydrogen systems, USA: Society of Automotive Engineers (2008).
(6) JARIS 001, Technical standards for containers for compressed-hydrogen vehicle fuel device, Japan: Japan Automobile
Research Institute (2004).
(7) Khan, M.T.I., Monde, M. and Setoguchi, T., “Hydrogen gas filling into an actual tank at high pressure and optimization of its
thermal characteristics”, Journal of Thermal Science, Vol.18 (2009), pp. 235-240.
(8) Werlen, E., Pisot, P., Barral, K. and Renault, P., “Thermal effects related to H2 fast filling in high pressure vessels depending on
vessels types and filling procedures: modelling, trials and studies”, European Energy Conference, September, 2003.
(9) Liu, Y.L., Zhao, Y.Z., Zhao, L., Li, X., Chen, H.G., Zhang, L.F., Zhao, H., Sheng, R.H., Xie, T., Hu, D.H. and Zheng, J.Y.,“Experimental studies on temperature rise within a hydrogen cylinder during refueling”, International Journal of Hydrogen
Energy, to be published.
(10) Lee, S.H., Kim, Y.G., Kim, S.C. and Yoon, K.B., “Temperature change of a type IV cylinder during hydrogen fueling process”,
Proceeding of The 3rd International Conference on Hydrogen Safety, Ajaccio (2009).
(11) Dicken, C.J.B., and Merida, W., “Modeling the transient temperature distribution within a hydrogen cylinder during refueling”,
Numerical Heat Transfer, Part A, Vol. 53 (2008), pp. 685-708.
(12) Suryan, A., Kim, H.D. and Setoguchi, T., “Comparative study of turbulence models performance for refueling of compressedhydrogen tanks”, International Journal of Hydrogen Energy, to be published.
(13) Setoguchi, T., Alam, M.M.A., Monde, M., and Kim, H.D., “Characteristics of turbulent confined jets during fast filling of H2tank at high pressure, International Journal of Aeroacoustics, Vol. 12, No. 5-6 (2013), pp. 455-474.
(14) Menter, F.R., “Zonal Two Equation k-Ȧ Turbulence Models for Aerodynamic Flows”, AIAA Paper 93-2906 (1993).
(15) Menter, F. R., “Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications”, AIAA Journal, Vol. 32 (1994),
pp. 269-289.
(16) Menter, F.R., Kuntz, M. and Langtry, R., “Ten Years of Industrial Experience with the SST Turbulence Model”, In: Hanjalic,
K., Nagano, Y. and Tummers, M., editors, Turbulence, Heat and Mass Transfer, 4, Begell House Inc., (2003), pp. 625-632.
(17) Sakoda, N., Shindo, K., Motomura, K., Shinzato, K., Kohno, M., Takata, Y. and Fujii, M., “Burnet PVT Measurements ofHydrogen and the Development of a Virial Equation of State at Pressure up to 100MPa”, International Journal of
Thermophysics, Vol. 33 (2012), pp. 381-395.
(18) Leachmann, J.W., Jacobsen, R.T., Penocello, S.G. and Lemmon, E.W., “Fundamental Equations of State for Parahydrogen,
Normal Hydrogen, and Orthohydrogen”, Journal of Physical Chemical Reference Data, Vol. 38, No.2 (2009), pp. 721-748.
(19) ANSYS Inc., Fluent user’s guide, USA, (2013).

