An assessment of k-ε turbulence models for gas distribution analysis

An assessment of k-ε turbulence models for gas distribution analysis 1 corresp first-author Corresponding author: saeedphysics96@gmail.com Corresponding author: saeedphysics96@gmail.com saeedphysics96@gmail.com 1 corresp Corresponding author: saeedphysics96@gmail.com Corresponding author: saeedphysics96@gmail.com saeedphysics96@gmail.com 1 1 1 Department of Engineering Physics, Tsinghua University, Beijing 100084, China Department of Engineering Physics, Tsinghua University, Beijing 100084, China This paper presents the gas distribution analysis by injecting air-fountain into the containment and simulations with the HYDRAGON code. Turbulence models of standard k - ε (SKE), re-normalization group k - ε (RNG) and a realizable k - ε (RLZ) are used to assess the effects on the gas distribution analysis during a severe accident in a nuclear power plant (NPP). By comparing with experimental data, the simulation results of the RNG and SKE turbulence models agree well with the experimental data on the prediction of dimensionless density distributions. The results illustrate that the turbulence model choice had a small effect on the simulation results, particularly, the region near to the air-fountain source. Turbulence Model Hydrogen Combustion Nuclear Power Plant Accident HYDRAGON Air-Fountain 1. Introduction 1 During a severe accident in an NPP, a substantial amount of hydrogen may be generated and released into the containment. The density of hydrogen is fourteenth times lower than air and has a higher flammability range (4vol.%–74 vol.%). These distinctive features indicate that hydrogen can disperse extremely faster in a severe accident. The released hydrogen, produced as a result of oxidation and core degradation, can form a flammable mixture in the presence of oxygen in the containment. Accidental ignition of this mixture may initiate a combustion process. This might challenge the containment reliability and can damage the related safety systems of the NPP. The hydrogen risk was first observed in the Three Miles Island (TMI) accident in 1979, when a significant amount of hydrogen was released in the containment and initiated burning. The study of hydrogen distribution in the containment has been a main subject of interest after the TMI accident. After the Fukushima Dai-ichi NPP accident in 2011, modeling of the gas behavior became a significant topic of research [ 1 1 - 5 5 1-5 ]. It is important to measure the hydrogen concentration in the containment in accidental scenarios so as to evaluate possible hydrogen risk and efficiency of the mitigation systems installed in the containment. For years, two thermal-hydraulics approaches have been used for this purpose, including the codes of lumped parameters (LP) such as CONTAIN and MAAP, and the second approach of computational fluid dynamics (CFD). They were effective for establishing the averaged concentration of hydrogen in a compartment, but could not provide the information on the local distribution inside a compartment. The hydrogen detonation or deflagration is reliant upon local distribution of the hydrogen-air mixture and its temperature. Therefore, multidimensional CFD codes were introduced to treat the hydrogen distribution at lower spatial scales. CFD has significantly increased the accuracy and systematic analysis of hydrogen releases [ 6 6 - 8 8 6-8 ]. Three-dimensional CFD codes such as GOTHIC, GASFLOW, TONUS, FLCAS, ADREA-HF, REACFLOW, CFD-ACE and CFX-5.7 were designed for industrial and containment analyses. GOTHIC can be used to do LP computations or comprehensive multidimensional analysis. GASFLOW was designed for the simulating of hydrogen/steam distribution and combustion analysis in complex nuclear reactor containment geometries. TONUS was developed for hydrogen risk analysis, while TOPAZ and PHOENICS were developed to simulate hydrogen dispersion and release rate. These early designed codes restricted turbulence modeling, in their mesh capabilities and multiprocessor parallel performance as compared to the recent commercial CDF codes such as FLUENT and CFX [ 9 9 - 12 12 9-12 ]. New three-dimensional CFD codes were developed to provide more information about hydrogen behaviors in NPPs [ 13 13 ]. In 2009, Prankul Middha et al . [ 14 14 - 16 16 14-16 ] improved CFD code FLCAS for performing hydrogen dispersion calculations. In 2010, Heitsch et al . [ 17 17 ] made an assessment of the CFD codes FLUENT and CFX in postulated accident scenario to persuade the ability of the hydrogen mitigation devices. These provided a comparable quantity of hydrogen mitigation, but still failed to provide the local hydrogen concentration [ 18 18 ]. In 2004, HySafe Network of Excellence was initiated for making an assessment of the accuracy of CFD codes and their prospective applicability by comparing simulation results with experimental data. A large-scale project of CFD validation is in progress in the HySafe network. A series of experiments have been selected as Standard Benchmark Exercise Problems (SBEPs) to test the accuracy of CFD codes in capturing appropriate physical phenomena that occur in hydrogen accident scenarios [ 19 19 ]. The development of computational tools for investigating hydrogen-gas mixtures in the containment and predicting the behavior of the gas species was a major concern for the nuclear safety analysis. The location and mechanism of burnable cloud formation, mixing and their transport in the containment must be understood for mitigating the hydrogen combustion risk. Turbulence modeling was an important component for simulating gas mixing and transport equations. However, extra computational efforts were required to solve these turbulence transport equations [ 20 20 ]. Different turbulence models standard such as SKE turbulence model, RNG turbulence model, RLZ turbulence model and k - ω turbulence model were used to simulate mixed gas behavior [ 21 21 ]. In this framework, a three-dimensional HYDRAGON code [ 22 22 , 23 23 ] was developed to evaluate the thermal-hydraulic behavior of different gas species. To obtain reasonable and precise results, the simulation was quite sensitive to the turbulence models applied to the HYDRAGON code. This study was to test the influence of turbulence models to the gas distribution analysis and to demonstrate the code thermal-hydraulic simulation capability during NPP accident. Published experimental data named “air-fountain in the erosion of gaseous stratification” were used to study the break-up of a gaseous stratification by injecting air-fountain into the containment [ 24 24 ]. Helium gas, rather than hydrogen, was used in the experiment to avoid the risk of hydrogen explosion [ 25 25 ]. For this workout, reported experimental data of the air fountain case were selected as a benchmark. 2. Numerical methodology 1 2.1. Governing equations 2 The governing mass conservation and momentum equations were used to predict behaviors of incompressible multicomponent gas species. For gas distribution problems, the governing equations must be in transient form. Therefore, the governing equations conducted in this research were the unsteady average Navier-Stokes (N-S) equations. The fluid properties (temperature, pressure, velocity, etc.) in the multicomponent fluid were solved by transport equations. The transport equation was solved for each component as follows [ 26 26 , 27 27 ]. <math display="block" id="M1"><mrow><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><msub><mi>Y</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>.</mo><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mi>α</mi></msub><mi>U</mi><msub><mi>Y</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mo>−</mo><mo>∇</mo><mo>.</mo><msub><mrow><mover accent="true"><mi>J</mi><mo stretchy="true">⇀</mo></mover></mrow><mi>i</mi></msub><mo>+</mo><msub><mi>S</mi><mi>i</mi></msub><mo>,</mo></mrow></math> where, ρ is the fluid density of the gas mixture, Y is the mass fraction of gas species, <math id="M2"><mrow><msub><mrow><mover accent="true"><mi>J</mi><mo stretchy="true">⇀</mo></mover></mrow><mi>i</mi></msub></mrow></math> is the diffusion flux of gas species, U is the fluid velocity vector, and S i is the conserved mass that can be converted by phase changes from one state (liquid) to another state (gas) or split into a component by chemical reactions. In order to simplify numerical computations, the phase changes and chemical reactions were not considered. Simulation of Eq.(1) for each gas specie gets the general continuity mass equation, because the summation of source term ( S i ) of Eq. (1) for each gas specie equals to zero. Summation of mass fraction for each mixture component i is determined as: <math display="block" id="M3"><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mi>i</mi><mrow><mi>N</mi><mi>C</mi></mrow></msubsup><mrow><msub><mi>Y</mi><mi>i</mi></msub></mrow></mstyle><mo>=</mo><mn>1.</mn></mrow></math> Diffusion flux for component i is: <math display="block" id="M4"><mrow><msub><mrow><mover accent="true"><mi>J</mi><mo stretchy="true">⇀</mo></mover></mrow><mi>i</mi></msub><mo>=</mo><mo>−</mo><mi>ρ</mi><msub><mi>D</mi><mi>i</mi></msub><mover accent="true"><mo>∇</mo><mo stretchy="true">⇀</mo></mover><msub><mi>Y</mi><mi>i</mi></msub><mo>−</mo><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><mi>s</mi><msub><mi>c</mi><mi>t</mi></msub></mrow></mfrac><mover accent="true"><mo>∇</mo><mo stretchy="true">⇀</mo></mover><msub><mi>Y</mi><mi>i</mi></msub><mo>.</mo></mrow></math> Conservation of momentum: <math display="block" id="M5"><mrow><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>U</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>.</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>U</mi><mi>U</mi></mrow><mo>)</mo></mrow><mo>=</mo><mo>−</mo><mo>∇</mo><mi>P</mi><mo>+</mo><mo>∇</mo><mo>.</mo><mi>τ</mi><mo>+</mo><mi>ρ</mi><mi>g</mi><mo>.</mo></mrow></math> Continuity mass equation: <math display="block" id="M6"><mrow><mfrac><mrow><mo>∂</mo><mi>ρ</mi></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mo>∇</mo><mo>.</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>U</mi></mrow><mo>)</mo></mrow><mo>=</mo><mn>0</mn><mo>,</mo></mrow></math> where <math display="block" id="M7"><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>N</mi><mi>C</mi></mrow></msubsup><mrow><msub><mi>S</mi><mi>i</mi></msub></mrow></mstyle><mo>=</mo><mn>0.</mn></mrow></math> Stress tensor: <math display="block" id="M8"><mrow><mi>τ</mi><mo>=</mo><msub><mi>μ</mi><mrow><mi>e</mi><mi>f</mi><mi>f</mi></mrow></msub><mrow><mo>(</mo><mrow><mo>∇</mo><mi>U</mi><mo>−</mo><mfrac><mn>2</mn><mn>3</mn></mfrac><mi>ρ</mi><mi>k</mi><mi>I</mi></mrow><mo>)</mo></mrow><mo>.</mo></mrow></math> Effective viscosity for turbulence: <math display="block" id="M9"><mrow><msub><mi>μ</mi><mrow><mi>e</mi><mi>f</mi><mi>f</mi></mrow></msub><mo>=</mo><mi>μ</mi><mo>+</mo><msub><mi>μ</mi><mi>T</mi></msub><mo>.</mo></mrow></math> Thermodynamic equation of gas: <math display="block" id="M10"><mrow><mi>p</mi><mo>=</mo><mi>z</mi><mrow><mo>(</mo><mrow><mi>ρ</mi><mo>,</mo><mi>T</mi></mrow><mo>)</mo></mrow><mi>ρ</mi><mfrac><mi>R</mi><mi>M</mi></mfrac><mi>T</mi><mo>.</mo></mrow></math> Density of the gas mixture: <math display="block" id="M11"><mrow><mfrac><mn>1</mn><mi>ρ</mi></mfrac><mo>=</mo><mfrac><mrow><msub><mi>Y</mi><mrow><mi>a</mi><mi>i</mi><mi>r</mi></mrow></msub></mrow><mrow><msub><mi>ρ</mi><mrow><mi>a</mi><mi>i</mi><mi>r</mi></mrow></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><msub><mi>Y</mi><mrow><mi>h</mi><mi>e</mi><mi>l</mi><mi>i</mi><mi>u</mi><mi>m</mi></mrow></msub></mrow><mrow><msub><mi>ρ</mi><mrow><mi>h</mi><mi>e</mi><mi>l</mi><mi>i</mi><mi>u</mi><mi>m</mi></mrow></msub></mrow></mfrac><mo>.</mo></mrow></math> Molecular viscosity: <math display="block" id="M12"><mrow><mi>μ</mi><mrow><mo>(</mo><mi>T</mi><mo>)</mo></mrow><mo>=</mo><msub><mi>c</mi><mn>1</mn></msub><msup><mrow><mrow><mo>(</mo><mrow><mfrac><mi>T</mi><mrow><msub><mi>c</mi><mn>2</mn></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mn>1.5</mn></mrow></msup><mo>×</mo><mfrac><mrow><msub><mi>c</mi><mn>2</mn></msub><mo>+</mo><mn>110.4</mn></mrow><mrow><mi>T</mi><mo>+</mo><mn>110.4</mn></mrow></mfrac><mo>,</mo></mrow></math> where D i is the mass diffusion coefficient of the gas species, NC is the number of components, SC t = 0.7 [ 28 28 ]is the turbulent Schmidt number; τ is the Reynold stress tensor, p is the pressure, μ is the molecular viscosity, μ T is the turbulent viscosity, g is the gravitational vector, M is the fluid molecular weight, and R is the universal gas constant. 2.2. Turbulence models 2 The two governing equations (mass and momentum) are not closed. Therefore, the turbulence models are close to the system of the non-linear Reynolds stresses term. One transport equation for the turbulent kinetic k and another for turbulent dissipation rate ε are included. Therefore, eddy viscosity turbulence models are included in this study to close the Reynolds average equation [ 29 29 , 30 30 ]. The turbulent kinetic energy is given as <math display="block" id="M13"><mrow><mi>k</mi><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><mrow><mover accent="true"><mrow><msup><mi>U</mi><mo>′</mo></msup><mo>.</mo><msup><mi>U</mi><mo>′</mo></msup></mrow><mo stretchy="true">¯</mo></mover></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> Dissipation rate: <math display="block" id="M14"><mrow><mi>ε</mi><mo>=</mo><mo>∇</mo><mo>.</mo><mrow><mo>(</mo><mrow><mi>ν</mi><msup><mi>U</mi><mo>′</mo></msup><msup><mi>U</mi><mo>′</mo></msup></mrow><mo>)</mo></mrow><mo>.</mo></mrow></math> Production of turbulent due to the viscous force: <math display="block" id="M15"><mrow><msub><mi>p</mi><mi>k</mi></msub><mo>=</mo><msub><mi>μ</mi><mi>t</mi></msub><mrow><mo>[</mo><mrow><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mi>i</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mi>i</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mi>j</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mo>]</mo></mrow><mo>.</mo></mrow></math> Production of turbulent due to the buoyancy effect: <math display="block" id="M16"><mrow><msub><mi>G</mi><mi>k</mi></msub><mo>=</mo><mo>−</mo><msub><mi>g</mi><mi>i</mi></msub><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><mi>ρ</mi><msub><mrow><mi>Pr</mi></mrow><mi>t</mi></msub></mrow></mfrac><mfrac><mrow><mo>∂</mo><mi>ρ</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mo>,</mo></mrow></math> where, Pr t =0.9 is the turbulent Prandtl number. 2.2.1. Standard k - ε turbulence model 3 The standard k - ε turbulence model is a semi-empirical model, valid only for fully turbulent flows (high Reynolds number flows). <math display="block" id="M17"><mrow><msub><mi>μ</mi><mi>t</mi></msub><mo>=</mo><msub><mi>C</mi><mi>μ</mi></msub><mi>ρ</mi><mfrac><mrow><msup><mi>k</mi><mn>2</mn></msup></mrow><mi>ε</mi></mfrac><mo>,</mo></mrow></math> SKE transport equation of <math id="M18"><mrow><mi>k</mi><mo>:</mo></mrow></math> <math display="block" id="M19"><mrow><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>k</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>k</mi><msub><mi>u</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mi>μ</mi><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><msub><mi>σ</mi><mi>k</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mfrac><mrow><mo>∂</mo><mi>k</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac></mrow><mo>]</mo></mrow><mo>+</mo><msub><mi>P</mi><mi>k</mi></msub><mo>+</mo><msub><mi>G</mi><mi>k</mi></msub><mo>−</mo><mi>ρ</mi><mi>ε</mi><mo>,</mo></mrow></math> SKE turbulent dissipation rate ε : <math display="block" id="M20"><mrow><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>ε</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>ε</mi><msub><mi>u</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mi>μ</mi><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><msub><mi>σ</mi><mi>ε</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mfrac><mrow><mo>∂</mo><mi>ε</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac></mrow><mo>]</mo></mrow><mo>+</mo><msub><mi>C</mi><mn>1</mn></msub><mfrac><mi>ε</mi><mi>k</mi></mfrac><mrow><mo>(</mo><mrow><msub><mi>P</mi><mi>k</mi></msub><mo>+</mo><msub><mi>G</mi><mi>k</mi></msub></mrow><mo>)</mo></mrow><mo>−</mo><msub><mi>C</mi><mn>2</mn></msub><mi>ρ</mi><mfrac><mrow><msup><mi>ε</mi><mn>2</mn></msup></mrow><mi>k</mi></mfrac><mo>,</mo></mrow></math> where P k is the turbulence production due to the mean velocity gradients, G k is the influence of buoyancy force, σ k and σ ε are the turbulent Prandtl numbers[ 30 30 , 31 31 ], and C μ =0, C 1 =1.44, δ k =1.0, δ ε =1.3 and Pr t =0.85 for the SKE turbulence model. 2.2.2. RNG k-ε turbulence model 3 The RNG k - ε turbulence model is similar to the SKE turbulence model with various modifications on ε equation. However, the modifications of the RNG k - ε turbulence model are expected to make it more accurate and reliable for a wide range of flow types than the SKE turbulence model. The turbulence kinetic energy equation k for RNG model is identical to one of the SKE models, but the turbulence eddy dissipation ε and the constant values are different. Moreover, C 1RNG in the ε equation for RNG model is defined by Eq. (22). It is expected that these modifications enhance the accuracy for predicting air-flows [ 29 29 , 33 33 - 35 35 33-35 ]. The eddy viscosity is calculated by <math display="block" id="M21"><mrow><msub><mi>μ</mi><mi>t</mi></msub><mo>=</mo><msub><mi>C</mi><mi>μ</mi></msub><mi>ρ</mi><msup><mi>k</mi><mn>2</mn></msup><mo>/</mo><mi>ε</mi><mo>.</mo></mrow></math> The RNG transport equation of k : <math display="block" id="M22"><mrow><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>k</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>k</mi><msub><mi>u</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mrow><mo>[</mo><mrow><msub><mi>σ</mi><mi>k</mi></msub><mrow><mo>(</mo><mrow><mi>μ</mi><mo>+</mo><msub><mi>μ</mi><mi>t</mi></msub></mrow><mo>)</mo></mrow><mfrac><mrow><mo>∂</mo><mi>k</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac></mrow><mo>]</mo></mrow><mo>+</mo><msub><mi>P</mi><mi>k</mi></msub><mo>+</mo><msub><mi>G</mi><mi>k</mi></msub><mo>−</mo><mi>ρ</mi><mi>ε</mi><mo>.</mo></mrow></math> The RNG turbulent dissipation rate ε : <math display="block" id="M23"><mrow><mfrac><mrow><mo>∂</mo><mo stretchy="false">(</mo><mi>ρ</mi><mi>ε</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mo stretchy="false">(</mo><mi>ρ</mi><mi>ε</mi><msub><mi>u</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mo>=</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mi>μ</mi><mo>+</mo><msub><mi>μ</mi><mi>t</mi></msub></mrow><mo>)</mo></mrow><msub><mi>σ</mi><mi>ε</mi></msub><mfrac><mrow><mo>∂</mo><mi>ε</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac></mrow><mo>]</mo></mrow><mo>+</mo><msub><mi>C</mi><mrow><mn>1</mn><mi>R</mi><mi>N</mi><mi>G</mi></mrow></msub><mfrac><mi>ε</mi><mi>k</mi></mfrac><mrow><mo>(</mo><mrow><msub><mi>P</mi><mi>k</mi></msub><mo>+</mo><msub><mi>G</mi><mi>k</mi></msub></mrow><mo>)</mo></mrow><mo>−</mo><msub><mi>C</mi><mn>2</mn></msub><mi>ρ</mi><mfrac><mrow><msup><mi>ε</mi><mn>2</mn></msup></mrow><mi>k</mi></mfrac><mo>,</mo></mrow></math> <math display="block" id="M24"><mrow><msub><mi>C</mi><mrow><mn>1</mn><mi>R</mi><mi>N</mi><mi>G</mi></mrow></msub><mo>=</mo><mfrac bevelled="true"><mrow><mn>1.42</mn><mo>−</mo><mi>η</mi><mrow><mo>(</mo><mrow><mn>1</mn><mo>−</mo><mfrac><mi>η</mi><mrow><mn>4.38</mn></mrow></mfrac></mrow><mo>)</mo></mrow></mrow><mrow><mrow><mo>(</mo><mrow><mn>1</mn><mo>+</mo><mn>0.012</mn><msup><mi>η</mi><mn>3</mn></msup></mrow><mo>)</mo></mrow></mrow></mfrac><mo>,</mo></mrow></math> <math display="block" id="M25"><mrow><mi>η</mi><mo>=</mo><mfrac><mi>k</mi><mi>ε</mi></mfrac><mi>S</mi></mrow></math> <math display="block" id="M26"><mrow><msub><mi>S</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mi>i</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mi>j</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>,</mo></mrow></math> <math display="block" id="M27"><mrow><mi>S</mi><mo>=</mo><msqrt><mrow><mn>2</mn><msub><mi>S</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><msub><mi>S</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></msqrt></mrow></math> where C μ = 0.0845, C 2 =1.68 and σ k = σ ε =1.39. 2.2.3. Realizable k-ε turbulence model 3 Realizable k - ε turbulence model is similar to SKE turbulence model with modifications on turbulent viscosity and others. The feature of the realizable k - ε turbulence model shall be more accurate in predicting the spreading rate of various type flows such as; jet flows, flows involving rotation, boundary layers, separation flows and flows with strong pressure gradient [ 36 36 - 40 40 36-40 ]. The RLZ transport equation of k : <math display="block" id="M28"><mrow><mfrac><mrow><mo>∂</mo><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>k</mi></mrow><mo>)</mo></mrow></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mrow><mo>(</mo><mrow><mi>ρ</mi><mi>k</mi><msub><mi>u</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mo>=</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mi>μ</mi><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><msub><mi>σ</mi><mi>k</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mfrac><mrow><mo>∂</mo><mi>k</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac></mrow><mo>]</mo></mrow><mo>+</mo><msub><mi>P</mi><mi>k</mi></msub><mo>+</mo><msub><mi>G</mi><mi>k</mi></msub><mo>−</mo><mi>ρ</mi><mi>ε</mi><mo>.</mo></mrow></math> The RLZ turbulent dissipation rate ε : <math display="block" id="M29"><mrow><mfrac><mrow><mo>∂</mo><mo stretchy="false">(</mo><mi>ρ</mi><mi>ε</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><mo stretchy="false">(</mo><mi>ρ</mi><mi>ε</mi><msub><mi>u</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac><mo>=</mo><mfrac><mo>∂</mo><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mrow><mo>[</mo><mrow><mrow><mo>(</mo><mrow><mi>μ</mi><mo>+</mo><mfrac><mrow><msub><mi>μ</mi><mi>t</mi></msub></mrow><mrow><msub><mi>σ</mi><mi>ε</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mfrac><mrow><mo>∂</mo><mi>ε</mi></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac></mrow><mo>]</mo></mrow><mo>+</mo><msub><mi>C</mi><mn>1</mn></msub><mi>ρ</mi><mi>S</mi><mi>ε</mi><mo>+</mo><msub><mi>C</mi><mrow><mn>1</mn><mi>ε</mi></mrow></msub><mfrac><mi>ε</mi><mi>k</mi></mfrac><msub><mi>P</mi><mi>k</mi></msub><mo>−</mo><msub><mi>C</mi><mn>2</mn></msub><mi>ρ</mi><mfrac><mrow><msup><mi>ε</mi><mn>2</mn></msup></mrow><mrow><mi>k</mi><mo>+</mo><msqrt><mrow><mi>υ</mi><mi>ε</mi></mrow></msqrt></mrow></mfrac><mo>,</mo></mrow></math> <math display="block" id="M30"><mrow><msub><mi>μ</mi><mi>t</mi></msub><mo>=</mo><msub><mi>C</mi><mi>μ</mi></msub><mi>ρ</mi><mfrac><mrow><msup><mi>k</mi><mn>2</mn></msup></mrow><mi>ε</mi></mfrac><mo>,</mo></mrow></math> <math display="block" id="M31"><mrow><msub><mi>C</mi><mi>μ</mi></msub><mo>=</mo><mfrac><mn>1</mn><mrow><mn>4.04</mn><mo>+</mo><msub><mi>A</mi><mi>s</mi></msub><msup><mi>U</mi><mo>*</mo></msup><mfrac><mi>k</mi><mi>ε</mi></mfrac></mrow></mfrac><mo>,</mo></mrow></math> where, σ ε =1.2, C 2 =1.9, C 1 =max|0.43, η /( η +5)| (with η=kS / ε ) and C μ =1/(4.04+ A S kU* / ε ), where A S =6 1/2 cos Φ , Φ =(cos −1 6 1/2 W )/3, <math id="M32"><mrow><mi>W</mi><mo>=</mo><msub><mover accent="true"><mi>S</mi><mo>¯</mo></mover><mrow><mi>i</mi><mi>j</mi></mrow></msub><msub><mover accent="true"><mi>S</mi><mo>¯</mo></mover><mrow><mi>j</mi><mi>k</mi></mrow></msub><msub><mover accent="true"><mi>S</mi><mo>¯</mo></mover><mrow><mi>k</mi><mi>i</mi></mrow></msub><mo>/</mo><msup><mover accent="true"><mi>S</mi><mo>⌣</mo></mover><mn>3</mn></msup><mo>,</mo><mtext> </mtext><mover accent="true"><mi>S</mi><mo>⌣</mo></mover><mo>=</mo><msqrt><mrow><mn>2</mn><msub><mi>S</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><msub><mi>S</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></msqrt><mo>,</mo><mtext> </mtext><msup><mi>U</mi><mo>*</mo></msup><mo>=</mo><msqrt><mrow><msub><mover accent="true"><mi>S</mi><mo>¯</mo></mover><mrow><mi>i</mi><mi>j</mi></mrow></msub><msub><mover accent="true"><mi>S</mi><mo>¯</mo></mover><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>+</mo><msub><mover accent="true"><mi>Ω</mi><mo>¯</mo></mover><mrow><mi>i</mi><mi>j</mi></mrow></msub><msub><mover accent="true"><mi>Ω</mi><mo>¯</mo></mover><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></msqrt><mo>,</mo><mtext> </mtext><msub><mover accent="true"><mi>S</mi><mo>¯</mo></mover><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mi>i</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mo>+</mo><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mi>j</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>−</mo><mfrac><mn>1</mn><mn>3</mn></mfrac><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mi>k</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>k</mi></msub></mrow></mfrac><msub><mi>δ</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>,</mo><mtext> and </mtext><msub><mover accent="true"><mi>Ω</mi><mo>¯</mo></mover><mrow><mi>i</mi><mi>j</mi></mrow></msub><mo>=</mo><mfrac><mn>1</mn><mn>2</mn></mfrac><mrow><mo>(</mo><mrow><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mi>i</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>j</mi></msub></mrow></mfrac><mo>−</mo><mfrac><mrow><mo>∂</mo><msub><mi>u</mi><mi>j</mi></msub></mrow><mrow><mo>∂</mo><msub><mi>x</mi><mi>i</mi></msub></mrow></mfrac></mrow><mo>)</mo></mrow><mo>.</mo></mrow></math> Unlike the standard and RNG k - ε turbulence models where C μ is constant, the RLZ C μ is a function in Eq.(29). Furthermore, the effect of mean rotation is added to the turbulent viscosity. The model constants are σ ε =1.2, C 1 ε =1.44 and C 2 =1.9. 3. The facility and test conditions 1 As shown in Fig. 1 Fig. 1 , the facility is a box of 0.92 m( l )×1.29 m( h )×0.92 m( w ), with a Φ 40 mm air inlet in the bottom and the outlet gap of 6 mm all over the bottom plate. The environmental conditions are room temperature and normal pressure. In the first 300 s, 9.1 g of helium was injected from the two facing horizontal nozzles into the air filled enclosure, and the simulation was initiated. After 60 s of helium injection, air was injected, and the exceeding gas flew out of the containment to keep its constant pressure. The gas density was monitored at P1=1.29 m, P2=1.09 m, P3= 0.96 m, P4=0.83 m, P5=0.69 and P6=0.49 m from the containment bottom ( Z -axis is upward). The monitoring results agreed within 10% with the data measured in Ref. [ 21 21 ]. The density sampling rate was 0.5 Hz in averaging time of 0.1 s. The parameters are summarized in Table 1 Table 1 1 1 Inlet diameter for helium (mm) 1 1 4 1 1 Inlet nozzle diameter for air (mm) 1 1 20 1 1 Containment sizes (m) 1 1 0.92( l )×1. 29( h )×0.92( w ) 1 1 Air-injection time (s) 1 1 300 1 1 Pressure (Pa) 1 1 101235 1 1 Temperature (°C) 1 1 20 1 1 Air-fountain velocity (ms −1 ) 1 1 2.803 1 1 Monitor elevation (P1/P2/···/P6) (m) 1 1 1.29/1.09/0.96/0.83/0.69/0.49 To reduce the simulation time and simplify the computation, a 3D quarter size of the whole geometry was considered. In the computation, an isotropic system with the adiabatic heat transfer was considered. A square inlet section of (8.8 mm)×(8.8 mm) was defined for the inlet air-injection. The model was simulated with structural Cartesian coarse mesh and the computational domain cells were 15 cells along the X -axis, 15 cells along the Y -axis and 30 cells along the Z-axis. The highest mesh density was applied to the region near the air fountain source to capture the rapid change in the density gradient during air-injection. During the air-injection, small time steps were used to solve the transient equations due to the high turbulence intensity near the air-fountain source. Therefore, long simulation time was required. The mesh sensitivity analysis was performed and an average mesh size of 4.5 cm×4.5 cm×5 cm and 1.8 cm×1.8 cm×3.1 cm was tested, but it was found that the meshes did not affect the simulation results, hence the use of the coarse mesh. Convective terms were discretized by using the pure up-wind difference scheme. The initial helium mass fractions were obtained from the benchmark experimental data by E. Deri et al ., [ 21 21 ] in Fig. 2 Fig. 2 . The velocity of air-injection was 2.803 m/s from 150 s to 450 s. The mass diffusion coefficient of helium in the mixture was 7.35×10 −5 m 2 ·s −1 [ 26 26 ]. Fully inlet turbulence flow was calculated with k in =0.001 U in 2 and ε in = C μ k in 3/2 / l in ,,,where l in =0.42 y p . By putting U in =2.803 ms −1 in Eq. (30), we have k in =0.00785 m 2 ·s −2 . 4. Results and discussion 1 Initially, the helium molecules continuously diffused inside the enclosure and homogeneously accumulated on the top of the containment enclosure, due to the lower density of the helium. A 150 s later, the air-injection was started with an initial velocity of 2.803 ms −1 into the enclosure. The helium-injection was stopped after 450 s of the simulation initiated. The enclosure atmosphere was mixed due to the air-injection and the exceeding gas flow-out of the enclosure through the outlet kept the thermodynamic pressure constant. The simulations were carried out in two steps according to the flow type in the facility. During the release period (i.e., Helium injection and air fountain), small steps (Δ t ) were needed to solve the transient equations due to the turbulence intensity near the air-fountain source. Therefore, long simulation time was required during air-injection. While in the later phase, variable time steps were adopted. Furthermore, the “time step maximal variation” was set at low levels. The development of the helium stratification layers, the uniform helium concentration layer was observed at 100 s, while the stratification breakup caused by the impinging of air-injection on the density interface was observed at 200 s and 400 s. Furthermore, stratification layers driven by molecular diffusion force were observed again in the mixed gas at 1500 s. The computations obtained from the turbulence models depict in order to penetrate a helium stratification layer completely; air-injection did not have enough momentum . The contour of the mixture density computed with various turbulence models at 200 s during the air injection. The height of the air-jet penetrated the helium density interface was controlled by the buoyancy force and the momentum force of the injected air. The impact of the upward air-fountain flow distance was varied with k - ε turbulence models. Flow simulated with RLZ turbulence model reached a higher distance than the SKE and RNG turbulence models and the horizontal spreading rate obtained by RLZ turbulence model was underestimated. The predictions for the density stratification breakup phenomena were generally better represented by SKE and RNG turbulence models. Six elevations (0.49, 0.69, 0.83, 0.96, 1.09 and 1.29 m) were nominated to test the two equation turbulence models and the simulated results were compared. All the simulation results presented in Figs. 3 Figs. 3 , 4 4 and 5 5 were dimensionless density ( ρ ref − ρ )/ ρ ref versus time, where ρ ref is density of the air injected. The physical simulation time of these figures was started after the end of helium-injection in the benchmark case. It was found that the dimensionless densities of the mixture predictions by the three turbulence models were different. Fig. 3 Fig. 3 illustrates the comparison of simulation results for the monitors P1 and P2 with the experimental data. It was observed that the simulation results obtained with SKE and RNG turbulence models were almost similar and captured the experimental trend better than the RLZ turbulence model at any region in the containment. Similarly, Fig. 4 Fig. 4 illustrates the comparison for the monitors P3 and P4 with the experimental data. The simulation results obtained with SKE and RNG turbulence models captured the experimental trend better than the RLZ turbulence model at any region in the containment. From the Fig. 5(e) Fig. 5(e) , unlike the first 5 monitors, the RNG turbulence model had better-estimated values than the SKE and RLZ turbulence models for the monitor P6. Monitor P6 was located in the region close to the air-jet source and the flow included strong streamlines. The simulation results obtained with the RLZ turbulence model differed locally; its results have improved for monitors P3, P4 and P5 before and after air-injected, but the simulation results overestimated during the air-injected, for monitors P4, P5. Consequently, at monitor locations P1, P2 (the regions where the turbulence was not as strong as the turbulence near the source regions) results were underestimated as shown in Fig. 3 Fig. 3 . In general, computation results obtained by the RLZ turbulence model were improved locally depending on the turbulence intensity, consequently the results were underestimated at the region far away from the source. 5. Conclusion 1 In this paper, the gas dynamic behavior has been studied by using HYDRAGON code. Therefore, the numerical results of the air-fountain simulations and various turbulence models were studied using HYDRAGON code. The mixing process of multicomponent gas with air jet impinging into the containment was mainly caused by molecular diffusion force, gravity and momentum driven force. Before the air-injection and long after the end of the air- injection, molecular diffusion was responsible for the mixing process. Stratification layers were formed in-between the region where the helium concentration varied, as the helium concentrated at the upper part of the containment. During the air-injection, gravity and momentum driven force was responsible for the gas mixing. The simulation results obtained by SKE and RNG turbulence models have generally yielded better agreement with the experimental data at different elevations. Only in the region near the air-fountain source, the RNG turbulence model had better estimation values than SKE and RLZ turbulence models. Overall, the RNG model offered better predictions near the air-injection source. The calculated results obtained by RLZ turbulence model in the region far away from jet source showed discrepancies with the experimental data and the results improved only in some regions at specific transient times. In general, for the whole region, RLZ turbulence model could not precisely estimate the gas mixture behavior, especially in the region far away from the air jet source. RLZ turbulence model showed better estimation at the region where the turbulent flow was fully developed. The turbulence models influence was found very small on the results for the region near air-injection source. In contrast, the turbulence models influence on the results was found to be larger in the region far away from air-injection source, which was the region where turbulence was not as strong as the region near the source. It was concluded that the two equation turbulence models better estimated at the region where the turbulence was fully developed. References D.C. Visser , M. Houkema , N.B. Siccama et al ., Validation of a FLUENT CFD model for hydrogen distribution in a containment . Nucl. Eng. Des . 245 , 161 ( 2012 ). doi: 10.1016/j.nucengdes.2012.01.025 http://doi.org/10.1016/j.nucengdes.2012.01.025 Validation of a FLUENT CFD model for hydrogen distribution in a containment IAEA , Mitigation of hydrogen hazard in severe accident in nuclear power plants , International Atomic Energy Agency , Vienna ( 2011 ). IAEA TECDOC No. 1661 Mitigation of hydrogen hazard in severe accident in nuclear power plants J. Zhang , M.A. Delichatsios , A.G. Venetsanos , Numerical studies of dispersion and flammable volume of hydrogen in enclosures . Int. J. Hydrogen Energ . 35 ( 12 ), 6431 ( 2010 ). doi: 10.1016/j.ijhydene.2010.03.107 http://doi.org/10.1016/j.ijhydene.2010.03.107 Numerical studies of dispersion and flammable volume of hydrogen in enclosures J. Deng ,  X. Cao , Hydrogen and steam distribution following a small-break LOCA in large dry containment . Nucl. Sci. Tech . 18 ( 3 ), 181 ( 2007 ). doi: 10.1016/S1001-8042(07)60043-8 http://doi.org/10.1016/S1001-8042(07)60043-8 Hydrogen and steam distribution following a small-break LOCA in large dry containment Z.J. Xiao , S.Z. Qiu , W. B. Zhuo et al ., The development and verification of thermal-hydraulic code on passive residual heat removal system of Chinese advanced PWR . Nucl. Sci. Tech . 17 ( 5 ), 301 ( 2006 ). doi: 10.1016/S1001-8042(06)60057-2 http://doi.org/10.1016/S1001-8042(06)60057-2 The development and verification of thermal-hydraulic code on passive residual heat removal system of Chinese advanced PWR T. F. Kanzleiter , K. O. Fischer , Multi-Compartment hydrogen deflagration experiments and model development . Nucl. Eng. Des . 146 , 417 ( 1994 ). doi: 10.1016/0029-5493(94)90347-6 http://doi.org/10.1016/0029-5493(94)90347-6 Multi-Compartment hydrogen deflagration experiments and model development R.E. Henry , C.Y. Paik , M.G. Plys , MAAP4-Modular accident analysis program for LWR plants code manual , FAI , 1994 . MAAP4-Modular accident analysis program for LWR plants code manual H.A. Olvera , A.R. Choudhuri , Numerical simulation of hydrogen dispersion in the vicinity of a cubical building in stable stratified atmospheres . Int. J. Hydrogen Energ . 31 ( 15 ), 2356 ( 2006 ). doi: 10.1016/j.ijhydene.2007.08.028 http://doi.org/10.1016/j.ijhydene.2007.08.028 Numerical simulation of hydrogen dispersion in the vicinity of a cubical building in stable stratified atmospheres A.G. Venetsanos , T. Huld , P. Adams et al ., Source, dispersion and combustion modelling of an accidental release of hydrogen in an urban environment . J. Hazard. Mater . 105 ( 1 ), 1 ( 2003 ). doi: 10.1016/j.jhazmat.2003.05.001 http://doi.org/10.1016/j.jhazmat.2003.05.001 Source, dispersion and combustion modelling of an accidental release of hydrogen in an urban environment Z. Cheng , V. Agranat , A.V. Tchouvelev et al ., PRD hydrogen release and dispersion, a comparison of CFD results obtained from using ideal and real gas law properties . First international conference on hydrogen safety , Pisa, Italy , Sept. 8-10, 2005 . PRD hydrogen release and dispersion, a comparison of CFD results obtained from using ideal and real gas law properties H. Wilkening , D. Baraldi , CFD modelling of accidental hydrogen release from pipelines . Int. J. Hydrogen Energ . 32 ( 13 ), 2206 ( 2007 ). doi: 10.1016/j.ijhydene.2007.04.022 http://doi.org/10.1016/j.ijhydene.2007.04.022 CFD modelling of accidental hydrogen release from pipelines R. Fotis , S. Sklavounos , Evaluation of hazards associated with hydrogen storage facilities . Int. J. Hydrogen Energ . 30 ( 13 ), 1501 ( 2005 ). doi: 10.1016/j.ijhydene.2005.06.004 http://doi.org/10.1016/j.ijhydene.2005.06.004 Evaluation of hazards associated with hydrogen storage facilities J. Kim , S.W. Hong , S.B. Kim et al ., Three-dimensional behaviors of the hydrogen and steam in the APR1400 containment during a hypothetical loss of feed water accident . Ann. Nucl. Energy . 34 , 992 ( 2007 ). doi: 10.1016/j.anucene.2007.05.003 http://doi.org/10.1016/j.anucene.2007.05.003 Three-dimensional behaviors of the hydrogen and steam in the APR1400 containment during a hypothetical loss of feed water accident P. Middha , O.R. Hansen . CFD simulation study to investigate the risk from hydrogen vehicles in tunnels . Int. J. Hydrogen Energ . 34 ( 14 ), 5875 ( 2009 ). doi: 10.1016/j.ijhydene.2009.02.004 http://doi.org/10.1016/j.ijhydene.2009.02.004 CFD simulation study to investigate the risk from hydrogen vehicles in tunnels P. Middha , O.R. Hansen , I.E. Storvik , Validation of CFD-model for hydrogen dispersion . J. Hazard. Mater . 22 ( 6 ), 1034 ( 2009 ). doi: 10.1016/j.jlp.2009.07.020 http://doi.org/10.1016/j.jlp.2009.07.020 Validation of CFD-model for hydrogen dispersion P. Middha , O.R. Hansen , J. Grune et al ., CFD calculations of gas leak dispersion and subsequent gas explosions: validation against ignited impinging hydrogen jet experiments . J. Hazard. Mater . 179 ( 1 ), 84 ( 2010 ). doi: 10.1016/j.jhazmat.2010.02.061 http://doi.org/10.1016/j.jhazmat.2010.02.061 CFD calculations of gas leak dispersion and subsequent gas explosions: validation against ignited impinging hydrogen jet experiments M. Heitsch , R. Huhtanen , Z. Techy et al ., CFD evaluation of hydrogen risk mitigation measures in a VVER-440/213 containment . Nucl. Eng. Des . 240 ( 2 ), 385 ( 2010 ). doi: 10.1016/j.nucengdes.2008.07.022 http://doi.org/10.1016/j.nucengdes.2008.07.022 CFD evaluation of hydrogen risk mitigation measures in a VVER-440/213 containment D. M. Prabhudharwadkar , K.N. Iyer , N. Mohan et al ., Simulation of hydrogen distribution in an Indian nuclear reactor containment . Nucl. Eng. Des . 241 , 832 ( 2011 ). doi: 10.1016/j.nucengdes.2010.11.012 http://doi.org/10.1016/j.nucengdes.2010.11.012 Simulation of hydrogen distribution in an Indian nuclear reactor containment D. Baraldi , A. Kotchourko , A. Lelyakin , An inter-comparison exercise on CFD model capabilities to simulate hydrogen deflagrations in a tunnel . Int. J. Hydrogen Energ . 34 ( 18 ), 7862 ( 2009 ). doi: 10.1016/j.ijhydene.2009.06.055 http://doi.org/10.1016/j.ijhydene.2009.06.055 An inter-comparison exercise on CFD model capabilities to simulate hydrogen deflagrations in a tunnel J. Xiao , J.R. Travis , How critical is turbulence modeling in gas distribution simulations of large-scale complex nuclear reactor containments? , Ann. Nucl. Energy . 56 , 227 ( 2013 ). doi: 10.1016/j.anucene.2013.01.016 http://doi.org/10.1016/j.anucene.2013.01.016 How critical is turbulence modeling in gas distribution simulations of large-scale complex nuclear reactor containments? E. Deri , M. Bucci ,  E. Studer et al ., Analytical and computational analysis of turbulent buoyant jets in the containment atmosphere . Proceedings of the 16th International Conference on Nuclear Engineering , 2008 , ICONE16-48237: 689698 . doi: 10.1115/ICONE16-48237 http://doi.org/10.1115/ICONE16-48237 Analytical and computational analysis of turbulent buoyant jets in the containment atmosphere A.A. Abdalla , J. Yu , M. Alrwashdeh , Application of some turbulence models to simulate buoyancy-driven flow . Proceedings of the 2014 22nd International Conference on Nuclear Engineering ,( 2014 ). doi: 10.1115/ICONE22-30060 http://doi.org/10.1115/ICONE22-30060 Application of some turbulence models to simulate buoyancy-driven flow M. Saeed , J. Yu , B. Hou et al ., Numerical distribution of hydrogen inside a compartment using HYDRAGON code . Proceedings of the 2016 24th International Conference on Nuclear Engineering, ICONE24-60610 . doi: 10.1115/ICONE24-60610 http://doi.org/10.1115/ICONE24-60610 Numerical distribution of hydrogen inside a compartment using HYDRAGON code E. Deri , B. Cariteau , D. Abdo , Air fountain in the erosion of gaseous stratification . Int. J. Heat Fluid Fl . 31 , 935 ( 2010 ) doi: 10.1016/j.ijheatfluidflow.2010.05.003 http://doi.org/10.1016/j.ijheatfluidflow.2010.05.003 Air fountain in the erosion of gaseous stratification V. Agarant , Z. Cheng . A. Tchouvelev . CFD modeling of hydrogen releases and dispersion in hydrogen energy station . Proceeding of WHEC-15 , Yokohama , 2004 . CFD modeling of hydrogen releases and dispersion in hydrogen energy station D.C. Wilcox , Formulation of the k-ω Turbulence Model Revisited . AIAA J . 46 ( 11 ), 2823 - 2838 ( 2008 ). doi: 10.2514/6.2007-1408 http://doi.org/10.2514/6.2007-1408 Formulation of the k-ω Turbulence Model Revisited Z.U.A. Warsi , Fluid dynamics: theoretical and computational approaches . CRC press , c 2006 , 3 rd edition. ISBN 9780849333972-CAT#3397 Fluid dynamics: theoretical and computational approaches J.P.S. Sanders , B. Sarh , H. Gokalp , Variable density effects in axisymmetric isothermal turbulence jet: a comparison between a first and second order turbulence model . Int J. Heat Mass Transf . 40 ( 4 ), 823 ( 1997 ). doi: 10.1016/0017-9310(96)00151-2 http://doi.org/10.1016/0017-9310(96)00151-2 Variable density effects in axisymmetric isothermal turbulence jet: a comparison between a first and second order turbulence model M. Latif , C. Masson , T. Stathopoulos et al ., Comparison of various types of k-ε models for pollutant emissions around a two-building configuration . J. Wind Eng. Ind. Aerodyne . 115 , 9 - 21 ( 2013 ). doi: 10.1016/j.jweia.2013.01.001 http://doi.org/10.1016/j.jweia.2013.01.001 Comparison of various types of k-ε models for pollutant emissions around a two-building configuration Z. Xie , Y. Yang , H. Gu et al ., Numerical analysis of turbulent mixed convection air flow in inclined plane channel with k-ɛ type turbulence model . J Nucl. Sci. Techol . 19 ( 2 ), 121 ( 2008 ). doi: 10.1016/S1001-8042(08)60036-6 http://doi.org/10.1016/S1001-8042(08)60036-6 Numerical analysis of turbulent mixed convection air flow in inclined plane channel with k-ɛ type turbulence model W.P. Jones , B.E. Lanuder , The prediction of laminarization with a two-equation model of turbulence . Int J. Heat Mass Transf . 15 ( 2 ), 301 ( 1972 ). doi: 10.1016/0017-9310(72)90076-2 http://doi.org/10.1016/0017-9310(72)90076-2 The prediction of laminarization with a two-equation model of turbulence G. Biswas , V. Eswaran , Turbulent flow, fundemental, experiments and modeling ( England) : Alpha Science International Ltd, Pangbourne, Boca Raton, Fla: CRC Press ; c 2002 . Turbulent flow, fundemental, experiments and modeling M. Andreani , R. Kapulla , R. Zboray , Gas Stratification break-up by a vertical jet: Simula-tions using the GOTHIC Code . Nucl. Eng. Des . 249 , 71 ( 2012 ). doi: 10.1016/j.nucengdes.2011.06.004 http://doi.org/10.1016/j.nucengdes.2011.06.004 Gas Stratification break-up by a vertical jet: Simula-tions using the GOTHIC Code F.S. Lien , E. Yee , H. Ji et al ., Progress and challenges in the development of physically-based numerical models for prediction of flow and contaminant dispersion in the urban environment . Int. J. Computat. Fluid Dyn . 20 , 323 ( 2006 ). doi: 10.1080/10618560600898528 http://doi.org/10.1080/10618560600898528 Progress and challenges in the development of physically-based numerical models for prediction of flow and contaminant dispersion in the urban environment V.O. Yakhot , S.A. Orszag , S. Thangam et al ., Development of turbulence models for shear flows by a double expansion technique . Phys. of Fluids. A 4 , 1510 ( 1992 ). doi: 10.1063/1.858424 http://doi.org/10.1063/1.858424 Development of turbulence models for shear flows by a double expansion technique T.H. Shih , W.W. Liou , A. Shabbir et al ., A new k-ε eddy viscosity model for high Reynolds number turbulent flows, Comput . Fluids . 24 , 227 ( 1995 ). doi: 10.1016/0045-7930(94)00032-T http://doi.org/10.1016/0045-7930(94)00032-T A new k-ε eddy viscosity model for high Reynolds number turbulent flows, Comput W.P. Jones , B.E. Launder , The prediction of laminarization with a two-equation model of turbulence , Int. J. Heat Mass Transf . 15 , 301 ( 1995 ). doi: 10.1016/0017-9310(72)90076-2 http://doi.org/10.1016/0017-9310(72)90076-2 The prediction of laminarization with a two-equation model of turbulence A. Suryan , H.D. Kim , T. Setoguchi , Comparative study of turbulence models performance for refueling of compressed hydrogen tanks . Int. J. Hydrogen Energ . 38 , 9562 ( 2013 ). doi: 10.1016/j.ijhydene.2012.07.055 http://doi.org/10.1016/j.ijhydene.2012.07.055 Comparative study of turbulence models performance for refueling of compressed hydrogen tanks S.M. Tauseef , D. Rashtchian , S.A. Abbasi , CFD-based simulation of dense gas dispersion in presence of obstacles . Journal of Loss Prevention in the Process Industries . 24 ( 4 ), 371 ( 2011 ). doi: 10.1016/j.jlp.2011.01.014 http://doi.org/10.1016/j.jlp.2011.01.014 CFD-based simulation of dense gas dispersion in presence of obstacles T.R. Marrero , E.A. Mason , Gaseous diffusion coefficients . J. Phys. Chem. Ref. Data . 1 , 3 ( 2009 ). doi: 10.1063/1.3253094 http://doi.org/10.1063/1.3253094 Gaseous diffusion coefficients The authors highly acknowledge the support of the National key Lab of Reactor System Design Technology Chengdu, China. Mr. Muhammad Saeed also acknowledges the Chinese Scholarship Council for the award of Doctoral study. 10.1007/s41365-017-0304-x 2016-11-23 2017-02-04 2017-02-08 2017-09-08 2017-10-01 2021-04-11 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2017 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-017-0304-x http://dx.doi.org/10.1007/s41365-017-0304-x . 2017 Nuclear Science and Techniques 10 28