CFD analyses of large-scale tests highlighting the effects of thermal radiation in steam-containing atmospheres

In simulations of heat transfer in large gaseous volumes under moderate temperatures, thermal radiation is traditionally neglected to simplify the computations. In recent years however, some small as well as large-scale experiments have demonstrated that thermal radiation may have a strong effect on the temperature field if the atmosphere contains an infrared absorbing gas such as vapor, even in small concentrations. This has in particular been shown in the H2P2 series of tests performed in the large-scale PANDA facility at the Paul Scherrer Institute. For these tests an air-helium-steam mixture was compressed with a helium injection, which caused the formation of a hot bubble. The initial vapor volume fraction varied between 0.1 % and 60%. In this paper, we report the results of CFD simulations conducted for the five H2P2 tests. The turbulence model k - ω SST was used, while thermal radiation was considered using the simple P1 model as a default, in conjunction with the Weighted Sum of Gray Gases Model (WSGGM) to specify steam absorptivity. For comparison purposes we also made use of the more advanced Monte-Carlo model for selected experiments. CFD Best Practice Guidelines were employed, in particular through preliminary grid and time-step sensitivity studies. The temperature and helium concentration fields matched the experimental data quite accurately, except for the test with very low vapor content (H2P2_1_2, 0.1% vol. steam), where the temperature was somewhat over-estimated. This may be due to the inadequacy of the radiation models at these extreme conditions. Two major conclusions can be drawn from this study: 1) ignoring thermal radiation can lead to large over-predictions of temperature, even when the steam content is low (order of 1%) and 2) the simple P1 model provides sufficient accuracy. Use of the more sophisticated Monte-Carlo model yielded similar and slightly more accurate results, albeit at a considerably larger CPU expenditure (5 – 7 times).


Background on thermal radiation in enclosures
During real (Fujisawa et al., 2021) or postulated nuclear plant accidents (Ercosam, 2014), large amounts of steam and non-condensable gases such as hydrogen might be discharged into the containment.The volume pressurization as well as hydrogen distribution depends among others on the condensation heat transfer at the walls, the prevailing temperature and the spreading, mixing and transport of the gases.Hydrogen concentration may lead to combustion or detonation, hence endangering containment integrity and challenging internal structures and safety systems.Proper mitigation of such risks requires accurate knowledge of the accident thermal-hydraulics phenomena.
In containment flows, the temperature magnitudes are relatively small, and therefore radiative heat transfer has historically been neglected in thermal-hydraulics analyses.However, recent research suggests that thermal radiation can be a significant heat transfer mode whenever convective currents are weak (Soucasse et al., 2012), as is the case in reactor containments where natural circulation dominates.Thermal radiation takes places between walls held at different temperatures, and/or whenever an infrared absorbing gas such as steam is present in the atmosphere (Kapulla et al., 2022;Liu et al., 2022b;Kapulla et al., 2021).Previous computational analyses of experiments such as in TOSQAN (Filippov et al., 2016) have provided the first insights into the importance of thermal radiation for the overall heat transport processes.The inclusion of radiative heat transfer to and from the steam in the mixture resulted in the smoothing of the temperature field and absence of local peaks, in agreement with experimental data.
To highlight the need for including radiative heat transfer at low temperatures and natural circulation conditions, Dehbi et al. (Dehbi et al., 2019) made a series of CFD analyses.The authors addressed in particular wall-to-wall heat transfer in an air-filled differentially heated cavity (Kalilainen et al., 2016) as well as tests in the technical scale THAI facility (Kumar et al., 2020) with an air-steam mixture.The authors showed that a significant improvement in the predictions is achieved if thermal radiation is considered.A similar conclusion was reached by Kumar et al. (Kumar et al., 2020) based on their simulations of buoyancy-driven air flow in the THAI facility.Additionally, selected analyses of tests held in the large PANDA facility and featuring the erosion of a stratified layer by a vertical vapor jet have also pointed to the necessity of including radiative heat transfer for better predictive capabilities (Schmidt et al., 2014;Liu et al., 2022c;Liu et al., 2022a;Andreani et al., 2020;Kelm and Müller, 2016;Andreani et al., 2019).

