About this Journal Submit a Manuscript Table of Contents
Abstract and Applied Analysis
Volume 2014 (2014), Article ID 583532, 8 pages
Research Article

Numerical Simulation of the Hydrogen Dispersion Behavior by a Parallel Characteristic Curve Method

1School of Engineering, Sun Yat-Sen University, Guangzhou 510275, China
2School of Pharmaceutical Sciences, Sun Yat-Sen University, Guangzhou 510006, China

Received 8 December 2013; Accepted 29 December 2013; Published 11 February 2014

Academic Editor: Hossein Jafari

Copyright © 2014 Qinghe Yao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


A parallel characteristic curve method is applied in domain decomposition system to simulate the dispersion behavior of hydrogen in this work. The characteristic curve method is employed to approximate the Navier-Stokes equations and the convection diffusion equation, and the feasibility of solving complex multicomponent flow problems is demonstrated by the numerical simulation of hydrogen dispersion in a partially open space. An analogy of the Boussinesq approximation is applied and numerical results are validated by comparing them with the experimental data. The dilution effect of ventilation is investigated. The transient behavior of hydrogen and the process of accumulation in the partially open space are discussed.

1. Introduction

Hydrogen is seen as a possible replacement to fossil fuels, which, besides being highly polluting, are fast being depleted. However, several obstacles must be overcome if hydrogen is to be used as a mainstream source of energy that safe storage is one of them. Because of the computation complexity of high Rayleigh number in hydrogen simulations, conventional numerical methods suffer from low convergence speed and poor stability and robustness. Because of the nonlinearity of these problems, product-type iteration methods are required [1, 2]. These methods occupy so much memory and computational time that they are difficult to be applied to large scale simulations. Based on the approximation of the material derivative along the trajectory of fluid particle, the characteristic curve method is natural in the simulation to physical phenomena and it is demonstrated to be unconditionally stable for a wide class of problems [37]. As a typical dispersion behavior of multicomponent flow, a methodology to simulate the hydrogen dispersion behavior efficiently is discussed in this work.

Safety is recognized as one of the most important issues in promoting hydrogen use. Hydrogen is odourless, colourless, and tasteless, so most human senses will not help detect a leak. It is no more or less dangerous than other flammable fuels, including gasoline and natural gas. Hydrogen is flammable and can behave dangerously under specific conditions. Accompanying a more extensive usage of hydrogen is the increased possibility of accidental release in the hydrogen infrastructure which comprises storage, bulk transportation and distribution, production, and utilization. However, hydrogen can be handled safely when simple guidelines are observed and the user has an understanding of its behavior. In particular, research on leaking hydrogen is important to prevent accidental ignition and to set the safety margin for leakage. Therefore, to ensure the safe usage of hydrogen it is necessary to predict and understand the characteristics of its leakage and dispersion. Over the last few decades, much research has been done on the evaluation of leak flow rate [8], dispersion behavior in residential areas [9], and the design of ventilation systems of hydrogen dispersion [10, 11]. Inoue et al. report the experimental data of a ventilation model [12]; Kanayama et al. report a numerical simulation to it by finite element method [13, 14]; however, the numerical results contain obvious oscillations and the maximum value of concentration is even higher than the value concentration at the inlet, which are contradicted with the experimental data and prevent it from being a better simulation.

The current study is to improve the simulation of hydrogen dispersion by a parallel characteristic curve method for its unconditional stability and to investigate the simulation of hydrogen dispersion behavior by doing this. It is expected to stabilize the numerical result and to have a better view about the recirculation zones by using a parallel characteristic curve method. To handle the problem caused by the nonlinear convective terms of flow problems, which result in the nonsymmetry of the stiffness matrix, a parallelized characteristic curve method for the domain decomposition method is used in this work. Compared with the traditional fashion is to employ some product-type methods as the iteration solver [14], the symmetry of the matrix enables computation problems with up to 30 million degrees of freedom (DOF) which can be solved [15]. In order to validate the solvability of dispersion behavior of hydrogen, the current computation results are compared with experimental results reported by Inoue et al. [12]. The transient dispersion behavior of hydrogen and several guidelines for safety in a ventilation model are discussed in this work.

