CFD benchmark on hydrogen release and dispersion in confined, naturally ventilated space with one vent

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

Available online at www.sciencedirect.com

ScienceDirect journal homepage: www.elsevier.com/locate/he

CFD benchmark on hydrogen release and dispersion in confined, naturally ventilated space with one vent S.G. Giannissi a,*, V. Shentsov b, D. Melideo c, B. Cariteau d, D. Baraldi c, A.G. Venetsanos a, V. Molkov b a

Environmental Research Laboratory, National Center for Scientific Research Demokritos, Aghia Paraskevi, Athens 15310, Greece b HySAFER Centre, University of Ulster, Newtownabbey BT37 0QB, UK c European Commission, Joint Research Centre (JRC), Institute for Energy and Transport (IET), Energy Conversion and Storage Technologies Unit, Westerduinweg 3, 1755 LE Petten, Netherlands d C.E.A. Saclay, D.E.N., D.M.2S., S.T.M.F, Laboratoire d'Instrumentation et Exp
erimentation en m
ecanique des Fluides et Thermohydraulique, 91191 Gif sur Yvette Cedex, France

article info

abstract

Article history:

A CFD benchmark was performed within the HyIndoor project, to study hydrogen release

Received 26 June 2014

and dispersion in a confined space with natural ventilation and one vent. Three experi-

Received in revised form

ments, performed earlier by CEA at their GAMELAN 1 m3 facility, were considered for this

28 November 2014

benchmark. In all three tests helium (instead of hydrogen for safety reasons) was released

Accepted 4 December 2014

vertically upwards at 60 NL/min from a 20 mm orifice near the centre of the enclosure. A

Available online 30 December 2014

different vent size was used for each test. Three HyIndoor partners (JRC, NCSRD and UU) participated in the benchmark, with three different CFD codes, (ANSYS Fluent, ADREA-HF

Keywords:

and ANSYS CFX) and three different turbulence models respectively (transitional SST,

Hydrogen safety

standard k-ε, dynamic Smagorinski LES). In general, good agreement was found between

Enclosure

predicted and measured helium concentrations. However, in the case of the vent with the

Natural ventilation

smallest vertical extension (vent c) all predictions overestimate the concentration at the

CFD

lower part of the enclosure at steady state.

Release and dispersion

Copyright © 2014, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.

Benchmark

Introduction The main objective of the HyIndoor project [1] is to develop guidelines for the safe indoor use of Hydrogen and Fuel Cell (HFC) systems for early markets, e.g. forklift refuelling and operation, back-up power supply, portable power generation.

The project addresses open safety issues that are related to the hydrogen behaviour in indoor releases like the accumulation of hydrogen in confined space, vented deflagrations and under-ventilated jet fires. More specifically, in Work Package 2 (WP2) the aim is to close knowledge gaps in the understanding of dispersion and to develop recommendations to avoid

* Corresponding author. Tel.: þ30 210 6503416. E-mail address: [email protected] (S.G. Giannissi). http://dx.doi.org/10.1016/j.ijhydene.2014.12.013 0360-3199/Copyright © 2014, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.

2416

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

accumulation of hydrogen in confined space by using experimental, analytical and numerical tools. Adding vents to a confined space is a fundamental element in designing safer infrastructures in order to reduce the amount and the concentration of the flammable mass that is formed if an accidental leak should occur. The shape of the vent and the vent areas are 2 key-parameters that significantly affect the gas behaviour inside the enclosure. Computational Fluid Dynamics (CFD) studies are foreseen in the project with several purposes: to provide indications for the set-up of the experiments, to contribute to the analysis of the experimental results, and to extend the range of investigation to parameters that cannot be included in the experimental programme where a limited number of tests can be performed. A crucial step that has to be performed before the use of CFD for the above purposes is to assess the accuracy of the CFD tools by comparison of simulation results with experimental measurements. In this context a series of CFD benchmarking activities was carried out in the project and in this paper the WP2 CFD validation activities are described. Since the CFD validation had to be performed at the beginning of the project when experiments had not been carried out in the project yet, it was necessary to identify experiments in another project. In the framework of the project ANR PANH DIMITRHY several dispersion experiments were performed by CEA. In the experiments helium was release in an 1 m3 enclosure with no vent, one vent and two vents [2,3]. Several injection conditions (e.g nozzle diameter, flow rate) and several vent sizes were investigated. Previous studies related to dispersion of a buoyant jet in enclosures without natural ventilation have been conducted [4e6]. More recently, several experiments regarding buoyant gas dispersion (helium or hydrogen) in a large scale enclosure (80 m3, 40 m3 and 20 m3) [7e9], and in a small scale enclosure (1 m3) [10] without ventilation have been performed. Experiments have also been performed, in order to study the effect of small ventilation of the enclosure [11]. Numerical studies have been presented on the cases in an enclosure without ventilation [12,13], with small forced

ventilation [14] and with natural ventilation [15e20]. Furthermore, in Ref. [3] a simple analytical model based on the homogeneous mixture hypothesis is compared with the experimental results. Among all the above CFD analyses, there are only 2 benchmarks with the participation of several partners and models. In Ref. [12], 10 partners performed simulations trying to reproduce a hydrogen release in a closed 20 m3 vessel. The intercomparison exercise showed a significant amount of scattering for the results for different turbulent models. In Ref. [15], 5 organizations carried out CFD simulations to describe the release of helium in 66 m3 garage-like facility with 2 vents. Also in that benchmark exercise, the scattering of results was relevant e.g. the values of the geometric variance (relative scatter) VG was in a range between 1.09 and 6.4. In this paper, the geometrical enclosure is a relative smaller volume (1 m3) due to the interest of the project for early market applications. For example in a forklift the fuel cell and the hydrogen storage will be enclosed in a small box that has to fit in the limited space that is available in the forklift for that purpose. The objective of the paper is to describe the results of the dispersion CFD benchmarking activities on the GAMELAN facility in the HyIndoor project. The 3 selected experiments for the CFD benchmark were performed by Cariteau et al. [3]. Helium is released in an enclosure with one vent through a 20 mm diameter nozzle and with release rate 60 NL/min. Three different vent sizes were examined. Three HyIndoor partners participated to the benchmark e Joint Research Centre (JRC), National Center for Scientific Research Demokritos (NCSRD) and University of Ulster (UU) e with different CFD codes and modelling approaches. The computational results are compared with the measurements and with each other and relevant conclusions about the accuracy of the CFD tools for hydrogen releases in confined areas and natural ventilation are drawn. Previous studies [19,20] have shown that helium can generally be used as a substitute to hydrogen and exhibit similar behaviour. For hydrogen dispersion in simple enclosures the maximum deviation between helium and hydrogen concentration was 15%. Therefore, any result and conclusion