Background on the HYMERES-2 PANDA tests highlighting thermal radiation effects
For the study of thermal radiation in enclosures, experiments such as those performed in the large-scale PANDA facility are preferred, since smaller scale experiments may suffer from inadequate scaling (Wang and Yan, 2021).Most of previous containment experiments have not been defined and conducted with the aim of evaluating radiative heat transfer.To remedy this, the HYMERES-2 H2P2 series (Paladino et al., 2022) are specifically designed to enhance the understanding of largescale thermal radiation phenomena and provide exhaustive data for lumped parameter and CFD code validation.
To draft the H2P2 test matrix, the conditions were defined to span the whole range where radiative heat transfer can be either neglected or is the dominant heat transfer mode.The experiments were performed in one of the PANDA vessels with a volume of 90 m 3 .
In this paper, we aim at validating CFD thermal radiation models against the extensive H2P2 test data.These tests involve simplified conditions such that convection is minimized, and hence, to the best of our knowledge, they constitute the first large scale containment experiments which allow a direct assessment of the thermal radiation models implemented in the CFD code used.
The paper is organized as follows: In Section 2, the experimental facility, test procedure and H2P2 test matrix are described.In Section 3, the CFD turbulence and radiation models are summarized.In Section 4, the physical modeling and numerical details of the simulations are provided.In Section 5, summary validation results and discussion are given.Finally, Section 6 provides conclusions on the work done.

The PANDA H2P2 experimental facility, procedure and test matrix
In this section, we give a brief description of the PANDA facility, the experimental procedure and outline the initial and boundary conditions for the five experiments conducted within the H2P2 series of the HYMERES-2 project.A more detailed discussion can be found in (Kapulla et al., 2022).The main parameters for the experiments are outlined in Table 1.The test matrix is spanned by i) the steam content, and ii) the initial gas and vessel temperature (T r : room temperature, T el : elevated temperature and iii) the initial pressure.
The experiment H2P2_1_2, starting at room temperature, T r , can be considered as the reference case with the lowest experimentally achievable steam content (<0.1 %).This experiment is complemented by test H2P2_2, which emphasizes the effects of a small increase of the steam content from 0.1 to 2 %.Experiments H2P2_3 and H2P2_4 demonstrate the importance of the initial vessel pressure on the temperature history, as both experiments were started at elevated temperatures and have the same steam content.Additionally, a comparison of H2P2_2 and H2P2_3 shows the influence of the initial temperature while keeping the same small steam content.Finally, experiment H2P2_5 is the closest to accident conditions, starting at elevated temperature and having a high initial steam content of 60 vol%.
A generic test procedure for all five experiments in terms of the helium flow rate (top) as well as the pressure (bottom) is provided in Fig. 1 and the experimental setup for H2P2 with the dimensions in mm is shown in Fig. 2. The initial temperature was adjusted according to the test matrix for the respective experiment and the initial gas composition consisted of dry air from the available compressed air supply.This dry air contains a minimal amount of humidity which was used and measured 'as is' for the reference experiment H2P2_1_2.For the other experiments, the humidity was adjusted by injecting steam (experiments H2P2_2 to H2P2_5).After the pre-conditioning phase, a helium stratification was built above the elevation of 6000 mm with a nominal helium content of X He ≈ 50 %, while the pressure was kept constant to avoid any preliminary temperature increase of the gas.The purpose of this stratification resides in 'isolating' the upward compression from the downward temperature increase, mimicking the motion of a piston.Fig. 3.
The main experimental procedure consisted in the injection of helium at room temperature with a flow rate of 10 g/s from the top through the man-hole of the vessel for a duration of 1200 s.This was accompanied by a transient gas temperature increase below the stratification.The temperature and pressure decayed after the helium injection was stopped, the latter phase being monitored for additional 1800 s.A circular disk with a diameter of 800 mm was placed 150 mm downwards of the helium injection to disperse the jet momentum and create a compression which is as homogeneous as possible, hence minimize the generation of any important convective currents.

