Erosion of a helium-rich layer by a steam jet in the presence of an inclined grid: Comparison of the predictions by URANS, STRUC-URANS and LES

Severe accidents in Light Water Reactors (LWRs) involve the production of hydrogen, which may accumulate in the upper regions of the containment and form flammable mixtures with the initial air inventory. To study phenomena related to light gas layer removal, the H2P1_10_2 HYMERES test (Paranjape et al., 2019, 2020) was conducted in the PANDA vessel with the aim of investigating the erosion of a helium-rich layer by a steam jet deflected by an inclined grid. In this paper, we make use of the test data to validate computational fluid dynamics (CFD) simulations of such a transient. The simulation models involve the standard k-ε URANS, the more advanced second generation STRUC-ε URANS, and a coarse Large Eddy Simulation (CLES). All simulations are conducted with the same numerical grid consisting of 3.5 million hexahedral cells. Results show that all models capture well the main physical quantities, namely: the helium erosion rate, the temperatures and the velocities downstream of the obstructing grid. Fine details such as the turbulent kinetic energy (TKE) are better reproduced by the CLES, although the latter miss somewhat the amplitudes due to lack of spatial resolution, especially in the near-field downstream of the grid region. The STRUC-ε model replicates qualitatively the flow field fluctuations in sections of the domain where the spatiotemporal resolution is adequate. Mesh refinement in the jet region shows a substantial quantitative improvement, demonstrating that the STRUC-ε model is superior to traditional URANS when spatial resolution is increased. Finally, the simulations show unmistakably that radiation modeling is necessary to achieve reasonable predictions of the temperature field, as the radiative flux can account for a substantial fraction of the heat transfer to the vessel walls, in particular when convective currents are weak.


