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

Authors: Qing-He Yao; Xin Pan
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

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.


---

## Body

## 1. 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–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. Section2 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.


## 2. Formulation

### 2.1. Governing Equations and Finite Element Scheme

Let∂Ω be the boundary of a three-dimensional polyhedral domain Ω, let H1(Ω) be the Sobolev space, and let L02(Ω) be the subspace of L2(Ω) functions with zero mean value. For incompressible, viscous, and laminar flow, the solving of the model can be summarized as finding (u,p)∈H1(Ω)3×L2(Ω) such that for any t∈(0,T), the following set of equations hold:
(1)∂tu+(u·∇)u-1ρ∇·σ(u,p)=1ρβ(Cr-C)g,∇·u=0,inΩ×(0,T),
where u is the gas mixture velocity [m/s]; p is the gas mixture pressure [Pa]; ρ is the gas mixture density [kg/m3]; g is the gravity [m/s2]; β is the concentration expansion coefficient determined in Section 3; C is the hydrogen mass concentration [mass%]; Cr is the reference mass concentration; σ(u,p) is the stress tensor [N/m2] defined by
(2)σij(u,p)≡-pδij+2μDij(u),Dij(u)≡12(∂ui∂xj+∂uj∂xi),i,j=1,2,3,
with the Kronecker delta δij and the laminar viscosity μ [kg/ms] of gas mixture.An initial gas mixture velocityu0 is applied in Ω at t=0. Dirichlet boundary conditions
(3)u=u^on⁡Γ1×(0,T)
and natural boundary conditions
(4)∑j=13σijnj=0on⁡∂Ω∖Γ1×(0,T)
are also applied, where Γ1⊂∂Ω; nj are the outward normal direction components to ∂Ω.The hydrogen concentrationC∈H1(Ω) is computed by solving
(5)∂tC+u·∇C-αΔC=SinΩ×[0,T],
where α is the diffusion coefficient [m2/s] and S is the hydrogen source [K/s]. An initial hydrogen concentration C0 is applied in Ω at t=0. Dirichlet and Neumann boundary conditions are set by
(6)C=C^on⁡Γ2×[0,T],(7)α∂nC=0on⁡∂Ω∖Γ2×[0,T],
respectively, where Γ2⊂∂Ω, ∂n is the outward normal derivate to ∂Ω.By using a characteristic finite element scheme [20], the material derivative in (1) and (5) at tn can be written as
(8)∂tM+(u·∇)M≈Mn-Mn-1(X1(un-1,Δt))Δt,
where X1(·,·) is a position function; M is u and C in (1) and (5), respectively.LetIh≡{K} denote a triangulation of Ω consisting of tetrahedral elements and let the subscript h denote the representative length of the triangulation. The finite element spaces are as follows:
(9)Xh≡{vh∈C0(Ω)3;vh|K∈P1(K)3,∀K∈Ih},Mh≡{qh∈C0(Ω);qh|K∈P1(K),∀K∈Ih},Vh(u^)≡{vh∈Xh;vh(P)=u^(P),∀P∈Γ1},Θh(C^)≡{θh∈Mh;θh(P)=C^(P),∀P∈Γ2},Vh≡Vh(0),Qh=Mh∩L02(Ω),Θh≡Θh(0).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 employed in this work and the finite element scheme for (1) reads as follows:
(10)Find{(uhn,phn)∈Vh(g)×Qh}n=1NTsuchthat(uhn-uhn-1(X1(uhn-1,Δt))Δt,vh)+a(uhn,vh)+b(vh,phn)+b(uhn,qh)+∑K∈IhτK(∇phn,-∇qh)K=(β(Cr-Chn)g,vh)+∑K∈IhτK(β(Cr-Chn)g,-∇qh)K.
Let (·,·) define the L2 inner product; the continuous bilinear forms a and b in (10) are introduced by
(11)a(u,v)≡22μρ(D(u),D(v)),(12)b(u,v)≡-(∇·u,q),
respectively. The stabilization parameter in this work is set as
(13)τK=min⁡{Δt2,hK2∥uhn-1∥∞,ρhk224μ},
where hK denotes the maximum diameter of an element K and ∥·∥∞ is the maximum norm.Similarly, the finite element scheme of (5) is to find {Chn∈Θh}n=1NT, such that
(14)(Chn-Chn-1(X1(uhn-1,Δt))Δt,θh)+(a∇Chn,∇θh)=(Sn,θh).
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].


### 2.2. A Balancing Domain Decomposition Method

