Modelling Considerations in the Simulation of Hydrogen Dispersion within Tunnel Structures

Abstract

The flow of leaked hydrogen gas in tunnel structures is simulated through a free, open source computational fluid dynamics (CFD) code for incompressible thermal convection flow. A one-fifth scale experimental model of a real tunnel is the target model to be simulated. To achieve this, studies on the effects of different boundary conditions, in particular, the wind speed, are carried out on smaller tunnel structures with the same hydrogen inlet boundary conditions. The results suggest a threshold/limiting value of air speed through tunnel. The target model computed with the most suitable boundary conditions shows some agreement with the experimental ones. A method to compute the buoyancy factor used in the code is also presented.

1. Introduction

The development of hydrogen-fuelled vehicles is currently underway. Car garages and road tunnels—the necessary infrastructure—are partially enclosed spaces, where leaked hydrogen might be constrained from dispersing freely into the atmosphere, and accumulate within the structure. Hydrogen-air gas mixtures have a low ignition energy and large flammability range—from 4 to 75% by volume concentration at room temperature—within which lies a range where detonation is possible, generally taken to be at 18–59% by volume concentration. As such, it is important that the behaviour of leaked hydrogen in partially enclosed spaces is investigated before the widespread use of hydrogen by the public is put in place. However, the inherent dangers in carrying out experiments involving hydrogen necessitate safety measures that can be costly and prohibitive. The solution to this is to use computational fluid dynamics (CFD).

The literature on hydrogen dispersion in partially enclosed spaces is extensive. Basic studies on simple geometries with a hydrogen inlet, and two vents have been carried out both experimentally and with CFD (e.g., [1–6]). Guidelines on the use of hydrogen in confined spaces have also been developed [7]. Works on tunnel-like structures, however, are few, for example, [8, 9]. Most studies, whether experimental or computational, involve overpressure calculations and concentrate on the eventual detonation of hydrogen within the tunnel for example, [10, 11].

This work involves the simulation of tunnel structures using a finite-element code developed in house. The code was written to solve incompressible thermal flow problems. The experimental model studied by Sato et al. [9], HT-5, is taken as the target model to be simulated. In this paper, the effects of the parameters—in particular the wind speed through the tunnel—and boundary conditions are investigated. Using the most appropriate parameters and boundary conditions, the target model is simulated and the results compared with the experimental ones.

We describe the tunnel geometries and the basic equations used in this work in the following section. In Section 3, we state the base boundary conditions used. We present our results and a discussion of these in Section 4. The results for boundary conditions other than the base settings are also presented and discussed. Finally, we state our conclusions in the last section.

2. Tunnel Geometries and Basic Equations

In this work, we limit our investigations to a study of tunnel geometry and ventilation (i.e., wind through the tunnel) speeds and their effects on the hydrogen distribution and flow within the tunnel. We explore the computational settings required to achieve converged results and whether these settings adversely affect the flow.

We use three models in our study—two representative models and a one-fifth scale mockup of a typical tunnel for road transport [9]. The mockup model (Model 3) is the only one where some results on the concentration distribution of hydrogen are available, and is thus our target model. Models 1 and 2 are representative tunnel models, used for parametric studies and to make qualitative, if not quantitative, observations about the flow. Through the results of Models 1 and 2, and trial runs of Model 3, the most appropriate settings and boundary conditions are applied to the computation of Model 3.

The representative tunnel models are shown in Figures 1 and 2. The difference in these models is the tunnel cross sectional area—0.89 m² for Model 1 and 1.75 m² for Model 2.

(a)

(b)

Figure 1 

Model 1 and Model 2. The tunnel roof and side walls are highlighted in (a); the hydrogen inlet is highlighted in (b).

(a)

(b)

Figure 2 

Schematic diagrams of Models 1 and 2 (not to scale).

The hydrogen inlet is a block of dimensions 0.2 × 0.2 × 0.1 m  (length × width × height), giving an inlet surface area of 0.04 m². The side from which wind blows is the tunnel entrance; the opposite end is the tunnel exit. The roof cross-section is a semicircle with a diameter equal to the tunnel width. Five sensors are placed downstream from the hydrogen inlet, on the lengthwise axis.

Model 3, shown in Figures 3 and 4, is similar to test HT-5 in [9]. The only difference is the length-test model HT-5 is 78.5 m long, whereas the model here is only 10 m. This difference was due to memory constraints. The tunnel cross sectional area of Model 3 is 3.74 m². The hydrogen inlet is of the same dimensions as Models 1 and 2. In [9], only the flow rate was given for HT-5. We use the same flow rate here.