Introduction
In the course of a severe accident in a light water nuclear reactor (LWR), hydrogen is produced as a result of the oxidation of the zirconium fuel cladding at high temperatures.Pressure spikes from the combustion of the accumulated hydrogen may pose a significant threat to the integrity and function of the containment.In the TMI accident in 1979, it has been estimated that about 350 kg of hydrogen burned, fortunately without causing damage to the containment function, and thus caused no significant radiological consequences to the environment or the public (Kljenak et al., 2012).A series of hydrogen detonations, on the other hand, occurred in the Fukushima accident, damaging the reactors buildings and contributing to the radioactivity releases to the atmosphere (ANS, 2012;TEPCO, 2012;IAEA, 2015).This accident highlighted the importance and renewed the interest worldwide for investigating and further understanding hydrogen transport and mixing in the containment (OECD/NEA, 2013, 2014;Gupta, 2015;Nishimura et al., 2015) and evaluating mitigation measures (Hedley et al., 2014).Hydrogen management has been included also in the targeted safety reassessment of the nuclear power plants undertaken in Europe following the Fukushima accident (ENSREG, 2011).
The threat of a hydrogen combustion-detonation is particularly relevant where its concentration increases, as is the case if hydrogen stratification is formed in the dome of the containment (Kuznetsov et al., 2015).On the other hand, fluid flows in the containment -inertia or buoyancy driven e.g. a buoyant steam jet or natural circulation due to temperature differences and steam condensation on surfaces -may compete the development of the stratified layer or break-up an existing one (Studer et al., 2012).A number of experimental campaigns have been conducted with the aim to simulate hydrogen accumulation and erosion in large-scale facilities such as MISTRA (Studer et al., 2012), THAI (Gupta et al., 2015) and PANDA (Paladino et al., 2013) or in smaller ones e.g. the MiniPanda facility (Ritterath, 2012).For the sake of safety, experiments in these facilities use a non-combustible light gasmost commonly heliumas a surrogate for the combustible/explosive hydrogen.
The most recent program aimed at developing a better understanding of hydrogen risk is the OECD HYMERES (HYdrogen Mitigation Experiments for REactor Safety) project (OECD/NEA, 2016).The project provides high quality experimental data gathered in the PANDA facility (volume up to 460 m 3 ) to allow proper validation for lumped parameter, as well as computational fluid dynamics (CFD) codes.The simulations conducted in the present work deal with a HYMERES benchmark test, which will be described later in the manuscript.
As regards modelling, CFD codes have, in principle, a promising potential to assist in understanding and predicting the transport, mixing, stratification and break-up of the light gas layer, due to their high resolution in comparison with integral or lumped parameters codes (Hoyes and Ivings, 2016).To ensure a successful implementation of CFD modelling in reactor containment flows and support safety analyses, it is necessary that the codes and models be extensively validated.Experimental data available from facilities, such as those mentioned above, are vital in validating and evaluating the models.In this context, a number of modelling works for simulation of light gas stratification and erosion experiments have been performed.In these works, both tailor-made codes oriented to containment flows, such as GOTHIC (Cosials et al., 2016) or containmentFOAM (Kampili et al., 2021), as well as multipurpose CFD codes (Visser et al., 2012;Hoyes and Ivings, 2016;Kelm et al., 2016;Sarikurt and Hassan 2017;Abe et al., 2018) have been used.A successful simulation of the stratification erosion may be uncertain, since, depending on the modeling approach e.g.turbulence model (Abe et al., 2015) and overall setup, different users may produce quite different results (Andreani et al., 2019).Moreover, the large scale of experiments and thus the large computational domains often involved, along with the long simulation times create a significant obstacle in evaluating and optimizing the models, because the computational cost for the required sensitivity tests may be extensive (Andreani et al., 2019).
Capturing the physics of a light gas layer formation/erosion in the containment is intimately linked to the appropriate modeling of turbulence in the flow field.Turbulence modeling for engineering flows falls mainly into three categories: URANS (Unsteady Reynolds-averaged Navier-Stokes), LES (Large Eddy Simulation) and hybrid LES-URANS.
In URANS, the Navier-Stokes equations are averaged in time, which engenders closure terms that need to be modeled.While the URANS approach works satisfactorily in simple flows and geometries, its predictions are less accurate and may actually fail in highly transient flows or complicated geometries with strong anisotropic turbulence.
At the other end of the spectrum, the Large Eddy Simulation (LES) (Pope, 2000) is based on the spatial filtering of the Navier-Stokes equations whereby the large energy-containing eddies are properly resolved, whereas the smaller eddies are modeled using sub-grid scale (SGS) models.LES is very efficient at predicting complex transient flows, but requires high spatial resolution and short time steps, leading to CPU times that are typically two orders of magnitude larger than URANS.
Because of this, so-called hybrid URANS-LES models have been proposed with the aim of capturing the main flow physics while keeping the computational expenditure at acceptable levels, using coarser grids and larger time steps than LES.Up to the last few years, most hybrid URANS-LES approaches could be found in three flavors (Gopalan and Jaiman, 2015): Blending, Segregated and Interfacing approaches.In blending approaches, the hybrid residual stress tensor is computed by a weighted sum of the individual LES and URANS residual stresses.In segregated models, URANS is activated in the wall boundary layer region, while LES is used in bulk regions where large eddies are predominant.However, the transition between the two regions is not continuous and requires additional modeling, hence engendering mathematical ambiguities.Interface models follow the same concept of separating the URANS and LES regions.However, unlike the segregated approach, a single velocity transport equation applies in both regions, forcing the velocity field to be continuous.An exhaustive review of early hybrid models is provided by Fröhlich and von Terzi (2008).
In recent years, a new URANS modeling approach, dubbed STRUC or Second-Generation URANS, was proposed (Lenci, 2016;Baglietto et al., 2017).As the name suggests, this is a pure URANS concept.STRUC models attempt to resolve a sizeable fraction of kinetic energy without significantly altering the spatial resolution or the time step compared to classical URANS.For example, in STRUC-ε model, source terms are added in the dissipation equation (Xu, 2020) with the aim of increasing the dissipation rate wherever the local spatial and temporal resolution is not fine enough to capture the physics.Hence, modeling is added implicitly, yielding increased LES-like fluctuations in the most dynamic regions of the computational space.The STRUC models have shown significant improvement in predictive capabilities in a wide range of 3D problems (Lenci, 2016;Feng et al., 2018;Feng et al., 2020;Xu 2020;García et al., 2020).
In this paper, we seek to compare and contrast the predictions of the HYMERES helium erosion benchmark using three turbulence models, namely the standard k-ε, the STRUC-ε and the coarse LES (CLES).In Section 2, we provide information about the HYMERES program, the experimental test section, and the definition of the HYMERES test transient.In Section 3, we give a short description of turbulence models used.In Section 4, we show detailed predictions of the models versus important experimental measurements.In Section 5, we discuss sensitivity of the predictions on the He-steam diffusion coefficient and radiation modeling.Finally, in Section 6 we summarize the key findings of this study.