The remaining sections are arranged as follows. Section 2 gives a brief description about the formulas and the parallel characteristic curve method. Section 3 describes the physical model, initial and boundary settings, and material properties. Numerical results and discussions are present in Section 4; finally, Section 5 gives concluding remarks in this research and points out the directions for the future study of dispersion behavior of hydrogen.

2. Formulation

2.1. Governing Equations

Let be the boundary of a three-dimensional polyhedral domain . is the Sobolev space of the first order and is the subspace of functions with zero mean value. Under the assumption that the flow field is incompressible, viscous, and laminar, the solving of the model can be summarized as finding such that, for any , the following set of equations holds: where is the velocity [m/s], is the pressure [N/m2], is the density (const.) [kg/m3], is the gravity [m/s2], is the concentration expansion coefficient [—], is the mass concentration [mass %], and is the stress tensor [N/m2] defined by with the Kronecker delta and the viscosity [kg/ms].

An initial velocity is applied in at . Dirichlet boundary conditions and Neumann boundary conditions are also applied, where are the outward normal direction to .

The governing equation of leaking hydrogen dispersion is to find such that where is the diffusion coefficient [m2/s] and is the source term [K/s]. An initial concentration is applied in at . Dirichlet and Neumann boundary conditions are set by respectively, where is the outward normal derivate to .

2.2. A Parallel Characteristic Curve Method

A characteristic curve method is applied in parallel manner to the nonlinear terms in (1) and (5), which works as follows.

Let be the time increment and the total step number, and thus stands for time step . Let be the current position of a particle, its position at can be approximated by where is the fluid velocity of position at and denotes a first order approximation; see [16]; therefore, the material derivative in (1) and (5) at can be written as In (1) and (5), is and , respectively.

A two-way coupling scheme is employed in this research to couple the solving process of (1) and (5), as can be seen in Figure 1.

Figure 1: The coupling scheme.

It should be noted that the element information term may cause some difficulty in the solving of governing equations (1) and (5), especially in this work when all the elements information is distributed on different processor elements (PEs) and the searching is done in parallel. Fortunately, the same is used for both equations; therefore, the searching results can be shared by them in one nonstationary step; see Figure 1.

The element searching algorithm requires a globalwise element information to determine the position of one particle in the previous time step. However, in the parallel system used in this work, the whole domain is split into several parts, one processor element (PE) works only on the current part and it does not contain any element information of other parts. Each part is further divided into many subdomains and the domain decomposition is performed by the PE in charge of the part. This parallelity causes a computational difficulty: for each time step, the particle is not limited within one part; see Figure 2; therefore, exchanging the data between different PEs is necessary, which demands the PEs to communicate in subdomainwise computation.

Figure 2: Searching over parts.

The number of total elements in one subdomain may be different, which means that some point to point communication techniques, such as MPI_Send/MPI_Recv or MPI_Sendrecv in MPICH, cannot be used in elementwise computation. In the previous research [17], a global variable to store all the old solutions is constructed. This method maintains a high computation speed but costs a huge memory usage. To reduce the memory consumption, a request-response system is used in [18]. In this work, the request-response system is further optimized by starting the request earlier in each loop; thus, less time is spent in communication among PEs.

3. Modeling

3.1. Geometry and Parameters

Hydrogen has a rapid diffusivity; it dilutes quickly into a nonflammable concentration when released. Therefore, to become a fire hazard, hydrogen must first be confined, but as the lightest element in the universe, confining hydrogen is very difficult. Industry takes these properties into account when designing structures where hydrogen will be used. The designs help hydrogen escape up and away from the user in case of an unexpected release.

A ventilation model with partially open space considered in this work is shown in Figure 3, used by several researchers [10, 12, 13] to assess the risk of an accident caused by a hydrogen leak. Hydrogen enters from an inlet at the base near one end of the model; ventilation is through a roof vent and a door vent near the opposite end. Four sensors are placed within the model; their locations are defined as the table embedded in Figure 3.