After the elimination of inner degree of freedom, the linear system of interface degree of freedom becomes(15)SuB=g,
where S is the Schur complement, uB is the interface degrees of freedom, and g 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 ill-conditioned 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
(16)SΔuB=r,
where r is the residual corresponding to an error correction ΔuB.LetV be the space of interface degrees of freedom in Ω. Let V(i) be the space of interface degrees of freedom in Ω(i) after nonoverlapping domain decomposition into N subdomains. The local space can be divided into
(17)V(i)=V0(i)⊕Z(i),
where Z(i) is the local coarse space including all the local potential singularities and V0(i) is the complement space of Z(i). The global coarse space W is constructed by
(18)W={v∈V∣v=∑i=1NN(i)D(i)Z(i)x,x∈V(i)}.
Note here that D(i) is the weighted function to exchange information between subdomains, which is a decomposition of unity on the space V and satisfies
(19)∑i=1NN(i)D(i)N(i)T=I,
where N(i) are the 0-1 matrices, mapping from Ω(i) to Ω. Inspired by previous research on advection and convection problems [26], the local coarse space in this work is constructed by
(20)Z(i)=∑j=1nodesNj(i)I,
where nodes denotes the total interface nodes; Nj(i) denotes the restriction operators from current point to the interface of Ω(i).The global space of interface degrees of freedomV can be decomposed in a similar manner as in (17):
(21)V=∑i=1N(I-P)N(i)D(i)V0(i)⊕W,
where P is the coarse projection operator onto W. The balancing domain decomposition preconditioned operator proposed by Mandel [17] is then of the form
(22)MBDD-1=(P+(I-P)QLS(I-P))S-1.
Here, QL 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:
(23)QL=∑i=1NN(i)(diag⁡(Si))-1N(i)T,
where Si is the local Schur complements.


## 2.1. Governing Equations and Finite Element Scheme

Let∂Ω be the boundary of a three-dimensional polyhedral domain Ω, let H1(Ω) be the Sobolev space, and let L02(Ω) be the subspace of L2(Ω) functions with zero mean value. For incompressible, viscous, and laminar flow, the solving of the model can be summarized as finding (u,p)∈H1(Ω)3×L2(Ω) such that for any t∈(0,T), the following set of equations hold:
(1)∂tu+(u·∇)u-1ρ∇·σ(u,p)=1ρβ(Cr-C)g,∇·u=0,inΩ×(0,T),
where u is the gas mixture velocity [m/s]; p is the gas mixture pressure [Pa]; ρ is the gas mixture density [kg/m3]; g is the gravity [m/s2]; β is the concentration expansion coefficient determined in Section 3; C is the hydrogen mass concentration [mass%]; Cr is the reference mass concentration; σ(u,p) is the stress tensor [N/m2] defined by
(2)σij(u,p)≡-pδij+2μDij(u),Dij(u)≡12(∂ui∂xj+∂uj∂xi),i,j=1,2,3,
with the Kronecker delta δij and the laminar viscosity μ [kg/ms] of gas mixture.An initial gas mixture velocityu0 is applied in Ω at t=0. Dirichlet boundary conditions
(3)u=u^on⁡Γ1×(0,T)
and natural boundary conditions
(4)∑j=13σijnj=0on⁡∂Ω∖Γ1×(0,T)
are also applied, where Γ1⊂∂Ω; nj are the outward normal direction components to ∂Ω.The hydrogen concentrationC∈H1(Ω) is computed by solving
(5)∂tC+u·∇C-αΔC=SinΩ×[0,T],
where α is the diffusion coefficient [m2/s] and S is the hydrogen source [K/s]. An initial hydrogen concentration C0 is applied in Ω at t=0. Dirichlet and Neumann boundary conditions are set by
(6)C=C^on⁡Γ2×[0,T],(7)α∂nC=0on⁡∂Ω∖Γ2×[0,T],
respectively, where Γ2⊂∂Ω, ∂n is the outward normal derivate to ∂Ω.By using a characteristic finite element scheme [20], the material derivative in (1) and (5) at tn can be written as
(8)∂tM+(u·∇)M≈Mn-Mn-1(X1(un-1,Δt))Δt,
where X1(·,·) is a position function; M is u and C in (1) and (5), respectively.LetIh≡{K} denote a triangulation of Ω consisting of tetrahedral elements and let the subscript h denote the representative length of the triangulation. The finite element spaces are as follows:
(9)Xh≡{vh∈C0(Ω)3;vh|K∈P1(K)3,∀K∈Ih},Mh≡{qh∈C0(Ω);qh|K∈P1(K),∀K∈Ih},Vh(u^)≡{vh∈Xh;vh(P)=u^(P),∀P∈Γ1},Θh(C^)≡{θh∈Mh;θh(P)=C^(P),∀P∈Γ2},Vh≡Vh(0),Qh=Mh∩L02(Ω),Θh≡Θh(0).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 employed in this work and the finite element scheme for (1) reads as follows:
(10)Find{(uhn,phn)∈Vh(g)×Qh}n=1NTsuchthat(uhn-uhn-1(X1(uhn-1,Δt))Δt,vh)+a(uhn,vh)+b(vh,phn)+b(uhn,qh)+∑K∈IhτK(∇phn,-∇qh)K=(β(Cr-Chn)g,vh)+∑K∈IhτK(β(Cr-Chn)g,-∇qh)K.
Let (·,·) define the L2 inner product; the continuous bilinear forms a and b in (10) are introduced by
(11)a(u,v)≡22μρ(D(u),D(v)),(12)b(u,v)≡-(∇·u,q),
respectively. The stabilization parameter in this work is set as
(13)τK=min⁡{Δt2,hK2∥uhn-1∥∞,ρhk224μ},
where hK denotes the maximum diameter of an element K and ∥·∥∞ is the maximum norm.Similarly, the finite element scheme of (5) is to find {Chn∈Θh}n=1NT, such that
(14)(Chn-Chn-1(X1(uhn-1,Δt))Δt,θh)+(a∇Chn,∇θh)=(Sn,θh).
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].