The HYMERES second phase benchmark
The HYMERES (HYdrogen Mitigation Experiments for REactor Safety) project is a multinational effort involving Paul Scherrer Institute (PSI) and Commissariat à l'Energie Atomique (CEA) Operating Agents, aimed at investigating the hydrogen risk phenomena in model containments.The project has three main components as detailed in OECD-NEA 2016.
Within the HYMERES project framework, the PANDA facility has been used to produce high quality experimental data in a number of areas.Examples of such studies include the erosion of a helium layer by a vertical air-helium jet (Kapulla et al., 2014), helium layer erosion by a horizontal steam jet (Paranjape et al., 2017), and natural circulation flow in a two-room type containment (Kapulla et al., 2018).
The HYMERES second phase benchmark experiment investigates the erosion of a helium-rich layer by a steam jet deflected by an inclined grid as shown in Fig. 1.In the initial phase, the PANDA facility is filled with steam at 1.3 bar, bringing the wall and gas temperatures to 103 • C to 108 • C. A helium-rich layer (25% helium by volume) is subsequently created in the region above elevation y = 6 m.The main phase is conducted by injecting superheated steam at 150 • C with a flow rate of 30 g/ s through a vertical tube of 0.2 m internal diameter.Venting of the vessel from the bottom is allowed in order to maintain the vessel pressure at 1.3 bar.An aluminum grid (side 0.96 m and unit element 28x28 mm 2 ) with an inclination of 17 • with respect to the horizontal is positioned between the helium layer and the outlet of the steam jet.This grid is intended to replicate effects an obstruction might have on the hydrogenrich layer erosion by steam jet in a hypothetical accident scenario.Parameters summarizing the transient are given in Table 1 and more details of the experiment can be found in Paranjape et al. (2019Paranjape et al. ( , 2020)).
The test lasts about 1600 s, i.e. until the helium-rich layer is fully eroded and the helium concentration in the vessel is roughly uniform.In the conditions of the test, no wall condensation takes place.
For the purposes of the benchmark, the PANDA facility was densely instrumented with thermocouples and helium gas concentration devices.In addition, a 2D Particle Image Velocimetry (PIV) system was used to measure velocities and turbulent kinetic energy (TKE) downstream of the grid.

Description of the turbulence models
We provide hereafter a summary of the three turbulence and P1 radiation models used in this work, as well as a short description of the meshing and numerical settings employed.

The k-ε model, including effects of buoyancy
The first model for turbulence resolution we use in this work is the URANS k-ε model.This model belongs to the well-known class of twoequation turbulence models, where turbulence is modelled by the solution of two additional transport equations for the turbulence scalar quantities.For the k-ε model, the transport equations for the turbulence kinetic energy (k) and the dissipation rate (ε) are solved.RANS models, including both the k-ε model and the k-ω model, i.e. the second representative model of the two-equation class, have been both used for buoyant jet and light gas erosion simulations (Abe et al., 2015;Abe et al., 2016;Kelm et al., 2016;Andreani et al., 2019).URANS models show the tendency to predict a wider turbulence transport in density interface and thus to overestimate the rate of erosion of the stratification (e.g.Abe et al., 2015).Inclusion of the effect of buoyancy on turbulence production may have a significant effect on the prediction of the erosion rate (Andreani et al., 2019).As shown in Visser et al. (2012), the two models performed similarly in the simulation of the THAI facility when the effect of buoyancy is included in the k transport equation (Visser et al., 2014).In these works, it was also found that a standard wall treatment with wall functions resulted in no significant differences compared to enhanced wall models, which allow for a full resolution of the boundary layer when the mesh is sufficiently fine.
In this work, the standard k-ε model is used with standard wall functions.ANSYS Fluent (ANSYS, 2021) offers the option to account for the generation of turbulence kinetic energy due to buoyancy through the addition of a relevant source term, G b , in the transport equation.For ideal gas mixtures, as assumed in this work, this term is written as follows: where g and ρ are, respectively, the gravity accelerationassumed to act in the vertical y direction-and the fluid density.The turbulent viscosity μ t is calculated from k and ε, and Pr t is the turbulent Prandtl number.By default, Fluent accounts only for the effect of the buoyancy in the production of k (the source term in the ε equation is zero).Here, we include the effect of buoyancy on the energy dissipation rate as well, through a non-zero source term (calculated from G b ) in the equation for ε, by enabling the relevant option in the code.This source term is calculated as follows: with C 1ε a constant with a default value equal to 1.44 and C 3ε = tanh(v/ u), where v and u are the components of the flow velocity parallel and perpendicular to the gravitational vector.The source term, hence, depends on the direction of the flow in relation to the gravity.For flows with a main direction aligned to the direction of gravity, C 3ε has a value of 1, while it is zero for a flow direction perpendicular to the direction of gravity (ANSYS, 2021).