Figure 3 

Model 3: geometry and mesh.

Figure 4 

Schematic diagram of Model 3.

We use ADVENTURE_sFlow Ver0.5b [12], a CFD code written for incompressible viscous flow and thermal convection problems, to carry out the computations. Developed in-house, it is available online as a free open source code. The code uses the finite-element method to discretize the coupled Navier-Stokes and advection diffusion equations, and the domain decomposition method to enable parallel processing. The code is used for gas dispersion problems through the application of an analogy that relates concentration and temperature. Details regarding the formulations underlying the code, including stabilization methods, can be found in [13]. The code is currently not equipped with any turbulence models. Here, we show the basic equations (the Navier-Stokes and advection-diffusion equations, (2.1)–(2.3)) and initial and boundary conditions used. In addition, we describe methods to compute the buoyancy term and the boundary condition at the hydrogen-air interface, both of which are not described in our previous work: 𝜕𝐮𝜕𝑡+(𝐮⋅𝜕∇𝐶)𝐮−∇2𝜈⋅∇𝐮⋅=𝐷0(𝐮)in+Ω∇𝑝×=(0−,𝑇𝛽𝐶),𝑔inΩ×(0,𝑇),((22..12))𝜕𝑡+𝐮⋅∇𝐶−𝑎Δ𝐶=𝑆inΩ×(0,𝑇).(2.3)

The symbols used are as follows. Ω is a three-dimensional polyhedral domain with boundary 𝜕Ω, 𝐮=(𝑢₁,𝑢₂,𝑢₃)^𝑇 is the velocity [m/s], t is time [s], 𝑣 is the kinematic viscosity coefficient [m²/s], p is the gas mixture gauge pressure (pressure) normalized by the density [m²/s²], 𝑔=(𝑔₁,𝑔₂,𝑔₃)^𝑇 is gravity [m/s²], 𝛽 is, in this case, the analogous coefficient of buoyancy [−], 𝐶 is the mass concentration of hydrogen [−], 𝑎 is the hydrogen diffusion coefficient in air [m²/s], 𝑆 is the source term [1/s], and 𝐷_𝑖𝑗 is the rate of strain tensor [1/s] defined by𝐷_𝑖𝑗(𝐮)≡12𝜕𝜕𝑥𝑢^𝑖_𝑗+𝜕𝑢^𝑗𝜕𝑥_𝑖,𝑖,𝑗=1,2,3.(2.4) The following boundary conditions are applied, where Γ_𝑢 and Γ_𝑐 denote the boundary with specified velocity and concentration, respectively,𝐶𝐮==𝐶𝐮̂oonnΓΓ_𝑢𝑐×(0,𝑇),×(0,𝑇),₃_𝑗=1𝑎𝜕𝜎𝐶_𝑖𝑗𝑛_𝑗=0on𝜕Ω−Γ_𝑢×(0,𝑇),𝐮=𝜕𝐧𝐮₀=0on𝜕Ω−Γ^𝑐×(0,𝑇),,𝐶=𝐶⁰inΩat𝑡=0,(2.5) where 𝑇 is the total time [s], 𝐮⁰ is the initial velocity [m/s], 𝐶⁰ is the initial concentration [−], 𝐮̂ is the boundary velocity [m/s], 𝐶 is the boundary concentration [−], and 𝜎(𝑢,𝑝) is the stress tensor normalized by the density [m²/s²] defined by𝜎_𝑖𝑗=−𝑝𝛿_𝑖𝑗+2𝜈𝐷_𝑖𝑗(𝐮)𝑖,𝑗=1,2,3,(2.6) with 𝛿_𝑖𝑗 being the Kronecker delta and 𝐧 being the unit normal vector.