Fig. 1 e The top view (left) and the side view (right) of the GAMELAN facility. w is for the length and h for the height of the vent respectively, which values for each vent are shown in Table 1. The experimental sensor probes' location is also indicated. The diagrams are not drawn to scale.

JRC

JRC used the ANSYS CFX 14.0 CFD code to simulate the test cases. The code has been validated against several experiments in similar situations. Some of these validations can be found in Refs. [15,24e26]. In Refs. [27,28] the comparison of the code with the results from other codes and/or engineering models is presented. Symmetry along the y axis is assumed and the computational domain was expanded outside the enclosure away from the vent on the top and on the side walls. The mesh is unstructured (tetrahedral) and consists of 150 308, 145 256 and

Table 1 e Experimental conditions. 

Vent Vent area D (mm) QHE T ( C) U dimension (mm2) (NL/min) (m/s) w  h (mm) Vent a Vent b Vent c

900  180 180  180 900  35

162 000 32 400 31 500

20 20 20

60 60 60

25.8 26.4 24.8

3.5 3.5 3.5

Table 2 e CFD parameters of each HyIndoor partner.

JRC

1st order implicit Standard keε Expanded outside Gamelan on the top and on the side walls (1.76  0.665  1.43 m) Yes

In this benchmark three HyIndoor partners participated with three different CFD packages, ANSYS Fluent [21], ANSYS CFX [22] and ADREA-HF [23], and they have simulated the three GAMELAN test cases. In all simulations the continuity equation for the mixture (air þ helium), the momentum conservation equation (NaviereStokes equations) for the mixture velocity and the conservation equation for the mass fraction of helium were solved. UU partner solved also the conservation equation of the mixture energy, whilst JRC and NSCRD partners assumed isothermal conditions. The main characteristics of the modelling approach and of the grid resolution of each partner are summarized in Tables 2 and 3 respectively. Details for each partner are given below.

Transient numerical scheme Turbulence model Computational domain-Total domain dimensions (HxWxL) Symmetry assumption

UU NCSRD

Benchmark description

ADREA-HF QUICK, 1st order upwind (momentum, species)

The experimental set up is a parallelepiped enclosure with a square base of 930 mm width and 1260 mm height (see Fig. 1). Three different vent sizes were examined. All three vents are located centred in the x dimension and 20 mm below the ceiling of the enclosure. Fig. 1 shows the top and side view of the facility. Helium is injected through a 20 mm diameter vertical nozzle 210 mm above the floor in quiescent atmosphere. The flow rate is controlled with mass flow regulators and it is 60 NL/min. Table 1 contains a summary of the experimental conditions for each vent case. 15 experimental sensors (katharometers) were deployed to measure the helium concentration (volume fraction). The sensors' location is depicted in Fig. 1 and it is the same for all vent cases.

ANSYS CFX 14.0 High Resolution scheme (2nd) based on boundedness principles [35] 2nd order SST transitional model (Gamma Theta Model) Expanded outside Gamelan on the top and on the side walls (14.465  8.5  10 m) Yes

Experimental description

CFD code Numerical scheme for convective terms

made in this work are considered to be valid for hydrogen releases too.

ANSYS Fluent 14.5 Bounded central differencing (momentum) 2nd order upwind (species, energy) Bounded 2nd order Dynamic LES Expanded outside Gamelan on the top and on the side walls (2.5  3  2.5 m) No

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

2417

2418

8 CVs (in the half source surface) 864 CVs (in the half source surface) Number of cells on source

UU-LES NCSRD-k-ε

Cartesian ~350 000 CVs (for all cases) (in the symmetric domain)

JRC-SST

Unstructured tetrahedral ~145 000 CVs (for all cases) (in the symmetric domain) Grid type Grid size

Table 3 e The grid characteristics of each HyIndoor partner.

Block-structured hexahedral 2 338 620 CVs (vent a) 1 047 750 CVs (vent b, c) (in the whole source domain) 68 CVs (in the whole source surface)

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

142 256 nodes for the case with vent a, vent b and vent c respectively. Mesh refinement is made in the source, along the jet region and in the vent region. For all cases the number of nodes along the source diameter is 46, and the number of cells on the source symmetric surface is 864. The helium inlet was set inside the pipeline and it was located on the floor (at a distance 210 mm from the pipe exit). 30 nodes were placed along the vertical direction from the floor to the pipe exit. For the turbulence the SST (Shear Stress Transport) transitional model (Gamma Theta Model) was used. In the SST turbulence model, the k-u model is applied in the near wall layers while the k-ε model is used in the free stream flow far from the walls. The transitional model is based on two transport equations, one for the intermittency and one for the transition onset criteria in terms of momentum thickness Reynolds number. The transitional model was extensively validated together with the SST turbulence model for a wide range of transitional flows [29e31]. The choice of a transitional model is justified by the change of regime of the flow after the gas comes out of the pipe. The threshold of laminar flows in channels was identified for a Reynolds number of 2100 [32]. In a pipe it is commonly accepted that the flow is laminar when the Reynolds number is below 1000 and it is fully turbulent when Re > 4000 [33]. However the full understanding of the transition from laminar to turbulent in jets is still missing, as described by Benintendi [34]. Reporting results from other authors, Benintendi [34] states that the onset of the flow instability have been found in jets for a wide range of Reynolds number, starting from 40 to values larger than 1000. Since the Reynolds number in the pipe is about 660, the flow can be considered as laminar inside the pipe and at the exit of the pipe. Nevertheless during the experiment, at a certain distance from the exit the flow undergoes a transition from laminar to turbulent. The flow transition has been also confirmed by the LES simulations that are presented in this work as shown in Figs. 5, 10 and 15. The advection scheme is “high resolution” and it is expressed as: r 4ip ¼ 4up þ bV4$D!

(1)