The STRUC-ε model
One of the major drawbacks of the k-ε model is its tendency to overpredict the eddy viscosity in truly 3D situations.To remedy this, the STRUC model was proposed by Lenci (2016) and Baglietto et al. (2017).The approach resides in decreasing the eddy viscosity wherever the local turbulence time scale is larger than the modelled time scale.The STRUC model has shown consistently better accuracy than classical URANS in a large number of benchmarks.
The original STRUC model (Lenci, 2016) however shows heightened sensitivity to the inlet boundary conditions in external flows, which is detrimental to model robustness and generality.The more recent STRUCT-ε model by Xu (2020) is an improvement over the original STRUCT model.It uses the identical k-ε equations for mass, momentum, species and energy conservation, but introduces one additional source term in the ε equation which reads: The dissipation rate equation contains the added source termC ε3 ρk|Π|, which has units of inverse square of time.Hence, this term becomes important only where the modeled time scale k ε is much larger than the resolved time scale.The latter is defined in terms of the second invariant |Π| of the velocity tensor, also referred to as the Q-criterion which is used to characterize flow structures: where Above, Ω ij and S ij are respectively the rotation rate and strain rate tensors.The influence of the model constant C ε3 on the results is weak as reported by Xu (2020).We keep therefore the value of 1.5 as suggested in this reference.
In essence then, resolution of turbulent structures is performed locally and seamlessly based on comparison between the modeled and resolved time scales.In areas where the two time scales are comparable, the TKE is more resolved via viscosity reduction.In addition, the model is independent of inlet conditions, and thus is more general than the original STRUC model by Lenci (2016).
This model was implemented in the ANSYS Fluent code within a User Defined Function (UDF).The buoyancy effect on the generation of turbulence kinetic energy is employed in this model as well, as described in the previous section.

The LES model
For density varying flows, the ANSYS Fluent (ANSYS, 2021) code introduces the Favre averaging, which operates as follows for any property: where the overbar denotes spatial filtering.Assuming the steam-helium mixture to be an ideal gas, the density is given by: Above, p is the pressure, T the temperature, Ỹi the mass fraction of the i th species, M i the molecular weight of the i th species and R the universal gas constant.
In the LES approach (Pope, 2000) the spatial filtering operation of the Navier-Stokes equations results in residual Reynolds stresses that need to be closed using sub-grid scale (SGS) models.The space-filtered equations for the conservation of mass and momentum read: The stress tensor due to molecular viscosity μ is defined as: On the other hand, the subgrid scale tensor has the following form: The deviatoric SGS tensor is defined according to the compressible Smagorinsky model: Owing to the small velocities here, the term τ kk is neglected.The Smagorinsky model is employed to define the SGS eddy viscosity as follows: Above, the rate of strain tensor of the resolved velocity field is defined as: The length scale Δ is taken to be the cubic root of the cell volume.The procedure by Germano et al. (1996) is adopted to compute the Smagorinsky parameter C s in a dynamic fashion, which consistently produces superior predictions compared to the original Smagorinsky model with fixed constant C s .The energy equation reads: h s and λ are respectively the sensible enthalpy and thermal conductivity of the mixture.The SGS enthalpy flux is modeled in a fashion similar to the residual stress tensor: A similar equation for the species transport equation is used (omitted here for the sake of conciseness).In both energy and species equations, the SGS Prandtl and Schmidt numbers are obtained dynamically using the prescriptions of Germano et al. (1996).