## 2.2. A Balancing Domain Decomposition Method

After the elimination of inner degree of freedom, the linear system of interface degree of freedom becomes(15)SuB=g,
where S is the Schur complement, uB is the interface degrees of freedom, and g 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 ill-conditioned 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
(16)SΔuB=r,
where r is the residual corresponding to an error correction ΔuB.LetV be the space of interface degrees of freedom in Ω. Let V(i) be the space of interface degrees of freedom in Ω(i) after nonoverlapping domain decomposition into N subdomains. The local space can be divided into
(17)V(i)=V0(i)⊕Z(i),
where Z(i) is the local coarse space including all the local potential singularities and V0(i) is the complement space of Z(i). The global coarse space W is constructed by
(18)W={v∈V∣v=∑i=1NN(i)D(i)Z(i)x,x∈V(i)}.
Note here that D(i) is the weighted function to exchange information between subdomains, which is a decomposition of unity on the space V and satisfies
(19)∑i=1NN(i)D(i)N(i)T=I,
where N(i) are the 0-1 matrices, mapping from Ω(i) to Ω. Inspired by previous research on advection and convection problems [26], the local coarse space in this work is constructed by
(20)Z(i)=∑j=1nodesNj(i)I,
where nodes denotes the total interface nodes; Nj(i) denotes the restriction operators from current point to the interface of Ω(i).The global space of interface degrees of freedomV can be decomposed in a similar manner as in (17):
(21)V=∑i=1N(I-P)N(i)D(i)V0(i)⊕W,
where P is the coarse projection operator onto W. The balancing domain decomposition preconditioned operator proposed by Mandel [17] is then of the form
(22)MBDD-1=(P+(I-P)QLS(I-P))S-1.
Here, QL 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:
(23)QL=∑i=1NN(i)(diag⁡(Si))-1N(i)T,
where Si is the local Schur complements.


## 3. Modeling

### 3.1. 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. Hydrogen enters from an inlet at a constant speed; ventilation is through a roof vent and a door vent near the opposite end. Four sensors are placed inside the partially open model, and their locations are indicated by red dotted lines.Figure 1
The ventilation model.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
(24)1ρ2(ρ-ρair)g=1ρβ(Cr-C)g,
where ρair is the density of air [kg/m3]; ρ is the mixture gas density [kg/m3] represented by
(25)ρ≡P[CRH2+(1-C)Rair]T,
where P is the absolute pressure of mixture gas [Pa]; T is the absolute temperature [K]; and RH2 and Rair are the gas constants [J/(kg·K)] of hydrogen and air, respectively. The dispersion behaviour is considered under 1 [atm] and 20 [°C] and material parameters are given by Table 1.Table 1
The material properties.
Parameters
Values
P
1.0
1
×
10
5 [m2/s]
T
293 [K]
R
H
2
4,122 [J/(kg·K)]
R
air
287 [J/(kg·K)]
ρ
air
1.209 [kg/m3]
μ
8
.
7
5
×
10
-
6 [kg/ms]
α
6.1
×
10
-
5 [m2/s]
g
(0, −9.8, 0) [m/s2]
S
0 [1/s]As can be seen from (24) and (25), the concentration expansion coefficient β keeps approximately independent of C with the parameters given in Table 1. This finding makes it very convenient to apply the current solver to simulate the hydrogen dispersion behaviour.