The computation of the analogous coefficient of buoyancy, 𝛽 in (2.1), is as follows. The ideal gas law may be written as𝑝=𝜌𝑅𝑇,(2.7) where 𝑝 is the pressure [kg m⁻¹ s⁻²] of the gas, 𝜌 the density [kg m⁻³], 𝑅 the specific gas constant [J kg⁻¹ K⁻¹], and 𝑇 the temperature [K]. Taking 𝐶 as the mass concentration of hydrogen in a hydrogen-air mixture and 𝑅_H2 and 𝑅_air as the specific gas constants for hydrogen and air, respectively, substituting these into (2.7) and rearranging gives𝜌=𝑝𝐶𝑅_H2+(1−𝐶)𝑅_air𝑇.(2.8) Noting that and 𝜌_H2=𝑝/𝑅_H2𝑇, 𝜌_air=𝑝/𝑅_air𝑇, and substituting these into the equation above gives 𝜌=1𝐶/𝜌_H2+(1−𝐶)/𝜌_air,(2.9) which, when rearranged, give1−𝜌^air𝜌=𝐶1−𝜌𝜌^air_H2,(2.10) multiplying throughout by 𝑔, the acceleration of gravity [m s⁻²], and equating with the buoyancy force (normalized by gas density) gives1−𝜌𝜌^air𝑔=𝐶1−𝜌𝜌^air_H2𝑔=−𝛽𝐶𝑔,(2.11) from which𝛽=𝜌^air𝜌_H2−1.(2.12) With 𝜌_H2=0.084 kg m⁻³ and 𝜌_air=1.21 kg m⁻³, (2.12) yields 𝛽=13.4. This method works for other gas combinations as well; in the case of helium, for example, using 𝜌_He=0.1785 kg m⁻³ yields 𝛽=5.77.

3. Boundary Conditions and Mesh Details

For the base (reference) settings, we use conditions similar to that used in our simulation of the so-called hallway model [13]. Γ_inlet,Γ_{t_ent}, Γ_{t_ext} and denote the boundary of the hydrogen inlet, the tunnel entrance, and the tunnel exit, respectively. These apply to all three models. At the inlet, hydrogen leaks at a constant rate in the vertical direction. The velocity and the concentration are as follows:𝑢𝑢₁=𝑢₃=0[m/s₂],𝐶==00.0.169o4r[0−.6]7[(m=/6s].94omnaΓss%_inlet),.(3.1)

The mass concentration 𝐶 (=0.0694, or 6.94 by mass%) represents the density ratio between hydrogen and air. The derivation of this value comes from the conservation of hydrogen mass flow [13, 14]. The hydrogen flow velocity of 0.67 m s⁻¹ and the inlet surface area of 0.04 m² gives a flow rate of 0.0268 m³s⁻¹, which is the flow velocity used in test HT-5 of [9]. At the tunnel exit, hydrogen and air are discharged outside freely. At the tunnel entrance, air enters at a constant velocity:³_𝑗=1𝑎𝜕𝐶𝜎_𝑖𝑗𝑛_𝑗=0m²/s²,𝜕𝑛=0[m/s]onΓ^t_ext𝑢₁,𝑢=variable[m/s₂],=𝑢₃=𝐶0=[m0/s[−]].onΓ_{t_ent},(3.2) The other boundaries are effectively non-slip walls, with no inflow or outflow of hydrogen and/or air:𝑢𝑎₁𝜕𝐶=𝑢₂=𝑢₃=0[m/s]𝜕𝑛=0[m/s]on𝜕Ω−Γ^inlet+Γ_{t_ent}+Γ_{t_ext}.(3.3) The initial conditions are as follows:𝑢𝐶₁==𝑢0₂[=−]𝑢₃=0[m/s]inΩ.(3.4) The numerical parameters are the kinematic viscosity of hydrogen (1.05 × 10⁻⁴ [m²/s]), the diffusion coefficient of hydrogen in air (6.1 × 10⁻⁵ [m²/s]), the analogous coefficient of buoyancy (13.4 [−]), gravity (0, −9.8, 0 [m/s²]), and the source term (0 [1/s]).

The same mesh density was used for all three tunnel models. This resulted in Model 1 having 100,432 elements (98,715 degrees of freedom), Model 2 with 206,462 (190,580), and Model 3 with 443,695 (391,285). The timestep used was 1 second; this timestep-mesh size combination was determined based on our experience with previous models. Models 1 and 2 were run to 500 s. Steady-state flow was achieved in most cases.

4. Results and Discussion

Model 1 was computed for air inlet velocities ranging from 0.7 to 3 m s⁻¹. Steady-state flow downstream of the hydrogen inlet was observed for all the computations, indicating (but not guaranteeing) that for these models and their respective boundary conditions and numerical parameters, stable flow was reached and therefore the meshsize-timestep combination used was appropriate.