where V4 is equal to average of the adjacent nodal gradients, r is the vector from 4up is the value at the upwind node, and D! the upwind node to the ip node. In the High Resolution Scheme a special nonlinear recipe is used for the b parameter at each node, and it is based on the boundedness principles used by Barth and Jesperson [35]. Further details on the numerical scheme can be found in ANSYS CFX manual [22]. The transient simulation was with an imposed time step at 0.1 sec. Smaller time steps (0.01, 0.001 sec) were tested with no appreciable differences in the results. For the boundary conditions, no-slip condition on all walls and on the ground was imposed. Symmetry conditions were set on the symmetry plane and its opposite, while all the remaining boundaries are treated as walls.

NCSRD NCSRD used the ADREA-HF CFD code to simulate the test cases. Earlier validation work using ADREA-HF code has been

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

summarized in Ref. [36]. Symmetry along y axis is assumed. The computational domain is expanded outside the enclosure and the boundaries were set away from the vent on the top and on the side walls. The geometry is meshed with a cartesian non-equidistant grid with refinement points on the source, the walls and the vent. The source (half surface) is consists of 8 computational cells and the grid size is 359 550, 391 790 and 303 080 cells for the case with vent a, vent b and vent c respectively. For the turbulence the standard k-ε model with extra buoyancy terms was used. For the discretization of the momentum convective terms and for the species convective terms the Quick and the 1st upwind schemes respectively were used. For time integration 1st order fully implicit scheme was employed and a CFL restriction equal to 10 was imposed. Initial temperature of released helium and air temperature in the domain were the ones that reported in the experimental description. Initial velocities were set to zero in the whole computational domain. Non-slip boundary conditions were applied to the solid surfaces. At the side outlet domain boundaries zero gradient was imposed to all variables except the helium mass fraction for which either a zero gradient boundary condition was applied if outflow occurs or a given value boundary condition (equal to the initial value) if inflow occurs. With this boundary condition is ensured that no helium flow back in the domain will occurred. At the top outlet domain the “constant pressure” boundary condition was prescribed. At the symmetric boundary (y ¼ 0) symmetry boundary conditions were imposed.

UU UU have simulated the test cases with the help of the ANSYS Fluent14.5 CFD software. The model details and methodology of mesh development along with validation against number of release rates and vents can be found in Ref. [33]. The computational domain is expanded outside the enclosure, and no symmetry was assumed. Two block-structured hexahedral computational grids were used: one for case with vent a (grid1) and one (grid2) for case with vent b and vent c. Grid1 consists of 2 338 620 control volumes (CVs) in the whole domain. Grid 2 has much coarser mesh outside of the domain and consists of total 1 047 750 CVs. The helium inflow boundary is a polygon inscribed in a circle with a crosssectional area consisting of 68 CVs. The polygon cross section area is equal to the inlet area. The inflow boundary is located inside the pipeline at a distance 5 cm from the pipe exit. There are 10 cells along the pipe axis from the inflow boundary to the pipe exit for the two grids. For the turbulence modelling dynamic LES was used with constant Sc t¼ 0.7. The Smagorinsky constant Cs is calculated dynamically. Initial value of the turbulent intensity at the inlet was set according to the equation I ¼ 0.16Re1/8 with length scale L ¼ 0.07d. These equations were chosen according to the ANSYS Fluent manual recommendations on the specification of initial turbulence on the inlet boundaries. The values of the initial turbulence obtained by the equation outlined above were validated on the number of previous laminar, transitional and turbulent flow simulations described in Ref. [33]. The transient scheme was the bounded second order

2419

formulation. CFL was kept below 40 at all times with time step 0.01 sec. It is worth noting here that the CFL number choice and therefore the time step in a particular series of simulations was done aiming for a reasonable simulation time that should not exceed a couple weeks on a 48 cores server. To improve the simulation time and accuracy of the results on a specific grid the best practice is to select a time step in a way to keep the CFL number below 40. This value can be smaller or higher subject to the complexity of a problem, the type of a flow, and a model chosen to carry out simulations. Generally a CFL number has its maximum value at the release exit and decreases quickly with a distance from the pipe if an outflow velocity is initially high and drops with the distance, see additional explanation on the CFL choice for example in Ref. [33]. For the spatial discretization bounded central differencing scheme was used in the momentum equation and second order upwind scheme in the species and energy equations. Initial temperature of released helium and air temperature in the domain were set equal to the specified in experimental data. Initial velocities were set to zero in the whole computational domain. Non-slip boundary conditions were applied to all solid surfaces. The “pressure outflow” condition was set at the domain boundaries with the same temperature as in the domain and gauge pressure equal to zero.

Statistical performance measures A useful tool to evaluate the performance of the dispersion models against the experimental dataset is the statistical analysis. The statistical measures often consist of one or more statistical parameter, and/or a graphical presentation. These measures compare the predicted values with the observed (measured) values at all available sensors. There are two types of measures: one type that indicates whether the model in general under- or over-predicts the measurements and one type that indicates the level of scatter. Usually, a pair of the two types is used, in order the statistical analysis to be complete. In this work four different statistical performance indicators were used, which were recommended by Hanna et al. [37,38] for evaluating air dispersion models; fractional bias (FB) and normalized mean square error (NMSE), geometric mean bias (MG) and geometric mean variance (VG). FB is the mean error that defines the residual of the observed (Co) and the predicted concentrations (Cp). The bias is normalized by the data-set averaged concentration. FB is an indicator if the model overall under-over predicts the concentration. With an ideal value equal to 0, positive values mean a model underprediction, while negative values mean an overprediction. FB should be accompanied by the NMSE, the variance counterpart of FB. NMSE indicates the scatter of the entire dataset and estimated the overall deviation between the observed and the predicted values. The ideal value is 0, and the model performs better if the NMSE has smaller values. MG and VG is the second pair of statistical indicators. MG measures the relative mean bias, and by taking the logarithm of the observed to predicted ratios rather than just the ratio, the asymmetry in the averaging process is over come, and MG

2420

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

value is less influenced by extreme ratios. VG measures the relative scatter. Both MG and VG have an ideal value of 1. MG values of 0.5e2.0 can be thought of as “factor of two” overpredictions and underpredictions in the mean, respectively [38]. VG value of about 1.6 indicates a typical factor of two scatter between the individual pairs of measured and predicted values. FB and MG are measures of mean bias and indicate only systematic errors, whilst NMSE and VG reflect both systematic and unsystematic (random) errors. MG and VG are more appropriate for dataset where both observed and predicted values vary by many orders of magnitude. However, MG and VG are strongly influenced by extremely low value, whilst FB and NMSE are strongly influenced by infrequently occurring high observed and predicted values. The equations for each measure are the followings. The overbar denotes the average over the entire dataset. FB ¼ 2