The P1 radiation model
As indicated in Andreani et al. (2019), where several models for A. Dehbi et al. simulation of the helium erosion by a steam jet in the PANDA facility were used, thermal radiation has a significant effect on temperature profiles and time evolution prediction.This was shown as well in a similar study performed by Dehbi et al. (2018) on a THAI facility test, whereby a Monte Carlo radiation model was employed.Here, we include radiation in heat transfer modelling via the P1 radiation model.P1 model belongs to the general category of P-N models, where the radiation intensity is expanded into a series of spherical harmonics.In P1, the simplest for the P-N models, the radiation heat flux is calculated using only four terms in the series (ANSYS, 2021).Radiation is included in heat transfer calculations by adding its contribution, which is a function of the radiation heat flux and the incident radiation gradient, to the volumetric source term in the heat transport equation.Radiation heat transfer between solid surfaces and between mixture and wall surfaces is accounted for in the calculations.It is assumed that the absorption coefficient of the steam-helium mixture is equal to 1 m − 1 .This approximate value was found to give acceptable results in earlier calibration studies.The wall emissivities are taken to be 0.3 based on measurements conducted in the PANDA facility.The underlying assumption of the P1 model, i.e. that the medium is optically thick, are valid in this study as the gas mixture consists of 75-100% by volume of steam.

Meshing, numerical settings, properties
A mass flow rate inlet boundary condition is set for the steam injection, along with the experimentally supplied temperature history.A pressure outlet of 1.3 bars is prescribed for the outlet flow.In addition, due to the high thermal inertia of the vessel, the wall temperatures are set as isothermal surfaces for the whole transient.Conjugate heat transfer is allowed to take place between the aluminum grid and the fluid.
The numerical mesh was provided in the context of an international benchmark exercise (Kelm, 2019) and consists of 3.5 million hexahedral cells.The mesh is block-structured with O-grids and details are shown in Fig. 2. The grid plate, having roughly 1000 channels, consists of 0.9 million cells, including the solid and flow regions.The objective of this work is to compare the performance of turbulence models under identical conditions.Therefore, the same mesh ("common mesh") was used for the k-ε, STRUC-ε as well as the LES models.
For the URANS k-ε mesh resolution, credit is taken from the GCI performed by Kelm et al. (2017) on simulations conducted on the same vessel under similar test conditions, i.e. helium layer erosion by a steam jet in open space with no grid.Three mesh resolutions were used: 0.5 M, 2.2 M (reference) and 3.0 M (fine).It was demonstrated that the GCI with respect to the global erosion rate was quite small, 0.33%.Hence, in open space without a grid, a mesh side of 2.2 M nodes is appropriate.In the present conditions, with 3.5 M cells (0.9 M for to the grid), a similar resolution is justified.We conclude therefore that for the k-ε simulation, the spatial resolution is adequate.However, the resolution is somewhat coarse for the STRUC-ε which requires finer meshes than traditional URANS (García et al., 2020).The resolution is evidently coarse in the LES computation.Hence, the LES computation will be referred to as "Coarse LES", and will be subsequently called CLES.To see this more clearly, a contour plot, based on the k-ε simulation, is shown in Fig. 3 for the ratio of the grid length and the Taylor microscale length in the grid region.The ratio is defined as: As can be seen in Fig. 3, the mesh size is up to 8 times larger than the Taylor length scale, clearly signifying that the resolution is coarse compared to optimum LES requirements.
Concerning wall spatial resolution, the y + is in the range of 0-7 for the grid, with a mean of around 1. In the pipe, the mean y+ is about 7, and in the vessel it ranges from 0 to 11, with a mean of about 2. Apart A. Dehbi et al. from mesh resolution, the simulations follow CFD best practice guidelines (ERCOFTAC, 2000) to the extent possible so that numerical and user errors are minimized.For the CLES simulation, 2nd order spatial discretization is used for all equations, except for the momentum equation where bounded central difference is employed for stability purposes.Second order implicit method is employed for the time integration, with a time step of 0.01 s.This ensures that the Courant number is lower than 1 for the great majority of cells, except in the grid region where it reaches 5 in a very small number of cells in the region of jet impingement.
For the URANS simulations, 2nd order spatial discretization is used for all equations except the k and ε equations in which 1st order discretization is used to facilitate convergence.First order temporal integration scheme is used, with a fixed time step of 0.05 s, which translates into a Courant number less than 20 in almost all the domain as recommended by the code developers (ANSYS, 2021) to ensure time-accuracy of the results.
In the ANSYS Fluent code, mixture properties (C p , λ, μ) are computed using the mass-weighted mixing laws.As will be clarified later, the diffusion coefficient of He in steam is rather uncertain in the literature.A best-fit default value of 0.0003 m 2 /s is used based on prior experiments conducted in the PANDA facility.and CLES.Worth noting is the surprisingly smaller TKE predicted by the STRUC-ε model compared to the k-ε model.This is due to inadequate mesh resolution, and will be addressed later on in the article.In the remaining sub-sections, we discuss more quantitatively the models predictions versus.

