Large Scale Simulation of Hydrogen Dispersion by a Stabilized Balancing Domain Decomposition Method

The dispersion behaviour of leaking hydrogen in a partially open space is simulated by a balancing domain decomposition method in this work. An analogy of the Boussinesq approximation is employed to describe the connection between the flow field and the concentration field. The linear systems of Navier-Stokes equations and the convection diffusion equation are symmetrized by a pressure stabilized Lagrange-Galerkin method, and thus a balancing domain decomposition method is enabled to solve the interface problem of the domain decomposition system. Numerical results are validated by comparing with the experimental data and available numerical results. The dilution effect of ventilation is investigated, especially at the doors, where flow pattern is complicated and oscillations appear in the past research reported by other researchers. The transient behaviour of hydrogen and the process of accumulation in the partially open space are discussed, and more details are revealed by large scale computation.


Introduction
With the development of hydrogen-fuelled vehicles, it becomes clearer and clearer that several obstacles must be overcome if hydrogen is to be used as a mainstream source of energy. As more usages of hydrogen are explored, the possibility of accidental release in the hydrogen infrastructure increases, which comprises storage, bulk transportation and distribution, production, and utilization. Hydrogen is flammable and can behave dangerously under specific conditions; however, hydrogen can be handled safely when simple guidelines are observed and the user has an understanding of its behaviour. Hydrogen is odourless, colourless, and tasteless, and most human senses will not help to detect a leak; therefore, to prevent accidental ignition and set the safety margin for leakage, it is necessary to predict and understand the characteristics of its leakage and dispersion. It is difficult to visualize the hydrogen dispersion by experiment in case of hydrogen leaks, because of its low kinematic viscosity and high diffusibility and risk. As such, clarifying the hydrogen dispersion with numerical simulation becomes important [1,2].
On hydrogen dispersion problems, the evaluation of leak flow rate [3], the dispersion behaviour in residential areas [4], and the design of ventilation systems [5,6] have been reported. Inoue et al. report the experimental data of a ventilation model [7], and Kanayama et al. report a numerical simulation to it by finite element method [1,8]; however, the numerical results contain obvious oscillations which prevent it to be a better simulation. Because of the computation complexity of high Rayleigh number in the modelling hydrogen dispersion, conventional numerical simulation methods suffer from low convergence speed, poor stability, and robustness [9,10]. These methods occupy too much memory and computational time to be applied to large scale simulations. By approximating the material derivative along the trajectory of fluid particles, the Lagrange-Galerkin method is reported to be unconditionally stable for a wide class of problems [11][12][13][14][15]. Moreover, the linear systems of Navier-Stokes equations and the convection diffusion equation are symmetrized and a balanced domain decomposition method is enabled to solve the interface problem of domain decomposition system [16].
The current study is to improve the simulation of hydrogen dispersion by a balancing domain decomposition method, which has shown its effectiveness for incompressible flow problems [17,18]. A variation of balancing domain decomposition is proposed and the efficiency is reported for simulation of hydrogen dispersion in this work. Compared with the traditional fashion of employing some product-type methods as the iteration solver [8], computation problems with up to 30 million degrees of freedom (DOF) can be solved on small Linux clusters by using the balancing domain decomposition [19]. 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. [7]. The transient dispersion behavior of hydrogen and several guidelines for safety in a ventilation model are discussed.
The remaining sections are arranged as follows. Section 2 gives a brief description about the formulas and the balancing domain decomposition method. Section 3 describes the physical model, decomposed mesh, initial and boundary settings, and material properties. Numerical results and discussions are presented in Section 4, and finally Section 5 gives concluding remarks in this research.