Co  Cp

(2)

Co þ Cp


NMSE ¼

Co  Cp

2

Co $Cp

3   Co 5 4 MG ¼ exp ln Cp

(3)

2

(4)

2

3  2 Co 5 4 VG ¼ exp ln Cp

(5)

Model acceptance criteria for the above statistical indicators were proposed by Chang et al. [39]. A “good” model would be expected to have mean bias ±30% of the mean, i.e. jFBj < 0:3 or 0.7 < MG < 1.3, and random scatter about a factor of two to three, i.e. NMSE < 1.5 or VG < 4. Besides the statistical parameters graphical presentation can also be used to evaluate the dispersion models. In the present work the scatter plot is used. Scatter plot is a diagram with the observed values versus the predicted values of concentration at each sensor. The points of the perfect model would be along the line y ¼ x. The closer to this line the points are the better the model's performance is.

the onset of flow instability and the transition from laminar to turbulent flow in free jets are still an open issue. The issue becomes even more complicate in the configuration that is under investigation in this paper: a jet in an enclosure with a vent. In that specific configuration, additional causes of flow disturbances are created by the interactions between the jet and the ceiling, by the interactions between the flow generated by the jet itself and the walls, and by the flows which go through the vent between the inside and outside of the enclosure. All those disturbances can trigger the onset of flow instability and cause an earlier transition compared to the configuration of a free jet. Since the flow regime affects significantly the transport of helium in the enclosure, the larger predicted concentration near to the ceiling and some inaccuracies in the prediction in the unsteady stages with the SST transitional model could be due to a slightly inaccurate prediction of the transition time and position. Moreover, high differences are observed in the number of cells covering the source surface and in the grid size, especially between the RANS (SST, k-ε) simulations and the LES simulation, whereas the results at steady state are quite close. Therefore, an ideal cell/CV/node number on the source surface and in the whole domain considering the computational cost should be investigated further for each modelling approach. The oscillations observed by the UU prediction are due to the LES model that was used for the turbulence modelling. The average value of the LES results seems to be close to the kε model prediction. The reason for the slight underestimation and delay in dynamics of the LES model for the lower sensor was mentioned earlier in Ref. [33] and is valid here. The decrease of the time step would probably improve the LES model simulations and decrease to some extent the underprediction of helium concentrations at lower sensors, but will result in longer simulation time.

Results and discussion In this section the results by all participated partners for each vent size separately are presented and discussed.

Vent a (900  180 mm) Vent a is the largest vent studied here with surface area 1620 cm2. Fig. 2 shows the predicted time histories and the measured ones at the sensor mast M4. All predictions are consistent with the experiment. However, both k-ε model and LES results, show better agreement than SST transitional model at the higher sensor, at which the SST transitional overpredicts the concentration levels. As explained in Section JRC,

Fig. 2 e Vent a. Comparison of the predicted helium concentration time histories and the experimental ones at sensor mast M4.

2421

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

In Fig. 3 the vertical concentration profiles at steady state and for M4 and M2 sensors are depicted. The measurements during steady state exhibit small amplitude perturbations, as discussed in Ref. [3], and a time averaged is taken as steady state value. The time average value that was taken in both experiment and predictions is the arithmetic mean of the values of fifty seconds around the steady state. For the case with vent a it is the arithmetic mean of the values at 350e400 sec. For now long the time averaged steady state values will be referred as the values at steady state. All concentration predictions are consistent with the experimental data. However, the k-ε model and the LES model overpredict the concentration at the lower sensors and underpredict it at the sensor below the vent (approximately at 1000 mm altitude). The overprediction of the bottom concentrations with the k-ε model can be attributed to the fact that k-ε model produces higher turbulent diffusivities in the region near the source where the flow is laminar. SST transitional model prediction is in very good agreement with experiment at the bottom part of the enclosure, whilst it overpredicts the concentration at the top part of the enclosure. In general, k-ε model and LES results are similar, and they are better than SST transitional model results at the 2 higher sensors, while SST transitional model gives better prediction at the lower sensors. As explained above, the overprediction of the SST transitional model in the upper region of the facility could be due to a small inaccuracy in the prediction of the time and position of the laminar turbulent transition. The homogeneous layer at the bottom of the enclosure below the vent is captured by all partners, as well as the stratification at the top. Table 4 displays the statistical measures for the concentration at steady state overall at 14 sensors for all partners. The values are very close to the ideal values for all partners. According to MG all three modelling approaches overall overpredict the concentration. Similarly, according to the FB, SST transitional model and LES model predictions

Fig. 3 e Vent a. Numerical versus experimental helium concentrations (by vol.) at M4 and M2 sensor masts and at steady state.

Table 4 e Statistical performance measures for case with vent a.

FB NMSE MG VG

Ideal value

SST (JRC)

k-ε (NCSRD)

LES (UU)

0 0 1 1

0.09511 0.289587 0.914907 1.011784

0.001646 0.489877 0.975351 1.038222

0.03421 0.167177 0.923228 1.018022

overestimate the concentration. For the k-ε model predictions FB shows underprediction, contrast to the MG indication. However, as it was mentioned in Statistical performance measure FB is strongly influenced by infrequently high under/or overpredictions. Therefore, the high underestimation of the k-ε model prediction at 940 mm (see Fig. 3) is probably responsible for the FB indication of underprediction. MG indication of overall overprediction seems to be more trustworthy in this case. Finally, k-ε model results have a greater scatter than SST transitional model and LES model results according to the NMSE and VG indicators. Fig. 4 is the scatter plot of the observed versus the predicted concentration at steady state for all partners and for 14 sensors. Different colours (in the web version) and markers in the points are for the different partner. The symbols with background represent the bottom sensors (<820 mm) and the symbols without background the top sensors (820 mm). The majority of the points is in the vicinity of the diagonal line. The diagonal line corresponds to a perfect model. The point with coordinates (5.2, 2.9) that deviates more from the diagonal line belongs to the k-ε model prediction and corresponds to the sensor which is located 940 mm above the floor. This point is located at the beginning of the steep concentration gradient in the vertical direction and its high deviation from measurements is also obvious in Fig. 3. The points that

Fig. 4 e Vent a. Scatter plot of the observed versus the predicted concentration for all models and for 14 sensors for vent a. The symbols with background represent the bottom sensors (<820 mm) and the symbols without background the top sensors (≥820 mm).