Figure 3: The ventilation model.

During hydrogen dispersion, the discrepancy in concentration is one of the main drive forces of flow motion, which give us the basic idea to use an analogue of Boussinesq approximation in thermal convection problems. According to the ideal gas equation of state, the concentration expansion coefficient keeps constant for hydrogen dispersion [13]. This finding makes it extremely convenient for us to apply the current solver to this problem. Using the notations defined in Section 2.1, the material properties used in this work are given by Table 1.

Table 1: The material properties.

3.2. Initial and Boundary Conditions

To be consistent with the experiment reported by [12], hydrogen leaks into the inlet at the speed of 0.02 [m/s] in the vertical direction, with a mass concentration of 6.94% (considering the difference between the density of air and hydrogen). At the roof, the hydrogen is discharged outside freely. At the door, the air is supposed to go in due to the pressure discrepancy between inside and outside of the door. All other boundaries are considered as endwalls.

Let , , and denote the boundary of inlet, roof and door, respectively; Dirichlet boundary conditions are set as follows:

Natural ventilation is imposed at both the roof and the doors. Neumann boundary conditions are set as follows:

Initial velocity and concentration are set as follows:

It can be seen that Rayleigh number in this research is and it is believed to be the main reason of oscillations in [13, 14].

4. Numerical Results

In order to validate the solvability of the current scheme on hydrogen dispersion behaviour, a comparison with the experimental data reported by (12) is made and the total computation time is set to 1500 s. To avoid the effect of the stationary judgment, the criterion of the stationary stage is set to and controlled by an element based norm, compare [19]. The convergence judgment of each time step is made by the Euclidean norm [20], and the criterion is .

The ADVENTURE_CAD and ADVENTURE_Metis [21] are used to create the geometry and mesh, and the local density of mesh around the sensors in Figure 3 is increased to have a better simulation. As an initial attempt of this research, an unstructured mesh with 1,124,294 elements is created for this work, as is shown in Figure 4.

Figure 4: The meshing.

The computation results of this work are compared with the experimental results reported by [12]. Figure 5 compares the computational hydrogen volumetric concentration values at the nodes closest to the four sensor positions against experimental values in [12]. The hydrogen concentration is measured in terms of volumetric concentration, where 0% indicates that the entire volume is occupied by air and 100% represents the opposite.

Figure 5: Comparison of numerical and experimental results.

It can be seen for Figure 5 that the current numerical results agree with the experimental data very well, especially sensor 2 and sensor 3, which are located at the top of the ventilation model; see Figure 3. The computational results at sensor 1 and sensor 4 are also closer to the experimental data and the trend fits that of the experimental data very well. Around the inlet and the door, the flow field is complicated and oscillations are also viewed at sensor 1 and sensor 4. Compared with numerical results in [13, 14], the current numerical results are more stable and closer to the experimental data and thus are more reliable.

Hydrogen has a wide flammability range (4–74% in air) and the minimum energy required to ignite hydrogen is very low (0.02 mJ, 10% of the minimum energy required to ignite gasoline vapour); as a result of this, the leakage of hydrogen in a confined space introduces the possibility of accidental ignition, which, in the worst case, may result in an explosion. Such safety issues necessitate special facilities if extensive physical experiments are to be carried out safely and frequently, further justifying the use of computational simulations. In Figure 6, the isosurface of volumetric concentration at 4% is shown, presenting the boundary of flammability inside the partially open space.

Figure 6: The isosurface of volumetric concentration at 4%.

Hydrogen leaks at a constant speed with a constant concentration at the inlet, as is shown by (12). It can be seen from Figure 6 that the 4% volumetric concentration isosurface is getting lower and lower toward the bottom of the ventilation model; however, after 500 s, the 4% volumetric concentration isosurface does not get lower obviously and the height of the isosurface is kept till the end of the computation. Thus, the reign below the 4% volumetric concentration isosurface is relatively safe. This finding could help in case if the hydrogen leakage occurs; the area around the place of accident can be classified by its safety, and instruction can be given differently to avoid the so-called “secondary accident.”