Physical and turbulence models
The PANDA vessel is modeled as a 15 • wedge with two vertical symmetry planes representing 1/24 of the total volume.The mixture is assumed to be an ideal gas, with properties computed using mixture rules and the kinetic theory of gases.The URANS k-ω SST model of turbulence is used, as it shows good predictive capabilities for buoyancy

Table 1
Test matrix for the five experiments within the H2P2 series.The important initial conditions are highlighted in bold.T r and T el refer respectively to room and elevated temperatures at the start of the test.driven flows, among others.In addition, the turbulence generation due to buoyancy is included in both the kinetic energy and dissipation equations.The inclusion of the buoyancy effect on turbulence production has been shown to be important in flow configurations similar to our problem (Visser et al., 2014).

Radiation and absorptivity models
To account for radiation heat transfer both inside the gas mixture and between the mixture and the wall, a conservation equation is necessary.For a medium that undergoes radiative absorption, emission and scattering at position r → and direction s → ′, the radiative transfer equation (RTE) can be written as follows (ANSYS, Inc., 2022): The following definitions are used to apply eq. ( 1): s →′ : scattering direction vector.
s: path length.a: absorption coefficient.n: refractive index.Different models propose ways for solving the RTE.This is addressed next.

The P1 model
In the ANSYS Fluent code (ANSYS, Inc., 2022), six models are proposed for the modeling of radiation, be it with or without a participating medium.In this study, we endeavor to use the relatively simple P1 model as a default because of its computational efficiency.In selected cases, we will compare the predictions of the P1 model against the more sophisticated Monte-Carlo model (MC).The P1 model reduces the RTE to a diffusion equation for optically thick media.While simplifying the treatment, the P1 model may overestimate the magnitude of radiative effects (Kapulla et al., 2023).The model is nonetheless worth using in our context in order to investigate its performance in typical containment conditions.
When dealing with gray radiation, the P1 radiation model specifies the radiation flux as follows: where C depicts the linear anisotropic phase function coefficient.It is taken here to be zero as isotropic scattering is assumed by default.G is the incident radiation whose transport equation reads: S G is a user-defined source term, if any.The two equations above can be combined to give: This expression is inserted in the energy equation to account for radiative sinks or sources.By default, the refractive index is set to 1.

The Monte-Carlo model
The MC model is a general model meant to address a wide range of problems ranging from optically transparent to optically thick atmospheres typical of combustion.MC solves the integral form of the RTE, and is able to tackle even complicated conditions, e.g.scattering from particles, non-gray-wall emissivity etc.It is deemed to be among the most accurate radiation models, but this comes at significant CPU costs.
In the MC model, a photon is tracked from its source until its weight decreases below a minimum threshold and thus subsequently removed from the accounting.The physical parameters are updated following a A. Dehbi and R. Kapulla photon "event" (e.g.collision with a surface, scattering or absorption).Hence a complete "history" of the photon is kept in memory.A large number of photon histories are required to obtain accurate estimates of the radiation effects.The emitted radiation is sampled by photon sources in bands that are analyzed independently, thus accounting for radiation frequency dependence.
In the RTE equation, the MC model considers the radiation field as a photon gas and takes the radiation intensity to be proportional to the differential angular flux of the photons.The absorption coefficient a is therefore akin to the probability per unit length for a photon to be absorbed at a certain frequency.As a result, the mean radiation intensity I is proportional to the distance traversed by the photon per unit time in a unit volume located at r → .In interactions with surfaces, the radiative heat flux is proportional to the photons flux impinging on that surface.Proper accounting of the histories of a photon sample yields the mean total radiative intensity.