Evolution of the helium layer
The primary objective of the HYMERES benchmark is the quantification of the time it takes for the steam jet to completely erode the helium layer.
We show on Fig. 8 the prediction of the helium molar fraction at the top-most point which is at elevation y = 8 m.All models predict well the history of the helium fraction, with a slight over-prediction of the erosion rate by the STRUC-ε model towards the end of the transient.We note also the more realistic fluctuating nature of the signal for the CLES and STRUC-ε models, whereas the signal is smoother for the classical kε model, as expected.
Somewhat below, at axial location y = 7 m, the predictions by all models are also quite good, and the signal fluctuations are more pronounced for the CLES.The STRUC-ε signal is somewhat damped, whereas the k-ε model shows no fluctuations at all due to the temporal filtering of the equations.For the off-axis measurements, the erosion rate is captured very well by all models in the region far from the jet (Fig. 9).In the jet area (Fig. 9 right), the erosion rate is slightly under-predicted.This can be traced to the model limitations in the case of k-ε.For the STRUC-ε and especially CLES, the under-prediction is most probably due to the lack of spatial resolution in the grid area, which yields lower levels of turbulence compared to the experimental data (as will be seen later), and hence less efficient mixing.

Velocity and turbulent kinetic energy
Experimental measurements for the gas velocity and velocity fluctuations are available for a region around the jet, downwind the obstacle grid (Paranjape et al., 2020) (Fig. 10).Measured data were obtained using PIV in periods equal to about 200 s around three different times throughout the experiment.Fig. 10 show the profiles of the velocity component perpendicular to the grid and the turbulent kinetic energy around 850 s.The profiles are shown along three lines parallel to the grid, at different distances from the top of the grid, namely, 20, 150 and 420 mm.For all three locations, the calculated velocity values and profiles are comparable for the different turbulence models.Some deviations can be observed, with the STRUC-ε model showing the narrowest jet with the highest maximum velocity magnitude.The higher velocity downstream the grid that is predicted by the STRUC-ε model is consistent with the velocity field in Fig. 4 and might partly be attributed to the higher velocity magnitude reaching the grid from below for the   STRUC-ε model.The lowest center velocity and the broadest jet is predicted by the CLES model, possibly due to the coarse computational mesh.CLES produces the lowest velocity also in the jet below the grid (Fig. 4).The velocity profile calculated with the k-ε model lies between those of the LES and STRUC-ε models.Discrepancies are observed between the calculated velocity and the experimental PIV measurements.Although the three models reproduce the shape of the velocity profile close to the grid, significant differences exist regarding the range of the velocity magnitude.In general, the peak values of the velocity are predicted satisfactorily, with accuracy that varies depending on the position.The minimum calculated velocity values are, however, significantly lower than the measured ones, throughout the whole window.Some discrepancies are also shown in the two other measurement lines further downstream, with all models failing to capture the rightwards shift of the jet as the distance from the plate increases.This is more clearly depicted in the more distant measurement line (Fig. 10  (d)).This may be related to a poor representation of the actual geometry by the computational mesh in the region of the grid.
In Fig. 11, the comparison between the predictions of the three models and the measurements for the turbulent kinetic energy is shown (for the CLES model, the resolved turbulent kinetic energy is presented).A decline of the turbulent kinetic energy is measured in the center of the jet, more clearly depicted in the line nearest to the grid plate (Fig. 11a).This seems consistent with the structure of the turbulent kinetic energy profile below the grid for all three models (see Fig. 5), however, none of the models captures closely this behavior downstream the grid, as it is shown in Fig. 11.The STRUC-ε and especially CLES models show a region with relatively constant average turbulent kinetic energy, whereas the k-ε model depicts a different profile, with the calculated magnitude increasing significantly towards the center of the jet.All models underestimate the turbulent kinetic energy in most of the measurement points.In the first line nearest the grid, the CLES and the k-ε models capture the order of magnitude of the turbulent kinetic energy, while the STRUC-ε model predicts about half the values of the other models.The differences between the calculated and the measured values become more significant further downstream, especially for the two URANS models, which may predict one order of magnitude lower turbulent kinetic energy values, compared with the experimental ones.Despite the coarse mesh, the CLES model captures significantly better the profile and magnitude of the turbulent kinetic energy, compared to the URANS models.In fact, the CLES model predicts some nonsymmetric generation of turbulence in the region of the measurement lines and agrees well with the measurements, especially at the more distant line (Fig. 11c).This higher turbulence is likely to promote a faster mixing and erosion of the density layer and may compensate the erosion delay due to the lower jet velocities calculated by the CLES model.This, in turn, might partly explain the steeper erosion slope history prediction and the better agreement with the measured erosion history slope in Figs. 8 and 9.