Governing Equations and Finite Element Scheme.
Let Ω be the boundary of a three-dimensional polyhedral domain Ω, let 1 (Ω) be the Sobolev space, and let 2 0 (Ω) be the subspace of 2 (Ω) functions with zero mean value. For incompressible, viscous, and laminar flow, the solving of the model can be summarized as finding ( , ) ∈ 1 (Ω) 3 × 2 (Ω) such that for any ∈ (0, ), the following set of equations hold: with the Kronecker delta and the laminar viscosity [kg/ms] of gas mixture.
An initial gas mixture velocity 0 is applied in Ω at = 0. Dirichlet boundary conditions and natural boundary conditions are also applied, where Γ 1 ⊂ Ω; are the outward normal direction components to Ω. The hydrogen concentration ∈ 1 (Ω) is computed by solving where is the diffusion coefficient [m 2 /s] and is the hydrogen source [K/s]. An initial hydrogen concentration 0 is applied in Ω at = 0. Dirichlet and Neumann boundary conditions are set by respectively, where Γ 2 ⊂ Ω, is the outward normal derivate to Ω.
By using a characteristic finite element scheme [20], the material derivative in (1) and (5) at can be written as where 1 (⋅, ⋅) is a position function; is and in (1) and (5), respectively. Let I ℎ ≡ { } denote a triangulation of Ω consisting of tetrahedral elements and let the subscript ℎ denote the representative length of the triangulation. The finite element spaces are as follows: Note that piecewise linear interpolations are employed for velocity and pressure, which do not provide a sufficient condition to connect the velocity and pressure space; therefore, the inf-sup condition [21] should be satisfied. In previous work [22], a penalty Galerkin least-squares (GLS) stabilization method for pressure [23] was employed and it is found difficult to be applied for the simulation of hydrogen dispersion; especially when the flow is very turbulent, the scheme becomes easy to diverge. A new stabilization technique for Boussinesq approximated saddle point problem is Journal of Applied Mathematics 3 employed in this work and the finite element scheme for (1) reads as follows: Let (⋅, ⋅) define the 2 inner product; the continuous bilinear forms and in (10) are introduced by respectively. The stabilization parameter in this work is set as where ℎ denotes the maximum diameter of an element and ‖ ⋅ ‖ ∞ is the maximum norm. Similarly, the finite element scheme of (5) is to find A week coupling of finite element scheme (10) and (14) is applied in this work and the element searching algorithm for Lagrange-Galerkin method only needs to be implemented once in a nonsteady loop [20].

A Balancing Domain Decomposition Method.
After the elimination of inner degree of freedom, the linear system of interface degree of freedom becomes where is the Schur complement, is the interface degrees of freedom, and is the relative RHS vector [16,24]. The Neumann-Neumann preconditioner might be a popular choice for domain decomposition preconditioning [18,25], as the original local Schur complement can be conveniently used as the local operator. However, it shows its drawback of lacking of mechanism to exchange information between subdomains, because of the singularities caused by the floating subdomains. The matrix eventually becomes illconditioned with the increase of the number of subdomains, indicating that it is not an efficient preconditioner for large scale problems. To prevent the propagation of error, Mandel [17] proposed to add a coarse problem to the original Neumann-Neumann preconditioner, which is generated by the attempt to guarantee the solvability of where is the residual corresponding to an error correction Δ . Let be the space of interface degrees of freedom in Ω. Let ( ) be the space of interface degrees of freedom in Ω ( ) after nonoverlapping domain decomposition into subdomains. The local space can be divided into where ( ) is the local coarse space including all the local potential singularities and ( ) 0 is the complement space of ( ) . The global coarse space is constructed by Note here that ( ) is the weighted function to exchange information between subdomains, which is a decomposition of unity on the space and satisfies where ( ) are the 0-1 matrices, mapping from Ω ( ) to Ω. Inspired by previous research on advection and convection problems [26], the local coarse space in this work is constructed by where nodes denotes the total interface nodes; ( ) denotes the restriction operators from current point to the interface of Ω ( ) . The global space of interface degrees of freedom can be decomposed in a similar manner as in (17): where is the coarse projection operator onto . The balancing domain decomposition preconditioned operator proposed by Mandel [17] is then of the form Here, denotes the local solver for the localized and balanced (16). The Neumann-Neumann algorithm is one of the options of the local solver; however, a diagonal scaling  preconditioning performs better in computation, which has been reported by many; see [19,27,28]. The expression of the local level preconditioner is as follows: where is the local Schur complements.