Specification of the absorptivity
The RTE requires the specification, among others, of the absorption coefficient, a, which is generally a function of temperature, pressure, gas composition and radiation frequency.The weighted-sum-of gray-gases model (WSGGM) is used to specify the absorption coefficient for the gas mixture.This approach is a reasonable middle ground between the simplistic gray gas model (constant for all radiation frequencies) and the more complicated model which involves all absorption bands.The WSGGM assumes that the total emissivity over distance is given by: In the above, a ε,i is the emissivity weighting factor of the i th gray gas, while the expression in the brackets represents the i th gray gas emissivity, κ i the absorption of the i th gray gas, p the sum of the partial pressures of all absorbing gases.The base case assumes a total pressure of one atmosphere.When this is not the case, scaling rules for κ i are employed.Value for a ε,i and κ i are obtained from (Coppalle and Vervisch, 1983).The mixture scattering coefficient σ s is taken to be zero as a default.
In the WSGGM, the absorptivity is estimated from the following relation: s is set to the mean beam length using the formulation of Siegel (Siegel and Howell, 1992): V is the fluid volume and A the total surface of the fluid boundaries.

Mesh independence
We strive in this work to use CFD Best Practice Guidelines (Casey and Wintergerste, 2000) to ensure that numerical errors are as small as possible.Second order accuracy is used for the momentum, energy and species transport equations, whereas first order accuracy is employed for the temporal advancement of the solution.Spatial resolution was checked first to ensure that the solution is independent of the grid used.Four mesh resolutions were used: the original (see Fig. 4), coarse, medium and fine grids, with cell counts of 64 K, 140 K, 280 K and 780 K, respectively.The original mesh is fully structured, with inflation layers at the wall to resolve the boundary layer and hence accurately predict the near-wall temperature gradients and accordingly the convective  heat transfer rates.Therefore, the y + values are of order 1 for the wallnearest cell.
The refinement was successively performed in the vessel section above y = 6 m because it is the region where the flow is most dynamic.Below that level, the flow consists of a plug being pushed downwards, similar to a piston motion, and thus increasing resolution there is wasteful.
The case H2P2_5 was simulated for the coarse, medium and fine meshes up to time t = 500 s.Comparison of the excess temperature and vertical component of velocity along the central axis is shown in Fig. 5.We define the excess temperature as the local temperature at any time minus the volume-averaged temperature of the gas mixture in the vessel at time zero.As can be seen, the medium and fine meshes provide virtually the same results, and hence the medium mesh is used for the rest of the simulations.

Time-step independence
Likewise, the sensitivity to the chosen time-step was investigated by running the simulation from time t = 500 s to t = 700 s with three time steps.As ANSYS Fluent is a fully implicit code, the developers recommend using a CFL of 20 or less in the most critical regions of the domain in order to guarantee time accuracy.This corresponds to time-steps around 0.12 s.It is demonstrated that the solution is insensitive to the time-steps in the range used, and as a result, the simulations were conducted with a Δt = 0.075 s, Fig. 6.

Initial and boundary conditions
The initial conditions (pressure, temperature, species volume fractions) are mapped into the numerical mesh following the test data.Likewise, the experimentally specified helium mass flow rate is imposed as a mass flow boundary condition at room temperature.Given the large thermal inertia of the vessel, a constant temperature is prescribed for the walls which are initially conditioned to have approximately the same temperature as the gas mixture.

Simulation results and discussion
Simulations were conducted for all five tests of the HYMERES-2 program.For the sake of conciseness, we restrict the presentation to the compression phase (t <1200 s) where radiation effects are most compelling.We present comparisons against the experimental data with and without radiation models.In the radiation simulations, the P1 model is used as a default.It is important to note that some of the axes scales in the figures are not shown or are presented in relative form, in conformity with the HYMERES-2 project guidelines.The project data will be made available to the public in 2024.For ease of comparison, the temperature scales in this section are identical for all tests.A. Dehbi and R. Kapulla

Qualitative description of the thermal and fluid fields
We first show some qualitative contours of the flow field for test H2P2_5 (60 % vol.steam content) midway through the compression sequence, at time 600 s.The helium molar fraction (red is 1 and blue is 0), which reflects a piston-like action of the injection process, with an upper region consisting of pure helium, and a lower region with no helium at all, is shown in Fig. 7.The magnitude of velocity is displayed in Fig. 8, with the helium jet penetrating approximately one third of the central vertical region, and otherwise a mostly stagnant volume.
The temperature evolution with and without radiation modeling in test H2P2_5 is shown next in Fig. 9. Drastically different temperature fields are predicted by the radiation and no-radiation models (note the different temperature scales).With no-radiation modeling, a hot bubble is formed and moves down the vessel due to the compression action.The upper atmosphere is cold because of the colder injected helium gas.On the other hand, when radiation is accounted for, the temperature varies only mildly in the vessel because of the homogenizing radiative heat exchange caused by the high steam content in the mixture.In the next section we discuss quantitatively the results of the five CFD simulations.

