Modelling Considerations in the Simulation of Hydrogen Dispersion within Tunnel Structures

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 onefifth 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.


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 .

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.89m 2 for Model 1 and 1.75 m 2 for Model 2.
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 2 .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 2 .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.
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: The symbols used are as follows.Ω is a three-dimensional polyhedral domain with boundary ∂Ω, u u 1 , u 2 , u 3 T is the velocity m/s , t is time s , v is the kinematic viscosity coefficient m 2 /s , p is the gas mixture gauge pressure pressure normalized by the density m 2 /s 2 , g g 1 , g 2 , g 3 T is gravity m/s 2 , β is, in this case, the analogous coefficient of buoyancy − , C is the mass concentration of hydrogen − , a is the hydrogen diffusion coefficient in air m 2 /s , S is the source term 1/s , and D ij is the rate of strain tensor 1/s defined by The following boundary conditions are applied, where Γ u and Γ c denote the boundary with specified velocity and concentration, respectively, where T is the total time s , u 0 is the initial velocity m/s , C 0 is the initial concentration − , u is the boundary velocity m/s , C is the boundary concentration − , and σ u, p is the stress tensor normalized by the density m 2 /s 2 defined by with δ ij being the Kronecker delta and n 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 where p is the pressure kg m −1 s −2 of the gas, ρ the density kg m −3 , R the specific gas constant J kg −1 K −1 , and T the temperature K .Taking C as the mass concentration of hydrogen in a hydrogen-air mixture and R H 2 and R air as the specific gas constants for hydrogen and air, respectively, substituting these into 2.7 and rearranging gives Noting that and ρ H 2 p/R H 2 T , ρ air p/R air T , and substituting these into the equation above gives which, when rearranged, give multiplying throughout by g, the acceleration of gravity m s −2 , and equating with the buoyancy force normalized by gas density gives With ρ H 2 0.084 kg m −3 and ρ air 1.21 kg m −3 , 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 −3 yields β 5.77.

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: C 0.0694 − 6.94 mass% .

3.1
The mass concentration C 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 −1 and the inlet surface area of 0.04 m 2 gives a flow rate of 0.0268 m 3 s −1 , 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:

3.2
The other boundaries are effectively non-slip walls, with no inflow or outflow of hydrogen and/or air: The initial conditions are as follows: The numerical parameters are the kinematic viscosity of hydrogen 1.05 × 10 −4 m 2 /s , the diffusion coefficient of hydrogen in air 6.1 × 10 −5 m 2 /s , the analogous coefficient of buoyancy 13.4 − , gravity 0, −9.8, 0 m/s 2 , 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.

Results and Discussion
Model 1 was computed for air inlet velocities ranging from 0.7 to 3 m s −1 .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 −1 , the concentration values increase suddenly and by almost an order of magnitude.From 2 m s −1 and above, the concentration values remain high throughout.Further computations were conducted with wind speeds between 1 to 2 m s −1 , and the change in trend was found to have occurred with a mere 0.005 m s −1 increase in wind speed 1.23 to 1.235 m s −1 .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 −1 are given in Figure 5.
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 −1 or below for this tunnel configuration.reached.Comparing these two, it can be seen that the flow pattern differs significantly: for wind speed 1.23 m s −1 , 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.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 −1 in Figure 7.
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 −1 .Two hydrogen inlet velocities were used-0.1 and 0.67 m s −1 .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 −1 was the value used in 9 .Model 3 with 0.1 m s −1 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 −1 , 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: Slip in the direction of mean flow i.e., the x-direction is implemented by removing the essential Dirichlet boundary condition u 1 0 from the tunnel walls: The second condition is attempted by setting C 0 at Γ t ext : 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 −1 at three settings-base, slip in the direction of mean flow, and slip in the direction of mean flow combined with C 0 at the tunnel exit-at 55 seconds.Figure 9  direction of mean flow and slip in the direction of mean flow combined with C 0 at the tunnel exit at 127 seconds.
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 C 0-slip conditions in the direction of mean flow and the hydrogen inlet velocity at 0.67 m s −1 -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 −1 hydrogen inlet velocity Figure 9 ; both models show the same flow patterns, but higher velocity magnitudes are observed with 0.67 m s −1 , 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.
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.

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 C 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.

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

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

Figures 6 aFigure 6 :
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.

Figure 8 :
Figure 8: Combined velocity vector and isosurface plots for Model 3 at t 55 s.From top to bottom: boundary conditions as per Section 3 base settings , slip in the x-direction, slip in the x-direction and C 0 at the tunnel outlet.The colour bar applies to all three figures.

Figure 9 :
Figure 9: Combined velocity vector and isosurface plots for Model 3 at t 127 s. a Slip in the x-direction, b slip in the x-direction and C 0 at the tunnel outlet.The colour bar applies both figures.

Figure 10 :
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 x-direction, and C 0 at the tunnel exit.