### 3.2. 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:
(26)(u1,u2,u3)=(0,0.02,0),C=6.94%,on⁡Γinlet,(u1,u2,u3)=(·,0,0),C=0%,on⁡Γdoor,(u1,u2,u3)=(0,·,0),α∂nC=0,on⁡Γroof,(u1,u2,u3)=(0,0,0),α∂nC=0,on⁡∂Ω∖(Γinlet+Γroof+Γdoor).
Initial gas mixture velocity and hydrogen mass concentration are set as follows:
(27)u0=(0,0,0)C0=0%.It can be seen that Rayleigh number in this research is very high (>1011), which is the main reason of oscillations in numerical results reported in [1, 8].


### 3.3. Domain Decomposition

The ADVENTURE_CAD and ADVENTURE_Metis [29] are used to create the geometry and mesh used in this work. Set D as the representative length of the unstructured mesh and set R(range,scale) as the refinement function; the local density of mesh around the sensors in Figure 1 is increased by R(0.4,1/3); the mesh density near edges is increased by R(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 R(0.1,1/4); especially the dispersion and ventilation path at these three places is enriched by R(0.2,1/5), as is shown by Figure 2(b).Mesh and domain decomposition.
(a)
(b)
(c)
(d)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 Figures2(c) and 2(d).A mesh sensitive study is done forD = 0.1, 0.05, 0.02, and 0.01 in this work and D=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.


## 3.1. 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. Hydrogen enters from an inlet at a constant speed; ventilation is through a roof vent and a door vent near the opposite end. Four sensors are placed inside the partially open model, and their locations are indicated by red dotted lines.Figure 1
The ventilation model.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
(24)1ρ2(ρ-ρair)g=1ρβ(Cr-C)g,
where ρair is the density of air [kg/m3]; ρ is the mixture gas density [kg/m3] represented by
(25)ρ≡P[CRH2+(1-C)Rair]T,
where P is the absolute pressure of mixture gas [Pa]; T is the absolute temperature [K]; and RH2 and Rair are the gas constants [J/(kg·K)] of hydrogen and air, respectively. The dispersion behaviour is considered under 1 [atm] and 20 [°C] and material parameters are given by Table 1.Table 1
The material properties.
Parameters
Values
P
1.0
1
×
10
5 [m2/s]
T
293 [K]
R
H
2
4,122 [J/(kg·K)]
R
air
287 [J/(kg·K)]
ρ
air
1.209 [kg/m3]
μ
8
.
7
5
×
10
-
6 [kg/ms]
α
6.1
×
10
-
5 [m2/s]
g
(0, −9.8, 0) [m/s2]
S
0 [1/s]As can be seen from (24) and (25), the concentration expansion coefficient β keeps approximately independent of C with the parameters given in Table 1. This finding makes it very convenient to apply the current solver to simulate the hydrogen dispersion behaviour.


## 3.2. 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:
(26)(u1,u2,u3)=(0,0.02,0),C=6.94%,on⁡Γinlet,(u1,u2,u3)=(·,0,0),C=0%,on⁡Γdoor,(u1,u2,u3)=(0,·,0),α∂nC=0,on⁡Γroof,(u1,u2,u3)=(0,0,0),α∂nC=0,on⁡∂Ω∖(Γinlet+Γroof+Γdoor).
Initial gas mixture velocity and hydrogen mass concentration are set as follows:
(27)u0=(0,0,0)C0=0%.It can be seen that Rayleigh number in this research is very high (>1011), which is the main reason of oscillations in numerical results reported in [1, 8].


## 3.3. Domain Decomposition

The ADVENTURE_CAD and ADVENTURE_Metis [29] are used to create the geometry and mesh used in this work. Set D as the representative length of the unstructured mesh and set R(range,scale) as the refinement function; the local density of mesh around the sensors in Figure 1 is increased by R(0.4,1/3); the mesh density near edges is increased by R(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 R(0.1,1/4); especially the dispersion and ventilation path at these three places is enriched by R(0.2,1/5), as is shown by Figure 2(b).Mesh and domain decomposition.
(a)
(b)
(c)
(d)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 Figures2(c) and 2(d).A mesh sensitive study is done forD = 0.1, 0.05, 0.02, and 0.01 in this work and D=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.


## 4. 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 of10-6; using an element-based H1×L2×H1 norm defined in [22], 10-4 is set as the criterion for nonsteady iterations.


### 4.1. Efficiency Test

In order to test the iteration efficiency of the preconditioners, different domain decomposition preconditioners are applied to solve the hydrogen dispersion problem atD=0.02 and Δt=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).Figure 3
Convergence comparison.It can be seen that balancing domain decomposition preconditioner described in Section2.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 Section2.2 (BDD) and the diagonal scaling preconditioning (DIAG) are compared together with a domain decomposition with no preconditioner (None).As is shown in Figure4, 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.Figure 4
Numerical scalability.


### 4.2. Validation and Discussion

The balancing domain decomposition preconditioner is employed in this work, A Linux cluster in Sun Yat-sen University is used for this computation, and 176 parts are created for all the PEs. In each part, 50 subdomains are created. See Table2 for the details of the cluster. The computation took about 24 hours to finish 50,000 nonsteady loops with D=0.02 and Δt=0.01.Table 2
Technical details.
Unit
Feature
CPU
Intel Xeon E5606 2.13 GHz (×44)
Memory
12 GB 1333 MHz LV RDIMM (×44)
Hard disk
2 TB 7.2 Krpm SAS (×10)
Network
Infiniband 40 Gb/s
OS
Centos 6
MPI
MVAPICH 1.9The computation results of this work are validated by comparing with the experimental results reported in [7]. Figure 5 compares the computational hydrogen volumetric concentration values at the nodes closest to the four sensor positions with experimental values in [7]. The hydrogen concentration is measured in terms of volumetric concentration, which is calculated by
(28)CV=CRH2CRH2+(1-C)Rair.Figure 5
Numerical results versus experimental results.In Figure5, 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 Figure6, the isosurface of volumetric concentration at 4% is shown, presenting the boundary of flammability inside the partially open space.Figure 6
Vectors and 4% volumetric concentration isosurface.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 Figure6, several characters of hydrogen dispersion in this model can be observed.(1)
Hydrogen flows vertically and then spreads to the upper roof vent after reaching the ceiling; therefore, the sensors near the ceiling (Sensor 2 and Sensor 3) keep a high level of hydrogen concentration. The concentration value of Sensor 1 is slightly higher than that of Sensor 4; both of them keep a low level 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 Figure7 and pressure isoline at 500 s is given by Figure 8.Figure 7
Pressure.Figure 8
Pressure isoline.It can be seen from Figure7 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 Figure9 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.Figure 9
Streamline (right view).


## 4.1. Efficiency Test

In order to test the iteration efficiency of the preconditioners, different domain decomposition preconditioners are applied to solve the hydrogen dispersion problem atD=0.02 and Δt=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).Figure 3
Convergence comparison.It can be seen that balancing domain decomposition preconditioner described in Section2.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 Section2.2 (BDD) and the diagonal scaling preconditioning (DIAG) are compared together with a domain decomposition with no preconditioner (None).As is shown in Figure4, 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.Figure 4
Numerical scalability.