Test H2P2_1_2: Room temperature with less than 0.1 % steam by volume
The first test we present is H2P2_1_2 where the gas mixture is initially at atmospheric pressure and temperature and the vapor content is as low as technically possible, that is less than 0.1 % by volume.This can be considered a base, dry test against which other tests with variable steam content can be compared.The temperature fields midway (600 s) and at the end (1190 s) of the compression transient are shown in Fig. 10a and b, respectively.While, as expected, there is only a slight difference between the radiation and no-radiation simulations, the temperature is somewhat over-estimated in the hot bubble region.Correspondingly, the pressure is also slightly over-predicted by about 5 % towards the end of the transient (Fig. 10d).The reason for this will be discussed in a later section of this paper.The helium fraction profile along the vertical axis (Fig. 10c) is accurately captured.

Test H2P2_2: Room temperature with 2 % steam by volume
The second test we present is H2P2_2 where the gas mixture is initially at atmospheric pressure and temperature and the vapor content is low at 2 % by volume.The excess temperature profiles along the vertical axis at t = 600 s and t = 1190 s are shown in Fig. 11a and b, respectively.Despite the low steam content, radiation effects are noticeable, and a much better agreement with the data is achieved by taking radiation into account.
As time evolves, and the peak temperature increases due to compression, the discrepancy for the no-radiation cases goes up, whereas the agreement between the radiation case and the data is maintained.We note in addition that the temperature profiles for the noradiation case are essentially the same as for the previous dry test H2P2_1_2, as expected.The helium concentration is displayed in Fig. 11c at the end of the compression phase at t = 1190 s.The agreement with the data is very satisfactory, and this is the case for all ensuing simulations.Hence the helium concentration plots will not be shown again.The pressure history is displayed in Fig. 11d.The maximum difference is about 4.8 % at the end of the transient.The no-radiation prediction of the maximum pressure is slightly above the radiation value, i.e. 2 %, because of the higher temperature.

Test H2P2_3: High temperature with 2 % steam by volume
We present next the simulation for test H2P2_3 which has similar initial conditions as test H2P2_2 except that the mixture temperature is initially higher.As a result, thermal radiation effects are expected to be stronger.We show the respective evolution of temperature profiles along the axis of the vessel in Fig. 12.As in the previous tests, the agreement with the experimental data is very good when radiation is modeled, and conversely, the no-radiation simulations show substantially larger discrepancies with the data than in H2P2_2 with lower initial temperature.The discrepancies grow bigger with time as the hot bubble temperature increases and become unacceptably high at the end of the compression phase.

Test H2P2_4: High temperature, high pressure, with 2 % steam by volume
Test H2P2_4 is similar to H2P2_3 except that the initial pressure is at 3 bar.The temperature predictions including radiation are in good agreement with the data as shown in Fig. 13.On the other hand, neglecting radiation still leads to the over-prediction of the temperature  A. Dehbi and R. Kapulla in the hot bubble region, but it is less drastic than at lower pressure (Test H2P2_3).This higher pressure in H2P2_4 decreases the available (added) negative buoyancy difference between the injected cold helium and the helium stratification in the vessel dome.This forces the injected colder helium to spread more horizontally, promoting better convective mixing and decreasing the magnitude of the hot bubble temperature.

Test H2P2_5: High temperature with 60 % steam by volume
The last simulation we show is that of test H2P2_5 characterized by a large initial steam volume fraction of 60 % and a high starting temperature, conditions which are typical of post-blowdown containments.Because of this, thermal radiation effects are strong, which is clearly seen in Fig. 14.Excellent agreement with the experimental data is achieved when radiation is included.Disregarding radiation results in large over-prediction of the temperature, and this gets worse as the compression progresses.