2422

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

Fig. 5 e Vent a. The predicted concentration contours on the symmetry plane at 400 sec: SST transitional model (left), k-ε model (center) and LES model (right) predictions.

lay above the diagonal line show a model overprediction, while the ones that lay below a model underprediction. So, the scatter plot indicates that the SST transitional model tends to overestimate the concentration at all sensors, while both k-ε model and LES model overestimate the concentration at the bottom sensors, and underestimate it at most of the top sensors. In Fig. 5 the predicted concentration contours by the SST transitional model, the k-ε model and the LES model on the symmetry plane and at 400 sec are illustrated. The predicted contours are similar for the SST transitional model and the k-ε model. However, in the SST transitional model contours it seems that more fresh air enters the enclosure, and dilutes the mixture more at lower heights close to the wall. This behaviour may be attributed to the fact that: the SST transitional model switches between the k-u model near walls and standard k-ε model far from the walls, while the k-ε model uses the standard k-ε model in the whole domain. The LES model predicted contours have similar behaviour with the other models, but the mixture seems less diluted on the top. Moreover, in the LES model contours image the produced by the turbulence eddies are distinct. This is happening because the LES turbulence model solves the governing equations for the instantaneous variables and not for the time average values like the RANS models. The formation of homogeneous mixture on the bottom of the enclosure is obvious in all three images. Fig. 6 shows the time histories of the experimental and predicted total helium mass inside the enclosure. The total mass is calculated based on the helium volume fraction of the sensors. The enclosure is divided in rectangular parallelepipeds. Each parallelepiped corresponds to a unique sensor height. It is assumed that the concentration in each parallelepiped is uniform and equal to the mean value of the enclosed sensors. According to Fig. 6 k-ε model prediction exhibits very good agreement with the experiment at steady state from 300 to 400 s, overestimating at the beginning, while LES model shows

good agreement at the beginning and slightly overestimating at the steady state. The SST transitional model predictions tend to overestimation throughout the whole period of the simulation, due to the small inaccuracy in the prediction of the time and position of the laminar turbulent transition.

Vent b (180  180 mm) In Fig. 7 the predicted and measured helium concentration time histories at M4 sensors are presented. The results by all partners are generally in good agreement with experiment for the majority of the sensors. However, the concentration level at the lowest sensor (100 mm altitude) is over-predicted by both the SST transitional model and the k-ε model at the early stage of the release, while LES model prediction is closer to the

Fig. 6 e Vent a. The total experimental and predicted helium mass inside the enclosure. The dashed line indicates the total released mass.

2423

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

Fig. 7 e Vent b. Comparison of the predicted helium concentration time histories and the experimental ones at sensor mast M4.

measurement. At steady state the predictions by all partners are in good agreement with the experiment. The SST transitional model tends to overestimate the concentration at the top. The k-ε model gives underprediction for the majority of the sensors. The LES model prediction is in better agreement with the experiment for the lower sensor, but overall it underestimates the concentration. The reason for the underestimation, as it was explained earlier could be the choice of the time step. Smaller time step could result in less underprediction, but also in higher computational time. Fig. 8 shows the verticalconcentration profiles at steady state for M4 and M2 sensors. As mentioned in Section Vent a (900  180 mm) the steady state value is the arithmetic mean of the values of fifty seconds around the steady state (for

vent b it is at 950e1000 sec). All three partners' predictions are consistent with the experiment. It is observed in both experiment and prediction that two homogeneous layers with different concentration levels are formed; one at the bottom and one with higher concentration levels at the top. The layers are separated by a high gradient. The SST transitional model prediction tends to overestimate the concentration at the lower sensors and to overestimate it at the top sensors. LES model results give underprediction at all sensors, whilst k-ε model results give underprediction at the lower sensors and overprediction at the top sensors. In general, k-ε model and LES model exhibit better agreement at the higher sensors than the SST transitional model. Table 5 shows the statistical measures for the concentration at steady state overall at 14 sensors for all partners. The values of the performance indicators are very close to the ideal value for all partners. According to FB and MG, the SST transitional model overall slightly overestimates the concentration at steady state, while k-ε model and LES model slightly underestimate it. LES model results seem to have the greatest scatter according the NMSE. The scatter plot in Fig. 9 reveals the overprediction of the SST transitional model at the majority of the sensors (points above the diagonal line), and the underprediction of both k-ε model and LES model at most of the sensors and at all sensors respectively. However, in general all predictions are in good agreement with the experiment. Fig. 10 shows the predicted by all partners concentration contours on the symmetry plane and at t ¼ 1000 sec. It is observed that the mixture is more diluted in the area close to the wall along the air entrainment zone in the k-ε model prediction. Other than that all contours are similar. Fig. 11 displays the time histories of the total experimental and predicted helium mass inside the enclosure. The SST transitional model prediction is in very good agreement with the experiment. The LES model underestimates the total mass from the early stage of the release, while the k-ε model underestimates it only at steady state.

Vent c (900  35 mm) Fig. 12 shows the predicted time histories and the measured ones at the sensor mast M4 and for test case with vent c. The agreement with experiments is acceptable in the upper part of the facility but in all simulations a poor agreement with the experiment is found at the lowest sensor located below the injection point (red (in the web version) line): the helium concentration levels are significantly overestimated. At that sensor and before the steady state is reached k-ε model predictions are close, while SST transitional model and LES model Table 5 e Statistical performance measures for case with vent b.

Fig. 8 e Vent b. Numerical versus experimental helium concentrations (by vol.) at M4 and M2 sensor masts and at steady state.

FB NMSE MG VG

Ideal value

SST (JRC)

k-ε (NCSRD)

LES (UU)

0 0 1 1

0.0181 0.769074 0.982843 1.001964

0.024635 0.573088 1.028072 1.001985

0.044486 0.893381 1.04755 1.003136

2424

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

Fig. 9 e Vent b. Scatter plot of the observed versus the predicted concentration for all partners and for 14 sensors for vent b. The symbols with background represent the bottom sensors (<820 mm) and the symbols without background the top sensors (≥820 mm).

are similar and give better prediction, especially at initial times (until ~600 sec). However, the steady state is reached at the same concentration level for both SST transitional model and k-ε model, while for LES model the steady state is reached at lower concentration level, closer to the experiment, but still highly overestimated. At the higher sensors SST transitional model and k-ε model results are similar and close to the measurements, while the LES model predicted concentration is underestimated for the majority of sensors. Fig. 13 depicts the vertical concentration profiles at steady state for the SST transitional model, the k-ε model and the LES model. Steady state refers to the arithmetic mean value of the