## 4.2. Validation and Discussion

The balancing domain decomposition preconditioner is employed in this work, A Linux cluster in Sun Yat-sen University is used for this computation, and 176 parts are created for all the PEs. In each part, 50 subdomains are created. See Table2 for the details of the cluster. The computation took about 24 hours to finish 50,000 nonsteady loops with D=0.02 and Δt=0.01.Table 2
Technical details.
Unit
Feature
CPU
Intel Xeon E5606 2.13 GHz (×44)
Memory
12 GB 1333 MHz LV RDIMM (×44)
Hard disk
2 TB 7.2 Krpm SAS (×10)
Network
Infiniband 40 Gb/s
OS
Centos 6
MPI
MVAPICH 1.9The computation results of this work are validated by comparing with the experimental results reported in [7]. Figure 5 compares the computational hydrogen volumetric concentration values at the nodes closest to the four sensor positions with experimental values in [7]. The hydrogen concentration is measured in terms of volumetric concentration, which is calculated by
(28)CV=CRH2CRH2+(1-C)Rair.Figure 5
Numerical results versus experimental results.In Figure5, 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 Figure6, the isosurface of volumetric concentration at 4% is shown, presenting the boundary of flammability inside the partially open space.Figure 6
Vectors and 4% volumetric concentration isosurface.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 Figure6, several characters of hydrogen dispersion in this model can be observed.(1)
Hydrogen flows vertically and then spreads to the upper roof vent after reaching the ceiling; therefore, the sensors near the ceiling (Sensor 2 and Sensor 3) keep a high level of hydrogen concentration. The concentration value of Sensor 1 is slightly higher than that of Sensor 4; both of them keep a low level 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 Figure7 and pressure isoline at 500 s is given by Figure 8.Figure 7
Pressure.Figure 8
Pressure isoline.It can be seen from Figure7 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 Figure9 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.Figure 9
Streamline (right view).


## 5. 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.



---

*Source: 686873-2014-07-16.xml*