Velocity profiles are given in Figure 7. Through these figures, several characters of hydrogen leakage inside the partially open space can be observed as follows.(1)Hydrogen flows vertically and then spreads to the upper roof vent after reaching the ceiling; therefore, the sensors near the ceiling keep a high level of hydrogen concentration; the sensors are sorted as by descending order of concentration value.(2)Air goes inside the model from the door vent, due to the changes in pressure field; that is to say, hydrogen concentration near the horizon of the door vent keeps a low level during the dispersion, as is shown in Figure 6.(3)The distribution of hydrogen concentration does not change a lot after 500 s; in some sense, a stationary state is achieved.(4)The flow field around the inlet and the door is very complicated due to the inflow of air, as the moving air is a driving force that cannot be neglected. Streamlines at 50 s are displayed in Figure 8 to demonstrate the flow pattern in the ventilation model.

Figure 7: Velocity vectors.
Figure 8: Streamlines.

As can be seen in Figure 8, the flow fields around the inlet and the door are affected by air also. The great distortion near these areas indicates that some turbulent flow appears during the computation, which is consistent with the oscillation in Figure 5.

The computation model contains 187,382 nodes and mesh size is about 0.03. It is analyzed parallelly by the scheme described in Section 2.2. The total DOF was 936,910. A Linux PC cluster with 20 PCs in Kyushu University was mainly utilized to carry out the computation. For each PC the configuration of hardware is as follows:CPU: Intel(R) Core(TM) i7 920@2.67GHz,Memory: 12 [GB].

Due to the unconditional stability of the parallel characteristic curve method, was set to 0.1 s. The computation took about 10 hours on this small PC cluster.

5. Conclusions

A ventilation model of hydrogen dispersion is numerically simulated in three dimensions by a parallel characteristic curve method in a domain decomposition system in this work. The main conclusions can be summarized as follows.(1)By using a parallel characteristic curve method, the results are more stable that the conventional computation result is gotten by product-type solvers; the current results are closer to the experimental data and are more reliable; the comparison of current computation results and the experimental data convinces us of the solvability of the current scheme in hydrogen dispersion problems.(2)The door plays an important role in preventing the accumulation of hydrogen around the bottom of the ventilation model. Air goes inside the model and the dilution effect appears.(3)Turbulent flow appears around the inlet and the door, and recirculation zones are found inside the model;(4)A stationary state is reached after about 500 s;

As is mentioned above, to investigate the effect of driven force by air more precisely, some stabilization techniques and a finer mesh may be required, which will be a future work to us.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publishing of this paper.