As the wind speed is increased, the computed concentration values tend to drop as expected, but from 1 to 2 m s⁻¹, the concentration values increase suddenly and by almost an order of magnitude. From 2 m s⁻¹ and above, the concentration values remain high throughout. Further computations were conducted with wind speeds between 1 to 2 m s⁻¹, and the change in trend was found to have occurred with a mere 0.005 m s⁻¹ increase in wind speed (1.23 to 1.235 m s⁻¹). The volumetric concentration profiles at 5 positions measured from the tunnel inlet (Figure 2(a)) of Model 1 for wind speeds 1.23 and 1.235 m s⁻¹ are given in Figure 5.

(a)

(b)

Figure 5 

Model 1: plots of volumetric concentration versus time for wind speed 1.23 and 1.235 m s−1.

Reverse flow of hydrogen gas (i.e., against the flow of air) was not observed in concentration isosurface plots of the computed models. This indicates that the critical velocity, where the wind speed is sufficient to ensure no reverse flow occurs, is 0.7 m s⁻¹ or below for this tunnel configuration.

Figures 6(a)–6(f) show combined velocity vector and concentration isosurface plots for wind speeds 1.23 and 1.235 m s⁻¹ at time = 100 s, when steady-state flow has already been reached. Comparing these two, it can be seen that the flow pattern differs significantly: for wind speed 1.23 m s⁻¹, air tends to flow through the hydrogen plume, whereas for the higher wind speed, air tends to flow around the plume. In other words, hydrogen-air mixing is greater at the slower speed than at the higher one. At higher wind speeds, hydrogen “brushes past” the hydrogen plume, and entrainment of hydrogen by air takes place. The existence of threshold values when using CFD to compute this type of flow has been identified. In addition, and as a result of this, the possibility of hydrogen and air mixing at different rates due to different velocities has also been shown.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

Model 1: combined velocity vector and concentration isosurface plots for wind speeds 1.23 (a, c, e) and 1.235 m s−1 (b, d, f) at 100 s elapsed time. The colour plot applies to all 6 figures. (a, b) Orthogonal view. (c, d) Side view. (e, f) Top (plan) view.

The volumetric concentration profiles at 5 positions measured from the tunnel inlet (Figure 2(b)) of Model 2 are given for air inlet velocities within the range of 0.7 to 0.95 m s⁻¹ in Figure 7.

(a)

(b)

(c)

(d)

Figure 7 

Model 2: volumetric concentration versus time for wind speeds 0.7 to 0.95 m s−1.

Some of the hydrogen concentration values obtained for Models 1 and 2 are very high—up to 0.5 (50%) for Model 1 and 0.11 (11%) for Model 2. However, it should be noted that these models have dimensions and boundary conditions that are unrealistic, and therefore the computed results are not truly representative of the concentrations that might occur in situations involving real tunnels.

The last model (Model 3) was run at a fixed wind speed of 0.26 m s⁻¹. Two hydrogen inlet velocities were used—0.1 and 0.67 m s⁻¹. The value for wind speed was chosen because (i) results of Model 2 hinted that larger cross sectional areas might allow for slower wind speeds, and (ii) 0.26 m s⁻¹ was the value used in [9]. Model 3 with 0.1 m s⁻¹ was used to study the effects of the boundary condition settings on the computation; the value of 0.1 was chosen because previous experience indicated that this was the “optimal” value for trial purposes given the other settings. The other hydrogen inlet value, 0.67 m s⁻¹, resulted in the same volume flow used in [9]. In this work, two boundary condition settings were explored: the relaxing of non-slip conditions in the mean direction of flow and the introduction of an additional boundary condition at the tunnel exit. Wall boundary conditions with no slip are expressed below:𝑎𝑢₁𝜕𝐶=𝑢₂=𝑢₃=0[m/s]𝜕𝑛=0[m/s]on𝜕Ω−Γ^inlet+Γ_tent+Γ_text.(4.1) Slip in the direction of mean flow (i.e., the 𝑥-direction) is implemented by removing the essential (Dirichlet) boundary condition 𝑢₁=0 from the tunnel walls:³_𝑗=1𝑢₂𝜎_1𝑗𝑛_𝑗=0m²/s²𝑎𝜕𝐶=𝑢₃=0[m/s]𝜕𝑛=0[m/s]on𝜕Ω−Γ^inlet+Γ_{t_ent}+Γ_{t_ext}.(4.2) The second condition is attempted by setting 𝐶=0 at Γ_{t_ext}:³_𝑗=1𝐶=𝜎0_𝑖𝑗𝑛m_𝑗=ass0%m²/s²onΓ_{t_ex}.(4.3) This setting was used in the hallway model in [13, 14] to simulate air inflow conditions. Both these conditions are believed to have little or no effect on the mean flow pattern.