Explaining the temperature inversion downstream of the baffle plate
In many tests (see hot test temperature Figures), there is a "temperature inversion" immediately downstream of the plate (at level close to 8 m).This trend is not predicted by the CFD simulations.The reason is as follows: it is known that for flows involving bluff bodies (Shirzadi et al., 2017), RANS models predict much longer wakes compared to experiments or Large Eddy Simulation (LES).In our particular problem, RANS predicts that the cold helium jet persists for a longer distance, whereas in reality, the wake is shorter and mixing between the cold wake fluid and hotter atmosphere results in higher temperature in the immediate region downstream of the plate.

Sensitivity of the predictions to radiation models: P1 vs Monte-Carlo
The adopted default P1 model of radiation is simple and computationally efficient.On the other hand, the Monte-Carlo model in ANSYS Fluent represents more closely the radiation physics but demands 9. Excess temperature field (in K) at different times in the compression transient.Above: No-radiation; Below: with radiation.significant computational effort.We show hereafter a sample comparison between the two models for tests where radiative effects are most prominent, namely H2P2_3 and H2P2_5.For the Monte-Carlo computation, 10 5 photon histories were tracked.This is a sufficiently large sample to produce accurate results.The temperature profiles at two snapshots in time are shown in Figs. 15 and 16, respectively.The predictions of both models are quite similar and close to the experimental data.For H2P2_5 with larger steam the Monte-Carlo model is slightly more accurate in the hot bubble region.Indeed, the P1 model slightly over-estimates the radiation effects, predicting more temperature homogeneity than the data, which is attributed to fact that the P1 model reduces the RTE to a diffusion equation as indicated in section 3.2.1.This confirms recent findings by Kapulla et al. (Kapulla et al., 2023) using a different CFD code.The Monte-Carlo profiles are somewhat "speckled".This could be remedied by significantly increasing the number of photons tracked, but it comes at a prohibitive CPU cost.
The agreement between the two models remains intact as time evolves.Although not shown here for the sake of brevity, similar results are obtained for other tests.We conclude therefore that the P1 model provides similar predictions as the MC model within the test conditions considered.The CPU expenditure is however 5 to 7 times larger for MC compared to the P1 model.

Peak pressure predictions
The difference in peak pressures between the CFD and data at the end of the compression is shown in Table 2.While the discrepancy is small (less than 6 %), the CFD prediction is systematically higher.One idealization of the CFD model accounts for a fraction of this over-estimation, and that is, the volume of the facility is actually 1.5 % higher than what is specified in the CAD files, due to "hidden" sub-volumes (connected pipework, flanges, etc.).Taking this into account, the maximum discrepancy would be of the order of 4 %, and is deemed satisfactory.A. Dehbi and R. Kapulla recommended that containment codes be updated to implement the P1 radiation model as a first step to ameliorate the accuracy of thermalhydraulic predictions in containment flows.

CPU expenditure
Future work should clarify the source of the overestimation of the hot bubble temperature at very low steam fractions (order of 0.1 %), as well as the small but systematic overprediction of the peak pressure.It is also recommended to simulate the H2P2 tests using a more advanced but less CPU-intensive radiation model such as Discrete Ordinates.A. Dehbi and R. Kapulla σ s : scattering coefficient.σ: Stefan-Boltzmann constant, 5.669 10 − 8 W/m 2 -K 4 .I: radiation intensity.T: temperature.Ф: Phase function Ω ′ : solid angle.

Fig. 1 .
Fig. 1.Generic test procedure for all five H2P2 experiments in terms of the helium mass flow rate (top) and the pressure (bottom).

Fig. 2 .
Fig. 2. One vessel of the PANDA facility used for the H2P2 series.Dimensions are in mm.The helium injection for the compression is done from the top.
Fig. 10.Test H2P2_1_2: a) Excess temperature at 600 s; b) Excess temperature at 1190 s; c) Helium molar fraction; d) Pressure time history.

Table 2
Difference in the peak pressure between the CFD predictions and the data.