This work was supported by the National Science Foundation of China (NSFC), Grants 11202248, 91230114, and 11072272. China Postdoctoral Science Foundation, Grant 2012M521646, and the Guangdong National Science Foundation, Grant S2012040007687.


  1. H. Kanayama, K. Komori, and D. Sato, “Development of a thermal convection solver with hierarchical domain decomposition method,” in Proceedings of the 8th World Congress on Computational Mechanics and the 5th European Congress on Computational Methods in Applied Sciences and Engineering, Venice, Italy, 2008.
  2. H. A. van der Vorst, Iterative Krylov Methods for Large Linear Systems, Cambridge University Press, 2003. View at MathSciNet
  3. M. Bercovier, O. Pironneau, and V. Sastri, “Finite Elements and characteristics for some parabolic-hyperbolic problems,” Applied Mathematical Modelling, vol. 7, no. 2, pp. 89–96, 1983. View at Scopus
  4. J. Douglas Jr. and T. F. Russell, “Numerical methods for convection-dominated diffusion problems based on combining the method of characteristics with finite element or finite difference procedures,” SIAM Journal on Numerical Analysis, vol. 19, pp. 871–885, 1982.
  5. O. Pironneau, “On the transport-diffusion algorithm and its applications to the Navier-Stokes equations,” Numerische Mathematik, vol. 38, no. 3, pp. 309–332, 1982. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  6. T. F. Russell, “Time stepping along characteristics with incomplete iteration for a Galerkin approximation of miscible displacement in porous media,” SIAM Journal on Numerical Analysis, vol. 22, no. 2, pp. 970–1013, 1985. View at Zentralblatt MATH · View at MathSciNet · View at Scopus
  7. O. Pironneau, Finite Element Methods for Fluids, John Wiley and Sons, 1989. View at MathSciNet
  8. R. W. Schefer, W. G. Houf, C. San Marchi, W. P. Chernicoff, and L. Englom, “Characterization of leaks from compressed hydrogen dispensing systems and related components,” International Journal of Hydrogen Energy, vol. 31, no. 9, pp. 1247–1260, 2006. View at Publisher · View at Google Scholar · View at Scopus
  9. H. A. Olvera and A. R. Choudhuri, “Numerical simulation of hydrogen dispersion in the vicinity of a cubical building in stable stratified atmospheres,” International Journal of Hydrogen Energy, vol. 31, no. 15, pp. 2356–2369, 2006. View at Publisher · View at Google Scholar · View at Scopus
  10. K. Matsuura, H. Kanayama, H. Tsukikawa, and M. Inoue, “Numerical simulation of leaking hydrogen dispersion behavior in a partially open space,” International Journal of Hydrogen Energy, vol. 33, no. 1, pp. 240–247, 2008. View at Publisher · View at Google Scholar · View at Scopus
  11. J. Choi, N. Hur, S. Kang, E. D. Lee, and K.-B. Lee, “A CFD simulation of hydrogen dispersion for the hydrogen leakage from a fuel cell vehicle in an underground parking garage,” International Journal of Hydrogen Energy, vol. 38, pp. 8084–8091, 2013.
  12. M. Inoue, H. Tsukikawa, H. Kanayama, and K. Matsuura, “Experimental study on leaking hydrogen dispersion in a partially open space,” Journal of Hydrogen Energy Systems Society of Japan, vol. 33, pp. 32–43, 2008.
  13. H. Kanayama, H. Tsukikawa, S. B. Ndong-Mefane, and O. Sakuragi, “Finite element simulation of hydrogen dispersion by the analogy of the Boussinesq approximation,” Journal of Computational Science and Technology, vol. 2, pp. 643–654, 2008.
  14. H. Kanayama, H. Tsukikawa, and I. Ismail, “Simulation of hydrogen dispersion by the domain decomposition method,” Japan Journal of Industrial and Applied Mathematics, vol. 28, no. 1, pp. 43–53, 2011. View at Publisher · View at Google Scholar · View at Scopus
  15. Q. H. Yao, H. Kanayama, M. Ognio, and H. Notsu, “Incomplete balancing domain decomposition for large scale 3-D non-stationary incompressible flow problems,” in Proceedings of the 9th World Congress on Computational Mechanics and 4th Asian Pacific Congress on Computational Mechanics, 2010.
  16. O. Pironneau, Finite Element Methods for Fluids, Masson, Paris, France, 1989. View at MathSciNet
  17. Q. Yao, H. Kanayama, H. Notsu, and M. Ogino, “Balancing domain decomposition for non-stationary incompressible flow problems using a characteristic-curve method,” Journal of Computational Science and Technology, vol. 4, pp. 121–135, 2010.
  18. Y. Qing-He and Z. Qing-Yong, “Investigation of the contamination control in a cleaning room with a moving AGV by 3D large-scale simulation,” Journal of Applied Mathematics, vol. 2013, Article ID 570237, 10 pages, 2013. View at Publisher · View at Google Scholar
  19. Q. Yao and Q. A. Zhu, “Pressure-stabilized Lagrange-Galerkin method in a parallel domain decomposition system,” Abstract and Applied Analysis, vol. 2013, Article ID 161873, 13 pages, 2013. View at Publisher · View at Google Scholar
  20. Q. Yao -H and Q. Zhu -Y, “Investigation of the contamination control in a cleaning room with a moving AGV by 3D Large-Scale simulation,” Journal of Applied Mathematics, vol. 2013, Article ID 570237, 10 pages, 2013. View at Publisher · View at Google Scholar
  21. http://adventure.sys.t.u-tokyo.ac.jp/.