Geometry and Parameters.
The dispersion behaviour of Leaking Hydrogen in a partially open space raises concern in industry. Hydrogen dilutes quickly into a nonflammable concentration when released because of its rapid diffusivity; therefore, it needs to be confined to become a fire hazard; however, as the lightest element in the universe, it is very difficult to be confined. These properties are taken into account when designing hydrogen structures, and these designs help hydrogen escape up and away from the user in case of an unexpected release. In order to assess the risk of an accident caused by a hydrogen leak, a ventilation model with partially open space used by several researchers [1,5,7] is considered in this work, as is shown in Figure 1. During hydrogen dispersion, the discrepancy in concentration is one of the main drive forces of flow motion. According to the Boussinesq approximation and (1), the buoyancy force is  and material parameters are given by Table 1.
As can be seen from (24) and (25), the concentration expansion coefficient keeps approximately independent of with the parameters given in Table 1. This finding makes it very convenient to apply the current solver to simulate the hydrogen dispersion behaviour.

Initial and Boundary Conditions.
To be consistent with the experiments reported by Inoue et al. [7] and the numerical experiments reported by Kanayama et al. [8] and Matsuura et al. [5], 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 end-walls and gradient zero conditions are applied at the roof and door vents. Let Γ inlet , Γ roof , Γ door denote the boundary of inlet, roof, and door, respectively; boundary conditions are set as follows: It can be seen that Rayleigh number in this research is very high (>10 11 ), which is the main reason of oscillations in numerical results reported in [1,8].

Domain Decomposition. The ADVENTURE CAD and ADVENTURE Metis
[29] are used to create the geometry and mesh used in this work. Set as the representative length of the unstructured mesh and set (range, scale) as the refinement function; the local density of mesh around the sensors in Figure 1 is increased by (0.4, 1/3); the mesh density near edges is increased by (0.3, 1/2) due to the boundary layer effect; see Figure 2(a); higher density of meshes appears around the dispersion and ventilation path. The flow fields around the hydrogen inlet, roof vent, and door vent are very complicated and turbulent; therefore, the mesh of these three places is enriched by (0.1, 1/4); especially the dispersion and ventilation path at these three places is enriched by (0.2, 1/5), as is shown by Figure 2(b).
A hierarchical domain decomposition is employed and the model is firstly divided into many parts and the processor element (PE) works only on the part under its charge. Every part is further divided into many subdomains and the domain decomposition is performed by the PE in charge of the part. A Linux cluster of 176 PEs is used in this work, and the decomposition of parts is demonstrated by Figures 2(c) and 2(d).
A mesh sensitive study is done for = 0.1, 0.05, 0.02, and 0.01 in this work and = 0.02 is found to be the best choice for this simulation. In this case, the mesh contains 12,712,960 tetrahedral elements and the global linear system contains 11,075,600 degrees of freedom in total.

Numerical Results and Discussion
The efficiency of new solver is evaluated in the first part of this section, and to validate the scheme, exact solutions, and available benchmark results classical computational models are compared in the second section. The CG convergence is judged by Euclidian norm with a tolerance of 10 −6 ; using an element-based 1 × 2 × 1 norm defined in [22], 10 −4 is set as the criterion for nonsteady iterations.