1120e1170 sec. It is shown that both SST transitional model and k-ε model predictions are consistent with the experiment for the large majority of sensors, while there is a significant over-estimation at the 2 bottom sensors. LES model prediction underestimates the concentration at higher sensors, and it exhibits better agreement with the experiment at the 2 bottom sensors compared to the other two models. Still, the measured concentration levels at the bottom are less than all predicted ones. The nearly homogeneous mixture on the top of the enclosure has been well predicted by all three partners. In Table 6 the values of statistical performance indicators are illustrated for all partners. The values have been calculated at steady state and for overall 14 sensors. According to the indicators the SST transitional model and the k-ε model overall overpredict the concentration, whilst the LES model underpredicts it. NMSE indicates great scatter by all models, with LES model to have the smallest one among them. This great value of scatter is due to the high overprediction of the concentration at the lower sensors (at 100 mm and 220 mm above the floor). These remarks are confirmed by the scatter plot in Fig. 14. Although the majority of the sensors is close to the diagonal line, there are two points of SST transitional model and k-ε model results (at 100 mm and 220 mm) and one point of LES model results (at 100 mm) that lay away from the line. The fact that serious overprediction at the lower part of the enclosure is observed in the case of vent c (the smallest vent) and at early stage of the release in the case of vent b, though the release conditions were the same in all vent cases, led us to the conclusion that the vent size is of great significance, and influences the prediction. The driving force for the outside air to enter the enclosure is the density difference between the enclosure and the exterior as helium fills in the enclosure. The bigger vent (vent a) allows more air to entrain, resulting in a more diluted mixture, as it can be concluded by the maximum concentration levels at steady state. As time progresses the dilution of the mixture leads to less incoming air until steady state is

Fig. 10 e Vent b. The predicted concentration contours on the symmetry plane at 1000 sec: SST transitional model (left), k-ε model (center) and LES model (right) predictions.

2425

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

Fig. 11 e Vent b. The total experimental and predicted helium mass inside the enclosure. The dashed line indicates the total released mass.

reached. Eventually, the stabilizing effect of buoyancy due to density gradient is more effective than the incoming fresh air kinetic energy, and therefore, homogenization is only possible in the lower part, as concluded in Ref. [3]. In the cases with vent b and c the volume flow rate of fresh air is less than in the case of vent a, since the surface of the vents is smaller. Although the total area of vent b is similar to the total area of vent c (32 400 mm2 and 31 500 mm2 respectively), they differ in height, which it is known that influences the volume rate [3,40]; the bigger the vertical extension of the vent the greater the volume rate is. Therefore, even though in both vent cases the helium accumulates in the upper part, in the case with vent b the incoming air kinetic energy is sufficient to produce a homogeneous layer at the bottom too. In the case with vent c

Fig. 13 e Vent c. Numerical versus experimental helium concentrations (by vol.) at M4 and M2 sensor masts and at steady state.

the dilution of the mixture is less, and the mixture stays unperturbed in almost the whole domain, similar to an entirely close box. The effect of that behaviour is not captured in the lower part of the facility by the CFD calculations. The exact reason for that is not obvious yet. The overestimated turbulent diffusivity could be a reason, since there is difference between the RANS models and the LES model results with the latter to show better agreement with the experiment in the lower region. Higher predicted values of turbulent diffusivity mean higher levels of diffusion. Due to the turbulent diffusion the hydrogen-rich mixture in the upper part of the enclosure flows down to the lower part drifted by the incoming fresh air. As a result high hydrogen concentration is predicted in the lower part of the enclosure and an almost homogeneous mixture is formed in the entire facility. Fig. 15 depicts the predicted concentration contours on the symmetry plane and at 1100 sec. The images are similar between the SST transitional model and k-ε model, while the LES model predicts a more dilute mixture on the top, as was showed by the previous diagrams, too. Moreover, LES model results seem more stratified on the bottom than the SST transitional model and k-ε model predictions that reproduce a nearly homogeneous layer on the bottom. As shown in Fig. 16, the k-ε model overpredicts the helium total mass inside the enclosure with a relative error 14% at steady state. At the early stage of the release the SST transitional model is consistent with the experiment; however, it

Table 6 e Statistical performance measures for case with vent c.

Fig. 12 e Vent c. Comparison of the predicted helium concentration time histories and the experimental ones at sensor mast M4.

FB NMSE MG VG

Ideal value

SST (JRC)

k-ε (NCSRD)

LES (UU)

0 0 1 1

0.07554 30.56763 0.90375 1.064903

0.09556 34.02919 0.885955 1.068611

0.057323 17.04979 1.036695 1.040506

2426

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

Fig. 14 e Vent c. Scatter plot of the observed versus the predicted concentration for all partners and for 14 sensors for vent c. The symbols with background represent the bottom sensors (<820 mm) and the symbols without background the top sensors (≥820 mm).

overestimates it at steady state, similar to the k-ε model prediction. The LES model underpredicts the total mass until almost 1100 sec, and then it is in very good agreement with the experiment.

Comparison between cases Both experiment and prediction showed that in case with large vent surface the air inflow rate is higher, resulting in dilution of the mixture and less helium accumulation in the enclosure. However, the vent surface is not the only factor

Fig. 16 e Vent c. The total experimental and predicted helium mass inside the enclosure. The dashed line indicates the total released mass.

that influences the air inflow rate. The vertical extension of the vent is also significant. Vent b and vent c have similar surfaces, yet they exhibit different behaviour due to the different vent shape/heights. With the thinnest vent (vent c), higher concentration levels are observed on the top part of the enclosure. According to the experiments for the vent with the highest vertical extension (vent a) a homogeneous layer at the bottom is formed quickly, while stratified layer is formed on the top of the enclosure. For the vent with small surface but the same vertical extension as vent a (vent b) there is a homogenization in the lower part and the upper part of the enclosure. The homogeneous upper layer is formed at the early stage of the release driven by buoyancy. As fresh air

Fig. 15 e Vent c. The predicted concentration contours on the symmetry plane at 1100 sec: SST transitional model (left), k-ε model (center) and LES model (down) predictions.

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

