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

NUCLEAR ENERGY SCIENCE AND ENGINEERING

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

Muhammad Saeed，

Ji-Yang Yu，

Aniseh Ahmed Atef Abdalla，

Xian-Ping Zhong，

Mahmood Ahmad Ghazanfar

Nuclear Science and Techniques

Vol.28, No.10

Article number 146

Published in print 01 Oct 2017

Available online 08 Sep 2017

DOI：10.1007/s41365-017-0304-x

36800

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 ModelHydrogen CombustionNuclear PowerPlant AccidentHYDRAGONAir-Fountain

1. Introduction

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-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-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-12].

New three-dimensional CFD codes were developed to provide more information about hydrogen behaviors in NPPs [13]. In 2009, Prankul Middha et al. [14-16] improved CFD code FLCAS for performing hydrogen dispersion calculations. In 2010, Heitsch et al. [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].

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

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

In this framework, a three-dimensional HYDRAGON code [22,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]. Helium gas, rather than hydrogen, was used in the experiment to avoid the risk of hydrogen explosion [25]. For this workout, reported experimental data of the air fountain case were selected as a benchmark.

2. Numerical methodology

2.1. Governing equations

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, 27].

\frac{\partial (ρ Y_{i})}{\partial t} + \nabla . (ρ_{α} U Y_{i}) = - \nabla . {\overset{⇀}{J}}_{i} + S_{i},

(1)

where, ρ is the fluid density of the gas mixture, Y is the mass fraction of gas species, ${\overset{⇀}{J}}_{i}$ 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:

\sum_{i}^{N C} Y_{i} = 1.

(2)

Diffusion flux for component i is:

{\overset{⇀}{J}}_{i} = - ρ D_{i} \overset{⇀}{\nabla} Y_{i} - \frac{μ_{t}}{s c_{t}} \overset{⇀}{\nabla} Y_{i} .

(3)

Conservation of momentum:

\frac{\partial (ρ U)}{\partial t} + \nabla . (ρ U U) = - \nabla P + \nabla . τ + ρ g .

(4)

Continuity mass equation:

\frac{\partial ρ}{\partial t} + \nabla . (ρ U) = 0,

(5)

where

\sum_{i = 1}^{N C} S_{i} = 0.

(6)

Stress tensor:

τ = μ_{e f f} (\nabla U - \frac{2}{3} ρ k I) .

(7)

Effective viscosity for turbulence:

μ_{e f f} = μ + μ_{T} .

(8)

Thermodynamic equation of gas:

p = z (ρ, T) ρ \frac{R}{M} T .

(9)

Density of the gas mixture:

\frac{1}{ρ} = \frac{Y_{a i r}}{ρ_{a i r}} + \frac{Y_{h e l i u m}}{ρ_{h e l i u m}} .

(10)

Molecular viscosity:

μ (T) = c_{1} {(\frac{T}{c_{2}})}^{1.5} \times \frac{c_{2} + 110.4}{T + 110.4},

(11)

where D_i is the mass diffusion coefficient of the gas species, NC is the number of components, SC_t= 0.7 [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

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, 30]. The turbulent kinetic energy is given as

k = \frac{1}{2} (\bar{U^{'} . U^{'}}),

(12)

Dissipation rate:

ε = \nabla . (ν U^{'} U^{'}) .

(13)

Production of turbulent due to the viscous force:

p_{k} = μ_{t} [\frac{\partial u_{i}}{\partial x_{j}} (\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}})] .

(14)

Production of turbulent due to the buoyancy effect:

G_{k} = - g_{i} \frac{μ_{t}}{ρ \Pr_{t}} \frac{\partial ρ}{\partial x_{i}},

(15)

where, Pr_t =0.9 is the turbulent Prandtl number.

2.2.1. Standard k-ε turbulence model

The standard k-ε turbulence model is a semi-empirical model, valid only for fully turbulent flows (high Reynolds number flows).

μ_{t} = C_{μ} ρ \frac{k^{2}}{ε},

(16)

SKE transport equation of $k :$

\frac{\partial (ρ k)}{\partial t} + \frac{\partial}{\partial x_{i}} (ρ k u_{i}) = \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{t}}{σ_{k}}) \frac{\partial k}{\partial x_{j}}] + P_{k} + G_{k} - ρ ε,

(17)

SKE turbulent dissipation rate ε:

\frac{\partial (ρ ε)}{\partial t} + \frac{\partial}{\partial x_{i}} (ρ ε u_{i}) = \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{t}}{σ_{ε}}) \frac{\partial ε}{\partial x_{j}}] + C_{1} \frac{ε}{k} (P_{k} + G_{k}) - C_{2} ρ \frac{ε^{2}}{k},

(18)

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,31], and C_μ=0, C₁=1.44, δ_k=1.0, δ_ε=1.3 and Pr_t=0.85 for the SKE turbulence model.

2.2.2. RNG k-ε turbulence model

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,33-35].

The eddy viscosity is calculated by