Efficiency Test.
In order to test the iteration efficiency of the preconditioners, different domain decomposition preconditioners are applied to solve the hydrogen dispersion problem at = 0.02 and Δ = 0.01. The iteration histories of the first nonsteady loop are compared and the result is shown in Figure 3, where "None" denotes no preconditioner is employed, "DIAG" is the diagonal scaling preconditioning, and "BDD" refers to the balancing domain decomposition preconditioner introduced in Section 2.2 by (19), (22), and (23).
It can be seen that balancing domain decomposition preconditioner described in Section 2.2 is more efficient than diagonal scaling preconditioning, and about 10 times of iteration loops are needed for diagonal scaling preconditioning to converge. The iteration does not converge within 5,000 loops if no preconditioning is employed. By using balancing domain preconditioning, the initial value of each iteration is more "correct" (less potential singularities); therefore, it converges faster than other preconditioners.
Another aim of balancing domain decomposition is to prevent the growth of condition numbers with the number of subdomains. The progress is assessed firstly by testing the so-called numerical scalability. By increasing the number of subdomains of a test problem with fixed mesh size (1,000,000 elements), the comparison of the numerical scalability of the balancing domain decomposition in Section 2.2 (BDD) and the diagonal scaling preconditioning (DIAG) are compared together with a domain decomposition with no preconditioner (None).
As is shown in Figure 4, with an increase in the number of subdomains, the number of loops needed for no preconditioning iteration (None) increases dramatically, so are the diagonal scaling iterations (DIAG). The increasing speed slows down when enough subdomains are created, which is due to the limit of test problem size. The iteration time of BDD does not change much for the fixed size test problem when the number of subdomains increases. However, memory usage increases when too many elements exist in one subdomain; therefore, a trade-off strategy is necessary for parameterization. The balancing domain decomposition preconditioner shows good convergence and is more suitable for large scale computations, especially for nonsteady problems with fixed stiffness matrices.   Table 2 for the details of the cluster. The computation took about 24 hours to finish 50,000 nonsteady loops with = 0.02 and Δ = 0.01. The computation results of this work are validated by comparing with the experimental results reported in [7]. positions with experimental values in [7]. The hydrogen concentration is measured in terms of volumetric concentration, which is calculated by In Figure 5, where 0% indicates that the entire volume is occupied by air and 100% represents the opposite, it can be seen that the current numerical results of all 4 sensors agree with the experimental data very well. Due to the acceleration and dilution of air, some oscillations are viewed at sensor 2 and sensor 3 (located at the top of the ventilation model) after 150 s. Compared with numerical results in [1,8,30], 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); the leakage of hydrogen in a partially open space introduces the possibility of accidental ignition, which may cause an explosion in the worst case. In Figure 6, the isosurface of volumetric concentration at 4% is shown, presenting the boundary of flammability inside the partially open space.
Hydrogen leaks at a constant speed with a constant concentration at the inlet, as is shown in (26). 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; after 500 s, the 4% volumetric concentration isosurface does not get lower obviously and the height of the isosurface does not change till the end of the computation. Thus, the reign below the 4% volumetric concentration isosurface is relatively safe. This finding could be very helpful for safety in case of the hydrogen leaks.
Through velocity vectors in Figure 6, several characters of hydrogen dispersion in this model can be observed.  concentration value as most hydrogen escapes up from the roof vent.
(2) Air goes inside the model from the door vent, due to the changes in pressure field; that is to say, hydrogen's concentration near the horizon of the door vent keeps a low level during the dispersion.
(3) The moving air is a driving force that cannot be neglected. It also affects the hydrogen to accumulate inside the partially open space.
Air flows in due to the pressure changes inside the model; to confirm this, the process of pressure change is demonstrated by Figure 7 and pressure isoline at 500 s is given by Figure 8.
It can be seen from Figure 7 that the pressure at the bottom of the model (Sensor 1 and Sensor 4) keeps negative value during the dispersion. Pressure at Sensor 1 and Sensor 2 is less stable than that at Sensor 1 and Sensor 4, indicating that the flow field around the vertical direction of the inlet is very complicated.
Streamline of right view is displayed in Figure 9 to investigate the flow field around the vertical direction of the inlet. As can be seen, the flow becomes turbulent when air goes inside and eddies appear at the corners and within the inlet.

Conclusions
A ventilation model of hydrogen dispersion is numerically simulated in three dimensions using balancing domain decomposition method in this work. By using a new pressure stabilization scheme and large scale computation, the results are more stable and closer than the conventional computation results. The current results are more reliable, and the comparison with experimental data convinces us of the solvability of the current scheme in hydrogen dispersion problems. 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. Turbulent flow appears around the inlet and the door, and recirculation zones are found inside the model. Safe reign inside this partially open model is classified, which could be very helpful for safety in case of the hydrogen leaks.