Mesh effects on the STRUC-ε predictions
The STRUC-ε model shows a significant under-prediction of the TKE, whereas one would expect a more accurate agreement with the data than for standard URANS k-ε.This stems from the fact that the mesh is quiet under-resolved in the grid region, as the STRUC model necessitates a finer spatial resolution compared to traditional URANS (García et al., 2020).Increasing the temporal resolution by a factor 5 (time-step of 0.01 s) gave predictions that were practically identical to the original computation, highlighting the fact that temporal resolution is sufficient for the coarse grid employed, and that a finer mesh resolution is required to better reproduce the experimental results near the grid.
To provide a direct evidence of this, the grid was refined in the region (x [-0.3 m-0.3 m], y [1.0 m-6.5 m], z [-0.3 m-0.3 m]).This encompasses the pipe region as well as the near field region of the grid where the jet is expected to be present.Refinement was performed by cutting each hexahedral cell into eight smaller cells.The resulting mesh has a total of 8.0 M cells, i.e. 2.3 times the original mesh.The simulation was conducted for the time interval 800 s to 900 s, with initial conditions at 800 s interpolated from the coarse mesh simulation.During this interval, the conditions are quasi-stationary, so that, for demonstration purposes, comparison between the coarse and fine mesh results is meaningful in the grid region, i.e. 20 mm and 150 mm away from the grid.The time step for the fine mesh simulation was 0.01 s, which results in a maximum CFL of order 5-10 in the jet impingement area, with a very small number of cell exceeding these values.
Figs. 12 and 13 show, respectively, the velocity and TKE plots for the coarse and fine meshes.The grid refinement ameliorates slightly the velocity predictions, but substantially the TKE results.As in previous investigations (García et al., 2020), it is found here that STRUC-ε provides superior predictions compared to classical URANS, but the grid needs to be refined to achieve this purpose.In the present situation, the grid refinement procedure was rather simplistic, and therefore it is most probable that an optimized grid will be substantially smaller than the 8.0 M cell grid used for this study.

Temperature field
The temperature upstream of the grid at elevation y = 5 m is shown in Fig. 14.The mean temperature is better predicted by the CLES.This stems from the fact that the CLES correctly predicts the jet break-up (see Fig. 6, left), resulting in an alternating pattern of cold and hot fluid chunks traversing the thermocouple.However, the rms of temperature is over-estimated by the CLES.This may be attributed to the substantial coarseness of the mesh compared to optimal resolution recommended for LES.On the other hand, both URANS models do not predict the jet Above the grid at axial location y = 6.5 m (Fig. 16), the models    close match of the experimental data is achieved by all models in this region where jet convection results in high fluid mixing rates specifically in the beginning of the simulation around 200 and towards the end of the simulation starting at 1200 s.At the time interval 800 s-1200 s, the models overestimate the temperature as shown in Fig. 17.The overestimation from the three models may be attributed to the previously mentioned underestimation of radiative effects.This overestimation occurs for a longer time than along the axis (Fig. 12), as the probe is farther away from the jet and at a higher helium density in the stratified layer.The turbulent fluctuations are present in the CLES data as well as the STRUC-ε, demonstrating the capability of the STRUC-ε model to reproduce such features.

Sensitivity studies on the diffusion coefficient and radiation effects
In this section, we look at the sensitivity of the results to the heliumsteam diffusion coefficient as well as to the radiation modeling.We use for such purpose the STRUC-ε as the base case model.