μ_{t} = C_{μ} ρ k^{2} / ε .

(19)

The RNG transport equation of k:

\frac{\partial (ρ k)}{\partial t} + \frac{\partial}{\partial x_{i}} (ρ k u_{i}) = \frac{\partial}{\partial x_{j}} [σ_{k} (μ + μ_{t}) \frac{\partial k}{\partial x_{j}}] + P_{k} + G_{k} - ρ ε .

(20)

The RNG turbulent dissipation rate ε:

\frac{\partial (ρ ε)}{\partial t} + \frac{\partial (ρ ε u_{i})}{\partial x_{i}} = \frac{\partial}{\partial x_{j}} [(μ + μ_{t}) σ_{ε} \frac{\partial ε}{\partial x_{j}}] + C_{1 R N G} \frac{ε}{k} (P_{k} + G_{k}) - C_{2} ρ \frac{ε^{2}}{k},

(21)

C_{1 R N G} = \frac{1.42 - η (1 - \frac{η}{4.38})}{(1 + 0.012 η^{3})},

(22)

η = \frac{k}{ε} S

(23)

S_{i j} = \frac{1}{2} (\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}}),

(24)

S = \sqrt{2 S_{i j} S_{i j}}

(25)

where C_μ= 0.0845, C₂=1.68 and σ_k=σ_ε=1.39.

2.2.3. Realizable k-ε turbulence model

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-40].

The RLZ transport equation of k:

\frac{\partial (ρ k)}{\partial t} + \frac{\partial}{\partial x_{i}} (ρ k u_{i}) = \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{t}}{σ_{k}}) \frac{\partial k}{\partial x_{j}}] + P_{k} + G_{k} - ρ ε .

(26)

The RLZ turbulent dissipation rate ε:

\frac{\partial (ρ ε)}{\partial t} + \frac{\partial (ρ ε u_{i})}{\partial x_{i}} = \frac{\partial}{\partial x_{j}} [(μ + \frac{μ_{t}}{σ_{ε}}) \frac{\partial ε}{\partial x_{j}}] + C_{1} ρ S ε + C_{1 ε} \frac{ε}{k} P_{k} - C_{2} ρ \frac{ε^{2}}{k + \sqrt{υ ε}},

(27)

μ_{t} = C_{μ} ρ \frac{k^{2}}{ε},

(28)

C_{μ} = \frac{1}{4.04 + A_{s} U^{*} \frac{k}{ε}},

(29)

where, σ_ε=1.2, C₂=1.9, C₁=max|0.43, η/(η+5)| (with η=kS/ε) and C_μ=1/(4.04+A_SkU*/ε), where A_S=6^1/2cosΦ, Φ=(cos⁻¹6^1/2W)/3,

$W = {\bar{S}}_{i j} {\bar{S}}_{j k} {\bar{S}}_{k i} / {\overset{⌣}{S}}^{3}, \overset{⌣}{S} = \sqrt{2 S_{i j} S_{i j}}, U^{*} = \sqrt{{\bar{S}}_{i j} {\bar{S}}_{i j} + {\bar{Ω}}_{i j} {\bar{Ω}}_{i j}}, {\bar{S}}_{i j} = \frac{1}{2} (\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}}) - \frac{1}{3} \frac{\partial u_{k}}{\partial x_{k}} δ_{i j}, and {\bar{Ω}}_{i j} = \frac{1}{2} (\frac{\partial u_{i}}{\partial x_{j}} - \frac{\partial u_{j}}{\partial x_{i}}) .$

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₂=1.9.

3. The facility and test conditions

As shown in 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.

Fig. 1

Sketch of the experimental design (in mm)

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]. The density sampling rate was 0.5 Hz in averaging time of 0.1 s. The parameters are summarized in Table 1

Parameters of the experiment.

Inlet diameter for helium (mm)	4
Inlet nozzle diameter for air (mm)	20
Containment sizes (m)	0.92(l)×1. 29(h)×0.92(w)
Air-injection time (s)	300
Pressure (Pa)	101235
Temperature (°C)	20
Air-fountain velocity (ms⁻¹)	2.803
Monitor elevation (P1/P2/···/P6) (m)	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] in Fig. 2.

Fig. 2

Helium mass fraction vs. sensor elevation obtained under the initial boundary conditions[21]

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⁻⁵ m²·s⁻¹ [26]. Fully inlet turbulence flow was calculated with k_in=0.001U_in² and ε_in= C_μk_in^3/2/l_in,,,where l_in=0.42 y_p.

By putting U_in=2.803 ms⁻¹ in Eq. (30), we have k_in=0.00785 m²·s⁻².

4. Results and discussion

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⁻¹ 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, 4 and 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

(Color online) Numerical results calculated by various models for the monitors of P1 (a) and P2 (b)

Fig. 4.

(Color online) Numerical results calculated by various models for the monitors of P3 (a) and P4 (b)

Fig. 5

(Color online) Numerical results calculated by various models for the monitors of P5 (a) and P6 (b)

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 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), 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. 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

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