Figure 8 compares the flow patterns of Model 3 at a hydrogen inlet speed of 0.1 m s⁻¹ at three settings—base, slip in the direction of mean flow, and slip in the direction of mean flow combined with 𝐶=0 at the tunnel exit—at 55 seconds. Figure 9 compares slip in the direction of mean flow and slip in the direction of mean flow combined with 𝐶=0 at the tunnel exit at 127 seconds.

(a)

(b)

(c)

Figure 8 

Combined velocity vector and isosurface plots for Model 3 at 𝑡=55s. From top to bottom: boundary conditions as per Section 3 (base settings), slip in the 𝑥-direction, slip in the 𝑥-direction and 𝐶=0 at the tunnel outlet. The colour bar applies to all three figures.

(a)

(b)

Figure 9 

Combined velocity vector and isosurface plots for Model 3 at 𝑡=127s. (a) Slip in the 𝑥-direction, (b) slip in the 𝑥-direction and 𝐶=0 at the tunnel outlet. The colour bar applies both figures.

We see from Figures 8 and 9 that the upstream flow pattern is the same. The similarity of the flow patterns for all three cases indicate that the changes made to the boundary conditions have little effect on the flow, as believed.

Figure 10 shows a graph of concentration versus time at 6 sensor positions (see Figure 4) for 𝐶=0—slip conditions in the direction of mean flow and the hydrogen inlet velocity at 0.67 m s⁻¹—this velocity is the same as the target model in [9]. Steady-state flow was not achieved, unlike Models 1 and 2. However, as mentioned previously, Model 3 has been shortened to 10 metres as compared to 78.5 m in [9]; due to this, we only simulate the early stages of the flow. Hydrogen is first detected in order of sensor distance from the hydrogen inlet, with the exception of sensors at 0 m and 1 m, where the sensor at 1 m picks up hydrogen first. This indicates that, unlike Models 1 and 2, reverse flow has occured, probably due to the low wind speed relative to the hydrogen inlet velocity. Reverse flow was also observed in Model 3 with 0.1 m s⁻¹ hydrogen inlet velocity (Figure 9); both models show the same flow patterns, but higher velocity magnitudes are observed with 0.67 m s⁻¹, as expected. The concentration values appear to stabilize at around 0.1 volume concentration. We find that overall, the concentration values obtained for the first 50 seconds are slightly higher than the values suggested in [9]. Longer run times are required to observe the development of concentration profiles within the tunnel; however, this can only be carried out if the tunnel length itself is lengthened or special boundary conditions, as yet undetermined, are applied at the tunnel entrance and exit.

Figure 10 

Hydrogen volumetric concentration versus time for Model 3 (hydrogen inlet velocity of 0.67 m s−1) at 6 sensor positions, having slip in the 𝑥-direction, and 𝐶=0 at the tunnel exit.

Due to the short run time (50 s) and shortened length of the tunnel, the results obtained here cannot be compared directly with the experimental results given in [9]. The concentration values shown in [9] were taken at 3 points in time when the flow through the tunnel was steady, as opposed to the results shown in Figure 10 (where the flow is clearly unsteady). Nonetheless, it is seen that the computed results, which average around 0.1, or 10%, compare somewhat favourably with the experimental results, which are around 0.03 to 0.07 (3–7%) in the first 5 meters downstream of the hydrogen inlet despite the difference in flow regimes.

5. Conclusions

In this work, we have investigated the effect of wind speed and tunnel geometry on the flow and dispersion of hydrogen within tunnel structures. Our work has suggested the existence of “threshold,” or perhaps “limiting” values with regard to CFD modeling of tunnel structures, in this case related to the wind speed. Large variations in entrainment and the rate of hydrogen-air mixing occur when these threshold values (i.e., wind speeds) are crossed. In addition, our findings imply that a change in tunnel cross-sectional area can affect change in the flow pattern within the tunnel, even if all other parameters are maintained. This has implications in the CFD modelling of such structures, as the appropriate computational settings and even boundary conditions have to be reconfirmed.

In this work, we have tested two methods: removal of the non-slip condition in the direction of mean flow, and setting 𝐶=0 at the tunnel outlet, and we found that both enable solutions to converge for longer simulation times without changing the mean flow pattern. More work must be done to verify that these findings are applicable to any tunnel structure that is analyzed with CFD. The possibility that other “threshold” values might exist for other variables at certain mesh densities should also be studied. We plan to investigate these in the future.