Sensitivity of the results to the helium-steam diffusion coefficient
Careful consideration was given to the specification of the He-Steam diffusion coefficient.In the literature, this coefficient is rather uncertain, with estimation at current conditions varying between 0.0001 m 2 /s and 0.0003 m 2 /s, depending on the correlation (Bird et al., 1960;Neufeld et al., 1972;Reid et al., 1988).A base-case value of 0.0003 m 2 /s was chosen based on experimental results of the OECD/SETH-2 tests conducted in the PANDA facility (Paladino et al., 2013).This chosen coefficient value was calculated using the (Bird et al., 1960) empirical correlation and was found to best match the experimental data acquired from the tests.The 0.0001 m 2 /s coefficient was computed from the kinetic theory and establishes the lower-end value in the work performed by (Paladino et al., 2013).In Fig. 18, we show the helium erosion rate for diffusion coefficient of 0.0003 m 2 /s and 0.0001 m 2 /s.As expected, the lower diffusion coefficient results in a sensibly slower erosion rate (≈ 200 s) compared to the base as well as experimental data.
Examining the helium erosion rate at a lower elevation at y = 7.478 m in Fig. 19, we find that the base-case STRUC-ε model underpredicts the molar fraction of helium between flow times 900 s and 1000 s.This may be attributed to a sensitivity of the radiation treatment, which becomes more apparent when the steam jet mixes with the helium-rich layer.Moving further down to y = 6.926 m as shown in Fig. 20, the effect of the diffusion coefficients diminishes starting at time 800 s, as turbulence becomes the driving mechanism for mixing.

Sensitivity of the results to the radiation modeling
In this section, we examine the effect of radiation modeling on the flow field.We first display on Fig. 21 the contour plot for the temperature field at time 850 s.It is clearly seen that omitting radiative heat transfer modeling (Fig. 21 Left) results in the formation of a hot fluid bubble in the upper section.This hot bubble is not observed experimentally, as can be seen in the plot on Fig. 22, where the temperature histories at two representative points are displayed.Thus, even the simple P1 radiation model yields acceptably accurate representation of the physics, while ignoring thermal radiation results in a substantial overestimation of the hot plume volume (see Fig. 21).
Finally, the total surface heat transfer rate with and without radiation is shown in Fig. 23.Negative values mean heat transfer from the fluid to the wall.In the first part of transient up to about 1100 s, convection is weak and hence neglecting radiation produces a substantial underprediction of the heat transfer from the fluid to the walls (up to a factor of 2).At the end of the transient, full mixing occurs, with high convective heat transfer coefficients to the walls.In this case, not including radiation results in higher total heat transfer to the walls   This sensitivity exercise shows clearly that in order to obtain accurate heat transfer and temperature field results, radiation must be modeled, especially when convective currents are weak.

Fig. 1 .
Fig. 1.Schematic of the PANDA vessel set-up and obstruction grid.
layer helium molar fraction at t = 0 0.25 A.Dehbi et al.

Fig. 2 .
Fig. 2. Details of the numerical grid: a) vertical mid-plane; b) side view of the grid region; c) top view of the grid region (reproduced from Kelm (2019)).

4.
Results of simulations by the k-ε, STRUC-ε AND CLES MODELS 4.1.Visualization of the flow field We begin by showing snapshots of the velocity, TKE, temperature, and helium concentration fields at time t = 850 s, which is slightly after the middle of the transient (Figs. 4 through 7).Qualitatively, the models show the expected behavior.The CLES displays intense fluctuations in regions of high shearing and turbulent mixing.At the other end of the modeling spectrum, the k-ε model shows rather smooth variations in the flow field.The STRUC-ε model features are somewhere in between k-ε

Fig. 3 .Fig. 4 .
Fig. 3. Mesh size in the grid region normalized by the Taylor length microscale (k-ε simulation).Left: grid region; Right: zoom on one channel in the grid.

Fig. 18 .
Fig. 18.Sensitivity of the He erosion rate to the helium-steam diffusion coefficient.

Fig. 19 .
Fig. 19.Sensitivity of the He erosion rate to the helium-steam diffusion coefficient at y = 7.478 m.

Fig. 20 .
Fig. 20.Sensitivity of the He erosion rate to the helium-steam diffusion coefficient at y = 6.926 m.
Fig. 21.Snapshot at 850 s showing the effect of radiation modeling.Left: no radiation; Right: P1 radiation model.

Fig. 22 .
Fig. 22.Effect of radiation modelling.Temperature history at two representative locations.Left: axial location; Right: off-axial location.

Fig. 23 .
Fig. 23.Total heat transfer rate to the wall with and without radiation modeling.

Table 1
Initial and boundary conditions of the HYMERES benchmark test.

Table 2
Comparison of CPU cost between the models.