enters the enclosure and flows towards the lower part of the enclosure, it homogenizes the near bottom area. Stratification is formed in the middle of the enclosure due to high density gradient. For the thinnest vent, the air that enters the enclosure is less, and permits the formation of stratification during the injection. As far as the modelling approaches are concerned, both RANS models (SST transitional and k-ε) are in very good agreement with the experiment for the vents with large vertical extension (vent a, b), while for the thinnest vent they overpredict the concentration on the lower part of the enclosure. Comparison between the RANS models shows that at steady state the SST transitional model overestimates the concentration at most of the sensors in vent a and b, while k-ε model underestimates it. In addition, SST transitional model captures better the stratification area between the upper and lower homogeneous layer in vent a and b. LES model exhibits good agreement with experiment in the vents with large vertical extension. For the thinnest vent it improves the prediction on the bottom layer compared to the RANS models, but there is still overprediction. However, it tends to underestimate the concentration on the top at latter times.

Conclusions In the framework of the HyIndoor project a CFD validation benchmark was performed to assess the accuracy of the CFD models that are used in the project to provide useful indications on the effect of some key parameters (size and shape of the vent) for the development of safety guidelines. Selected GAMELAN experiments were simulated with different numerical models and the numerical results were compared to each other and with the experiment. Helium was released to an almost cubic enclosure with one opening. Three different vent sizes were tested. The simulations produced overall good results, especially for the larger vent, i.e. vent a (900  180 mm). The partners that participated in this benchmark are JRC with the ANSYS CFX 14.0 CFD software (SST transitional model), NCSRD with the ADREA-HF CFD code (standard k-ε model) and UU with ANSYS Fluent 14.5 CFD software (LES model). The results in case with vent a (900  180 mm) are in good agreement with the experiment. In the case with vent b (180  180 mm) the helium concentration is overestimated by SST transitional model and k-ε model at the sensor below the injection point and before the steady state is reached; whilst it is slightly under-predicted by the LES model at the early stage of the release. However, the prediction is improved with time for all partners and the steady state is well predicted at most of the sensors. Any measured homogeneous layer is predicted satisfactorily in both test cases and by all models. The stratification region is captured better by the LES model in vent a (900  180 mm), and by the transitional SST transitional model in vent b (180  180 mm), whilst the k-ε model underpredicts the concentration in this region in both cases. In the case with vent c (900  35 mm) SST transitional model and k-ε model predictions are similar and consistent with the experiment except at the sensor below the injection

2427

point. In that sensor the highest deviation from measurement is found by all simulations. In that sensor LES model results show better agreement with the experiment, but still the concentration at steady state is over-predicted. The formation of homogeneous layer on the top of the enclosure and at steady state is well reproduced by all models. However, the LES model predicts a more dilute mixture on the top of the enclosure and a weakly stratified layer on the bottom, while SST transitional model and k-ε model predict a nearly homogeneous layer on the bottom. The overprediction of the helium concentration at the bottom in the case with the smaller vents (vent b and c), and especially, in the case with vent c, can be attributed to the overestimated turbulent diffusivity, which leads to more diffused results. For those cases, the LES model (UU) produces better results than the RANS models (JRC, NCSRD) in the lower section of the facility. In general, it can be concluded that all simulations are in good agreement with the experiment for all vent cases at the majority of the sensors. However, in the cases with smaller vents, and especially, the case with the very small vertical extension (vent c e 900  35 mm), the simulations produce the least good results with high overprediction on the bottom of the enclosure. This behaviour can be attributed to the overestimated turbulent diffusivity in the predictions. The statistical performance evaluation of the different CFD approaches revealed an overall satisfactory agreement between experiments and simulations. The value of the relative scatter VG is included in a range that goes from 1 to 1.06 and that is a very good result. A factor of 2 scatter would imply a VG ~ 1.6 [39]. The statistical analysis shows some tendencies like an overall over-prediction for all models for the case with the larger vent (vent a e 900  180 mm). In the case with the smallest vent (vent c e 900  35 mm) the SST transitional model and the k-ε model overall overpredict the concentration, while the LES model underpredicts it. In this case the statistical indicators showed the largest scatter among all vent cases. In the small vent case but with the higher vertical extension (vent b e 180  180 mm), SST transitional model overpredicts the concentration, while k-ε model and LES model under-predict it. The overall good agreement between simulations and experiments provided positive indications that the predictive capabilities of the CFD models are accurate enough to justify the use of CFD codes to investigate the effect of some relevant parameters like vent shape and size in facilities whose geometrical features are equal or similar to the those of the GAMELAN facility in terms of size and shape. Therefore the use of validated CFD models can contribute to the project main objective, providing useful indications for the development of guidelines for the safe use of hydrogen in indoor configurations.

Acknowledgements The work was performed within the framework of the project HyIndoor (“Pre-normative research on safe indoor use of fuel cells and hydrogen systems”, Grant agreement no: 278534),

2428

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

funded by the Fuel Cells and Hydrogen Joint Undertaking, a public-private partnership between the European Commission, fuel cell and hydrogen industries represented by the NEW Industry Grouping and the research community represented by the Research Grouping N.ERGHY, which the authors gratefully acknowledge. The simulated experiments have been performed by CEA within the framework of the project ANR PANH DIMITRHY, which the authors would also gratefully acknowledge.

references

[1] HyIndoor. Pre normative research on the indoor use of fuel cells and hydrogen systems. Project website: http://www. hyindoor.eu/. [2] Bernard-Michel G, Gentilhomme O, Trochon J, Tviazmina E. Rapport DEN/DANS/DM2S/STMF/LIEFT RT/2011/A. CEA; 2011. [3] Cariteau B, Tkatschenko I. Experimental study of the effects of vent geometry on the dispersion of a buoyant gas in a small enclosure. Int J Hydrogen Energy 2013;38(19):8030e8. [4] Baines WD, Turner JS. Turbulent buoyant convection from a source in confined region. J Fluid Mech 1969;37:51e80. [5] Germeles AE. Forced plumes and mixing of liquids in tanks. J Fluid Mech 1969;71:601e23. [6] Cleaver RP, Marshall MR, Linden PF. The build-up of concentration within a single enclosed volume following a release of natural gas. J Hazard Mater 1994;36:209e26. [7] Lacome JM, Jamois D, Perrette L, Proust CH. Large-scale hydrogen release in an isothermal confined area. Int J Hydrogen Energy 2011;36:2302e12. [8] Gupta S, Brinster J, Studer E, Tkatschenko I. Hydrogen related risks within a private garage: concentration measurements in a realistic full scale experimental facility. Int J Hydrogen Energy 2009;34:5902e11. [9] Shebeko YN, Keller VD, Yeremenko OY, Smolin IM, Serkin MA, Korolchenko AY. Regularities of formation and combustion of local hydrogeneair mixtures in a large volume. Chem Ind 1988;21e24:728e31. [10] Cariteau B, Tkatschenko I. Experimental study of the concentration build-up regimes in an enclosure without ventilation. Int J Hydrogen Energy 2012;37:17400e8. [11] Cariteau B, Brinster J, Tkatschenko I. Experiments on the distribution of concentration due to buoyant gas low flow rate release in an enclosure. Int J Hydrogen Energy 2011;36:2505e12. ~ as JM, Crespo A, [12] Gallego E, Migoya E, Martı´n-Valdepen Garcı´a J, Venetsanos A, et al. An intercomparison exercise on the capabilities of CFD models to predict distribution and mixing of H2 in a closed vessel. Int J Hydrogen Energy 2007;32:2235e45. [13] Venetsanos AG, Adams P, Azkarate I, Bengaouer A, Brett L, Carcassi MN, et al. On the use of hydrogen in confined spaces: results from the internal project InsHyde. Int J Hydrogen Energy 2011;36:2693e9. [14] Venetsanos AG, Papanikolaou E, Cariteau B, Adams P, Bengaouer A. Hydrogen permeation from CGH2 vehicles in garages: CFD dispersion calculations and experimental validation. Int J Hydrogen Energy 2010;35:3848e56. [15] Papanikolaou EA, Venetsanos AG, Heitsch M, Baraldi D, Huser A, Pujol J, et al. HySafe SBEP-V20: numerical studies of release experiments inside a naturally ventilated residential garage. Int J Hydrogen Energy 2010;35:4747e57. [16] Swain MR, Swain MN. Passive ventilation systems for the safe use of hydrogen. Int J Hydrogen Energy 1996;21:823e35.

[17] Swain MR, Grillot ES, Swain MN. Experimental verification of a hydrogen risk assessment method. Chem Health Saf 1999:28e32. [18] Barley CD, Gawlik K. Buoyancy-driven ventilation of hydrogen from buildings: laboratory test and model validation. Int J Hydrogen Energy 2009;34:5592e603. [19] Swain MR, Grilliot ES, Swain MN. Risks incurred by hydrogen escaping containers and conduits. In: Proceedings of the U.S. DOE hydrogen program review; 1998. NREL/CP-570-25315. [20] Swain MR, Grilliot ES, Swain MN. Dispersion of hydrogen clouds. Proceedings of the U.S. DOE hydrogen program review. NREL/CP-570e28890. [21] ANSYS FLUENT code, http://www.ansys.com/Products/ SimulationþTechnology/FluidþDynamics/Fluidþ DynamicsþProducts/ANSYSþFluent. [22] ANSYS CFX code, http://www.ansys.com/Products/ SimulationþTechnology/FluidþDynamics/ FluidþDynamicsþProducts/ANSYSþCFX. [23] ADREA-HF code, http://www2.ipta.demokritos.gr/pages/ ADREA-HF.html. [24] Wilkening H, Baraldi D, Heitsch M. CFD simulations of light gas release and mixing in the battelle model-containment with CFX. Nucl Eng Des 2008;238(3):618e26. [25] Papanikolaou E, Heitsch M, Baraldi D. Validation of a numerical code for the simulation of a short-term CO2 release in an open environment: effect of wind conditions and obstacles. Int J Hazard Mater 2011;190(1e3):268e75. [26] Papanikolaou E, Baraldi D, Kuznetsov M, Venetsanos A. Evaluation of notional nozzle approaches for CFD simulations of free-shear under-expanded hydrogen jets. Int J Hydrogen Energy 2012;37(23):18563e74. [27] Markert F, Melideo D, Baraldi D. Numerical analysis of accidental hydrogen releases from high pressure storage at low temperatures. Int J Hydrogen Energy 2014;39:7356e64. [28] Venetsanos AG, Papanikolaou E, Hansen OR, Middha P, Garcia J, Heitsch M, et al. HySafe standard benchmark problem SBEP-V11: predictions of hydrogen release and dispersion from a CGH2 bus in an underpass. Int J Hydrogen Energy 2010;35(8):3857e67. [29] Menter FR, Langtry RB, Likki SR, Suzen YB, Huang PG, € lker S. A correlation based transition model using local Vo variables part 1 e model formulation. ASME-GT2004-53452. Vienna, Austria: ASME TURBO EXPO; 2004. [30] Langtry RB, Menter FR, Likki SR, Suzen YB, Huang PG, € lker S. A correlation based transition model using local Vo variables part 2 e test cases and industrial applications. ASME-GT2004-53454. Vienna, Austria: ASME TURBO EXPO; 2004. [31] Langtry RB, Menter FR. Transition modeling for general CFD applications in aeronautics. AIAA paper 2005-522. 2005. [32] Reynolds O. An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels. Proc R Soc Lond 1883;35(224e226):84e99. [33] Molkov V, Shentsov V. Numerical and physical requirements to simulation of gas release and dispersion in an enclosure with one vent. Int J Hydrogen Energy 2014;39(25):13328e45. [34] Benintendi R. Laminar jet modelling for hazardous area classification. J. Loss Prev Process Ind 2011;24:123e30. [35] Barth TJ, Jesperson DC. The design and application of upwind schemes on unstructured meshes. AIAA Paper-890366. 1989. [36] Venetsanos AG, Papanikolaou E, Bartzis JG. The ADREA-HF CFD code for consequence assessment of hydrogen applications. Int J Hydrogen Energy 2010;35(8):3908e18.

i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 2 4 1 5 e2 4 2 9

[37] Hanna SR. Confidence limits for air quality model evaluations as estimated by bootstrap and jackknife resampling methods. J Atmos Environ 1989;23:1385e98. [38] Hanna SR, Chang JC, Strimaitis DG. Hazardous gas model evaluation with field observations. J Atmos Environ 1993;27A:2265e85.

2429

[39] Chang JC, Hanna SR. Air quality model performance evaluation. Meteorol Atmos Phys 2004;87:167e96. [40] Linden PF, Lane-Serff GF, Smeed DA. Emptying filling boxes: the mechanics of natural ventilation. J Fluid Mech 1990;212:309e35.