The Tailored CFD Package ‘containmentFOAM’ for Analysis of Containment Atmosphere Mixing, H₂/CO Mitigation and Aerosol Transport

Abstract

The severe reactor accident at Fukushima Daiichi Nuclear Power Plant (2011) has confirmed the need to understand the flow and transport processes of steam and combustible gases inside the containment and connected buildings. Over several years, Computational Fluid Dynamics (CFD) models, mostly based on proprietary solvers, have been developed to provide highly resolved insights; supporting the assessment of effectiveness of safety measures and possible combustion loads challenging the containment integrity. This paper summarizes the design and implementation of containmentFOAM, a tailored solver and model library based on OpenFOAM^®. It is developed in support of Research & Development related to containment flows, mixing processes, pressurization, and assessment of passive safety systems. Based on preliminary separate-effect verification and validation results, an application oriented integral validation case is presented on the basis of an experiment on gas mixing and H₂ mitigation by means of passive auto-catalytic recombiners in the THAI facility (Becker Technologies, Eschborn, Germany). The simulation results compare well with the experimental data and demonstrate the general applicability of containmentFOAM for technical scale analysis. Concluding the paper, the strategy for dissemination of the code and measures implemented to minimize potential user errors are outlined.

1. Background and Motivation

The severe reactor accident at Fukushima Daiichi (Japan, 2011) confirmed the need to understand the flow and transport processes inside the containment and connected buildings in more detail to design to assess the combustion risk. Over several years, CFD models, mostly based on commercial CFD solvers (e.g., [1]), have been developed to provide detailed insight and thus support the assessment of the effectiveness of safety measures and possible combustion loads that may challenge containment integrity. In the recent years, the open source CFD package OpenFOAM^® has been attracting the attention of academic and industrial users due to a wide choice of available models, numerical methods, and in particular its complete source code access. The latter is especially desirable, as it allows for a thorough review of existing code capabilities and paves the way for further model improvements. Even though tremendous efforts have been put into further development and validation of OpenFOAM^® for nuclear safety analysis, the works are always limited to specific and distinct phenomena, but not at an integral analysis of an accident sequence. Nevertheless, in other fields, e.g., fire safety research, such integrated simulation tools based on OpenFOAM^® exist. The fireFOAM code, developed by FM Global, aims at modeling fire spreading and suppression. It is comprised of robust models for the analysis of turbulent combustible gas and smoke transport, pyrolysis, thermal radiation, sprinklers, and fire suppression systems [2]. It is successfully used to disseminate knowledge and thus ensure improved design of fire mitigation and safety concepts.

Inspired by this, Forschungszentrum Jülich designed, implemented, and began validating a tailored solver containmentFOAM along with relevant physical and safety system model libraries to analyze pressurization, hydrogen transport, formation of combustible gas mixtures, and safety system performance during an accident in a nuclear power plant [3]. A generic application scenario is a loss-of-coolant accident (LOCA) as described e.g., in [4]. In such a scenario with core damage, the containment is the last remaining barrier against a release of radioactive material to the environment. In the long term, its pressurization (design pressure typically ~6 bar) is governed by the condensation of steam and the amount of non-condensable gases (O₂, N₂, H₂, CO, CO₂) in the containment atmosphere. In the short term, the local formation of a flammable gas mixture can result in dynamic combustion loads, which may challenge the integrity of the containment, its internal structures, and safety systems. To assess both pressurization paths, the 3D distribution of steam and non-condensable gases needs to be analyzed. Accidental flows and mixing processes inside the containment and connected buildings are wall bounded and dominated in the initial phase of a severe accident by inertial forces and later primarily driven by buoyancy forces. The density difference causing buoyant flows are either due to mass concentration gradients of the multicomponent gas mixture or temperature gradients due to conjugate heat transfer (convection and wall condensation), thermal radiation, local heat sinks, and sources (e.g., decay heat by aerosol transport, bulk condensation). Besides these physical phenomena, safety systems, e.g., passive auto-catalytic Hydrogen recombiners (PAR) or containment coolers interact with the surrounding containment atmosphere and introduce density differences. Vice versa, their efficiency is determined by the composition and temperature of the surrounding gas. Furthermore, new flow paths may open due to system conditions, e.g., pressure differences in case of doors or blow-out panels and thus, this system feedback needs to be modelled. All these characteristics are to be considered in the containmentFOAM package in a well-balanced level of detail to enable a numerically efficient and representative analysis of a full-scale containment flows and phenomena.

This paper summarizes the development outline and validation strategy of containmentFOAM, and highlights the status and capabilities of the solver (Section 2.1) and model library (Section 2.2). On this basis, the verification and validation strategy (Section 3.1) is presented and a recent application-oriented validation case is discussed: A PAR performance test, conducted under accidental conditions in the THAI facility is analyzed (Section 3.2). Section 4 summarizes the considerations for dissemination of containmentFOAM and measures implemented to minimize potential user errors, in particular a Graphical User Interface (GUI) for guided and standardized case setup according to well-tested procedures. Concluding, a summary and outlook is given in Section 5.

2. The containmentFOAM Package

There are different forks of OpenFOAM^® available. For the present work, the foundation version (openfoam.org) is selected as the development basis. The working version is built on OpenFOAM-6, while the transfer to the recent version (OpenFOAM-8/dev) is in progress. The containmentFOAM package consists of primarily two major parts, the tailored solver (Section 2.1) and model library (Section 2.2), which are summarized in the following. Besides, development and validation involve several pre- and post-processing tools, which are omitted for the sake of brevity.

2.1. Solver and Numerical Methods

The application of CFD to containment analysis must address large, multi- compartmented geometries, long transient durations, and many interacting physical phenomena. Despite the drastic increase of computing capabilities and a significant advancement of numerical methods, scale resolving simulations are still prohibitively demanding. Thus, an unsteady Reynolds Averaged Navier Stokes (U-RANS) approach is chosen, which does not require a highly resolved mesh and requires only two additional transport equations for modeling turbulence. Therefore, it can be directly scaled-up to containment scale without changing the representation of turbulent transport and dependent physics in the models. The multi-species solver development is discussed in detail in [5]. The governing equations are as follows.

Continuity equation:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \overrightarrow{U})={\dot{S}}_m^{‴}$$

(1)

Momentum conservation:

$$\frac{\partial (\rho \overrightarrow{U})}{\partial t}+\nabla \cdot \left(\rho \overrightarrow{U}\otimes \overrightarrow{U}\right)=-\nabla p+\nabla \cdot \tau +\rho \overrightarrow{g}+{\dot{S}}_u^{‴}$$

(2)

with the full buoyancy model $\rho g$ and the shear stress tensor:

$$\tau =\rho \left(\nu +{\nu }_t\right)\left[\nabla \overrightarrow{U}+{\left(\nabla \overrightarrow{U}\right)}^T-\frac{2}{3}\delta \nabla \cdot \overrightarrow{U}\right].$$

(3)

Species transport equations:

$$\frac{\partial \left(\rho Y_i\right)}{\partial t}+\nabla \cdot \left(\rho \overrightarrow{U}{Y}_i\right)=\nabla \cdot \left[\rho \left(D_{i,m}+\frac{{\nu }_t}{S{c}_t}\right)\nabla Y_i\right]+{\dot{S}}_y^{‴},$$

(4)

where Y_i is the mass fraction of the of species i and D_i,m is the diffusion coefficient of the species i in the mixture.

Heat transport is described by the total enthalpy equation,

$$\begin{array}{l}\frac{\partial \left(\rho h_{tot}\right)}{\partial t}+\nabla \cdot \left(\rho \overrightarrow{U}{h}_{tot}\right)=\rho \overrightarrow{U}\cdot g+\frac{\partial p}{\partial t}+\nabla \cdot \left[\left(\lambda +\frac{c_p{\nu }_t}{{\mathrm{Pr}}_t}\right)\nabla T\right]\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}+{\displaystyle \sum _{i=1}^n\nabla \cdot \left[\rho h_i\left(D+\frac{{\nu }_t}{S{c}_t}\right)\nabla Y_i\right]}+{\dot{S}}_h^{‴}\end{array}$$

(5)

with the total enthalpy $h_{tot}=h+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\left|\overrightarrow{U}\right|}^2$ and potential energy $\rho \overrightarrow{U}\cdot \overrightarrow{g}$. Furthermore, λ denotes the thermal conductivity and cp the isobaric specific heat capacity of the mixture.

In multi-component flows, the diffusion of enthalpy due to species diffusion (4^th term on the RHS) must be considered to obtain a consistent temperature field and was one of the major extensions to the underlying reactingParcelFOAM and chtMultiRegionFOAM solvers. The simple gradient diffusion hypothesis (SGDH) is used to model the turbulent heat and mass fluxes, where Sc_t = Pr_t = 0.9.

The model library (see Section 2.2) is partly implemented by means of sink and source terms in the governing equations, e.g., the volumetric source terms ${\dot{S}}_m^{‴}$, ${\dot{S}}_u^{‴}$, ${\dot{S}}_Y^{‴}$ and ${\dot{S}}_h^{‴}$ are utilized to implement models such as bulk condensation, thermal radiation, or to couple fluid flow with particle transport (see also [3]).

Conjugate heat transfer is modeled by solving an additional enthalpy equation,

$$\rho \frac{\partial h}{\partial t}=\nabla \cdot \left({\lambda }_s\nabla T\right)$$

(6)

for each structure and by enforcing a consistent heat flux and temperature values at the structure-gas interface (see also [6]).

OpenFOAM^® provides a variety of numerical schemes, methods, and linear solvers. To minimize user effects in further verification and validation of containmentFOAM, accurate discretization schemes, numerical solution parameters, and solver settings are defined as default values (refer to [3]). Unless mentioned otherwise, they are employed for verification and validation (V&V) cases. The default discretization schemes are second order accurate in time and space. Generally, V&V cases are run on high quality hexahedral meshes, i.e., with low non-orthogonality or skewness. However, the actual containment geometry is quite complex and only unstructured meshes, having significant amount of non-orthogonal cells, can be employed. Hence, future works with containmentFOAM will be with low quality unstructured meshes to confirm the model’s robustness and validity.

All the discretized equations are solved with the stabilized preconditioned bi-conjugate gradient (PBiCGStab) solver, except the pressure equation with the preconditioned conjugate gradient (PCG) solver. For large cell counts (>1 Mio. nodes), the Geometric Algebraic Multi-Grid (GAMG) solver proved to be more efficient than PCG. OpenFOAM^® uses a segregated approach for the solution of the governing equations, using the PIMPLE algorithm. In particular, the time step is chosen to maintain the mean CFL ≤ 1, while a local maximum CFL ≤ 10 is tolerated if outer loop convergence is satisfied. It must be stated that this criterion also allows to avoid under relaxation within a time step. Tests revealed that this can often lead to a reduced computational time, compared to the use of larger time-steps with under relaxation.

2.2. Model Library

Containment flows and mixing processes are essentially a multi-physics problem. All physical phenomena and technical systems must be considered in the full-scale application to represent a realistic accident sequence. This implies that both flow regime and local boundary conditions are changing frequently with transient progression and consequently it is hardly possible to identify the ‘best’ suitable model. Furthermore, containment flows involve multiple scales ranging from the physical scales up to the system size and accident duration. Consequently, a major requirement is a well-balanced modeling approach with a reliable representation of the physics as well as numerical efficiency.

To cope with these requirements, a baseline set of models for containment application has been defined at Forschungszentrum Jülich, for application in the CFD software ANSYS CFX (ANSYS Inc., Canonsburg, PA 15317, USA) [1]. It has been subject to systematic validation against separate effect tests and a variety of application relevant technical scale coupled effect tests (see Section 3.1). The containmentFOAM model library builds on this experience. As a first step, the baseline model has been transferred. To ensure comparability with previous work, available OpenFOAM^® models have been revised and extended (e.g., wall treatment or turbulence model) and missing models such as condensation were implemented. Furthermore, the baseline model is continuously extended to include for example radiative heat transfer or aerosol transport. The baseline model is summarized in

Table 1. Baseline model for containment applications.

Phenomenon	Model
Governing Equations (Section 2.1)	unsteady Reynolds-averaged Navier-Stokes equations, total energy equation (Non-unity Lewis number formulation), transport equation for species mass fraction
Turbulent Transport (Section 2.2.1)	k-ω based Shear Stress Transport (SST) model, including production and dissipation terms for buoyancy turbulence based on SGDH; Turbulent heat and mass transport: SGDH, Prt = Sct = 0.9
Multi-Species Transport (Section 2.2.2)	temperature dependent transport properties and specific heat capacity, Wilke mixture for transport properties, effective binary diffusion model, pressure and temperature dependent diffusion coefficients according to Fuller model*
Wall Treatment and Condensation (Section 2.2.3)	continuous wall functions wall condensation: single phase ‘diffusion-layer’ approach bulk condensation: single phase homogeneous equilibrium approach
Thermal Radiation (Section 2.2.4)	Monte Carlo RTE solver, importance sampling with 40–160 photon histories per cell, spectral models of participating media: gray (Planck mean absorption coefficient) and SNBCK, εwall = 0.6,
Aerosol Transport and Decay Heat (Section 2.2.5)	Euler-Lagrangian framework: Drag, lift and thermophoretic force, continuous random walk model Euler-Eulerian framework (in progress): Mixture model with corrections for particle inertia, Brownian diffusion for sub-micron particles, sectional modeling for evolution of aerosol size distribution due to steam condensation, agglomeration, de-agglomeration
System Models (Section 2.2.6)	porous media, conditional AMI for burst discs & doors coupling with mechanistic PAR model REKODIREKT

.

Further development and extension of the baseline model is ongoing. In the following sections, the baseline model and status and direction of ongoing developments are summarized.

2.2.1. Turbulence Transport in Buoyancy Affected Flows

Based on the previous validation experience [1], the standard k-ω SST turbulence model was adopted (SST-2003 [7]). It is the de-facto industrial standard and allows the use of logarithmic wall functions or integrate up to the wall.

To consider the production and dissipation of turbulence due to buoyancy, several models have been developed and validated [8]. Results obtained within the recent OECD/NEA PANDA CFD benchmark (Abe et al. [9]) highlighted again the need to consider the effect of buoyancy on turbulence production and dissipation to predict the erosion of a stratification accurately. Neglecting buoyancy turbulence production and dissipation leads to considerably faster erosion of a stable stratified gas layer. The simplest option is to model the buoyancy production in the turbulent kinetic energy by means of the Simple Gradient Diffusion Hypothesis (SGDH) [8]. The additional production term P_k,b is added as a volumetric source term in the transport equation of turbulent kinetic energy k:

$$P_{k,b}=-\frac{{\nu }_t}{{\sigma }_{\rho }}\left(\overrightarrow{g}\cdot \nabla \rho \right)$$

(7)

Here the turbulent Schmidt number is defined as ${\sigma }_{\rho }$ = 1.0. It is noted that in a stably stratified atmosphere ($\partial \rho /\partial y<0$), the source term is negative, while for an unstably stratified flow ($\partial \rho /\partial y>0$), excess turbulent kinetic energy is produced.

Correspondingly, a source term in the turbulent eddy frequency equation is defined:

$$P_{\omega ,b}=\frac{1}{{\nu }_t}\left[\left({\gamma }_1+1\right)C_3\cdot \max \left(P_{k,b},0\right)-P_{k,b}\right]$$

(8)

with ${\gamma }_1$ = 5/9 and $C_3$ = 1.

Ongoing work considers the Generalized Gradient Diffusion Hypothesis (GGDH) and a direction dependent dissipation term [8], which is expected to improve predictive capabilities in flows, where density gradients are not aligned with the gravity vector. Besides, low Reynolds number damping terms are investigated to enable more consistent prediction of near wall production and dissipation of turbulent kinetic energy, which is of high relevance, e.g., deposition of aerosols due to turbophoresis or mixed convection heat and mass transfer.

2.2.2. Multi-Species Transport

The major challenge in multi-species transport is the description of the species properties and the modeling of the mixture properties. In containmentFOAM, the fluid is assumed to be an ideal mixture of ideal gases. The specific heat capacity of the mixture is obtained by the mass weighted sum of individual species’ specific heat capacities. For each species, a temperature-dependent specific heat capacity is described by polynomials of the 4^th order, fitted to the NIST Chemistry Webbook [10] and VDI heat atlas [11] data.

Species transport properties $\phi$, dynamic viscosity $\eta$, and thermal conductivity are described by polynomials of 2^nd order, which were fitted to the NIST Chemistry Webbook data [10]. The mixture transport properties are obtained with Wilke’s mixture model [11], since it is accurate for gas mixtures involving light gases like Hydrogen or Helium.

$${\phi }_m^{}={\displaystyle \sum _i\frac{x_i{\phi }_i}{{\displaystyle \sum _j{x}_i{A}_{ij}}}}\mathrm{with} A_{ij}=\frac{{\left[1+{\left({\eta }_i/{\eta }_j\right)}^{0.5}{\left(M_j/M_i\right)}^{0.25}\right]}^2}{\sqrt{8\left(1+\left(M_i/M_j\right)\right)}}.$$

(9)

Here X denote molar fractions, and M the molar mass of the species i.

Figure 1 compares exemplarily the thermal conductivity of a H₂-N₂ mixture, calculated with Wilke’s model and the default mass weighted summation model available in OpenFOAM-6 against experimental data 12]. It is evident that the thermal conductivity is significantly underpredicted if the default model is employed.

To model multi-component diffusion, an effective binary diffusion model is employed. The pressure and temperature dependent species’ binary diffusion coefficients are calculated according to Fuller, Schettler, and Giddings. The effective binary diffusivity of the species ‘i’ in the mixture $D_{i,m}^{}$ is calculated according to [11]:

$$D_{i,m}^{}=\frac{1-x_i}{{\displaystyle \sum _{k\ne i}\frac{x_k}{D_{k,i}}}}$$

(10)

In containment typical flows, molecular transport becomes significant if a stagnant fluid volume is present (e.g., dead-end compartments or an atmospheric stratification) or when a low-Re wall treatment is employed.

2.2.3. Wall Treatment and Condensation

Containment structures are heat sinks driving the buoyant flows and transport processes. Considering the large difference in scale between the building and the flow boundary layers, a numerically efficient model for the near wall transport processes needs to be implemented in containmentFOAM. To be applicable on both fine meshes used in validation and coarse meshes used for an actual containment application, the underlying formulation is scalable and independent of the dimensionless wall distance y⁺. Thus, the wall fluxes are modeled by continuous dimensionless velocity, temperature and, in the case of wall condensation, by a species boundary layer profile. For velocity, the 4th order Spalding profile [13] is used. The thermal and species boundary layer is modelled using the Kader profiles [14]. As described by Menter [15], a zero gradient boundary condition for the turbulent kinetic energy is used and the eddy frequency is described continuously by blending the analytical solutions for the viscous and turbulent sublayers.

Wall condensation is modelled by means of a diffusion-layer model, which allows computation of the local condensation rate based on cell values (control volumes adjacent to a condensing interface) and does not require the definition of bulk parameters like correlation-based models. This approach assumes that the species boundary layer in the gas phase is the dominant transport resistance. Thus, the liquid condensate phase can be omitted. This implies that the boundary layer, including the viscous sublayer, is resolved (y⁺ ~ 1) or modelled properly to account for all relevant effects, such as wall normal mass transfer and near-wall buoyancy. The condensation mass flux per unit area ${\dot{m^{″}}}_{cond}$ can be described as:

$${\dot{m}}_{cond}^{″}=\frac{1}{1-Y_{H_{2}O,w}}\left({\overrightarrow{J}}_{H_{2}O,w}\cdot {\overrightarrow{n}}_w\right)$$

(11)

where ${\overrightarrow{J}}_{H_{2}O,w}\cdot {\overrightarrow{n}}_w$ is the wall normal diffusive flux of steam and $Y_{H_{2}O,w}$ is the steam mass fraction at the wall, determined at saturation condition.

In containmentFOAM, the condensation rate is treated as an additional outgoing flux through the wall face of a cell. Similarly, all transported quantities associated with the condensed mass are removed from the domain directly through the advection term. For further details on the model, its implementation and validation against the SETCOM (Separate Effect Tests for Condensation Modeling) forced convection database are summarized in [16]. SETCOM is a flow channel with square cross-section (edge length of 0.45 m). One side of the 6 m long channel consists of a cooled wall, while the other walls are adiabatic. Figure 2 exemplarily summarizes the wall fluxes and boundary layer profiles obtained by the implemented wall treatment for four different meshes with a boundary layer resolution of y⁺ ~ 1(resolved) to y⁺ ~ 80 (modeled) for the SETCOM test C_(U_in ~ 3 m/s, T_in ~ 70 °C, X_H2O,in ~ 30 vol.%, T_wall ~ 10 °C). As expected, the implemented wall functions provide a grid-independent solution of the wall shear stress, the sensible wall heat flux and the condensation rate once the boundary layers are developed (~1 m downstream the inlet). Comparing the dimensionless boundary layer profiles extracted at channel length of x = 5 m, a consistent representation can be determined for the different boundary layer resolutions. Compared to the measured total (i.e., sensible and latent) wall heat fluxes, again a reasonable and mesh insensitive representation by the model can be concluded for forced convection heat and mass transfer.

Ongoing work is directed towards a consistent wall treatment, with the ability to model mixed convection heat and mass transfer (Kelm et al. [17]). Bulk condensation is modeled as a single-phase phenomenon, while fog density ${\rho }_{fog}$ (i.e., fog mass per unit volume) is transported as a passive scalar using a drift flux approach.

$$\frac{\partial \left({\rho }_{fog}\right)}{\partial t}+\nabla \cdot \left[{\rho }_{fog}\left(U+U_{drift}\right)\right]=\nabla \cdot \left[\left(D_{Brownian}+\frac{{\nu }_t}{S{c}_t}\right)\nabla {\rho }_{fog}\right]+{\dot{S}}_{fog}^{‴}.$$

(12)

The drift velocity $U_{drift}$ includes effects of gravity and drag, while Brownian and turbulent diffusion are considered in the diffusive term.

The condensation/evaporation rate ${\dot{S}}_{fog,m}^{‴}$ is determined by the deviation from saturation condition, the available amount of steam and fog as proposed by Vyskocil et al. [18]:

$${\dot{S}}_{fog,m}^{‴}=\frac{{\dot{S}}_{fog,e}^{‴}}{h_{lat}},$$

(13)

The source term in the energy equation ${\dot{S}}_{fog,e}^{‴}$ is formulated as:

$${\dot{S}}_{fog,e}^{‴}=\frac{1}{\mathsf{\Delta }t}\left\{\begin{array}{c}\min \left[\rho c_p\left(T_{sat}(p_{steam})-T\right),\hspace{0.33em}\rho Y_{steam}{h}_{lat}\right] \mathrm{if} \mathrm{condensing}: T_{sat}(p_{steam})>T \\ \max \left[\rho c_p\left(T_{sat}(p_{steam})-T\right),\hspace{0.33em}-{\rho }_{fog}{h}_{lat}\right] \mathrm{if} \mathrm{evaporating}: T_{sat}(p_{steam})<T\end{array}\right..$$

(14)

Here, $\mathsf{\Delta }t$ is a time scale to reach saturation and $h_{lat}$ the latent enthalpy.

Ongoing work addresses the tighter integration of fog transport by means of a mixture model into the solver and enabling interaction with other physical models, such as radiation or aerosol transport.

2.2.4. Radiative Heat Transfer

Surface to surface and gas radiation heat transport in humid atmospheres impacts structure and gas temperatures [6,19], even at small temperature differences (<50 K). The radiation heat transfer changes the gas temperature field, which affects the buoyancy force and thereby the mixing process. Issues related to the thermal stratification and pressurization have been identified within the ERCOSAM project (e.g., [20]). Furthermore, condensation rates are affected by a change in humidity due to changes in the gas temperature and saturation conditions due to changes in the structure temperature [21]. Thus, gas radiation affects the water-steam balance.

Modeling thermal radiation transport in a containment is challenging due the complex geometry, the variation of absorption coefficient with wavelength, and the numerical error of radiative solvers. In OpenFOAM^®, there are two in-built radiation models, P1 and fvDOM. The P1 model is the lowest-order of spherical harmonics method, assuming the radiation intensity is near-isotropic. In other words, the P1 model cannot resolve the directional dependency of radiative heat transfer. Moreover, the P1 model is only applicable for an optically thick media. The finite volume discrete ordinate method (fvDOM), however, solves angular radiation intensity for a set of quadrature directions. Nevertheless, it is known to suffer from false scattering and ray effects [22] and is also computationally expensive [23].

For this reason, the Radiative Energy Absorption Distribution (READ) Monte Carlo method [24] with the statistical narrow band correlated-k (SNBCK) non-gray gas model [25] has been implemented into the containmentFOAM framework to solve the radiative heat transfer equation (RTE) accurately. The Monte Carlo (MC) method is a statistical method and simulates radiative heat transfer by tracking many photon bundles (histories) through arbitrary complex 3D geometric configurations; undoubtedly, it is computationally expensive. In this paper, the photon tracking module is based on the Lagrangian library in OpenFOAM^® using the CFD mesh. However, if many photons are tracked on mesh partitioned along multiple processors, the MPI communication time becomes intolerably high, since the photon information is transferred to the associated processor when passing a processor boundary patch. To overcome this issue in containmentFOAM, MC radiation transport is computed on a global mesh, while the number of photon histories is distributed among the processors. Consequently, there is no communication among the processors during the photon tracking and no CPU imbalance between processors.

Before coupling the MC method to the CFD solution, the MC RTE solver with the SNBCK database is verified against analytical 2D test cases by Goutiere et al. [26]. The basic configuration is a 2D rectangular box (1 m width, 0.5 m height) with a wall temperature of T_w = 0 K, i.e., the black walls are not emitting themselves (compare Figure 3 top left). Different gas mixtures (absorptivity) and temperature fields are defined and the volumetric radiative source terms as well as wall heat fluxes, obtained from the SNBCK model are compared against a reference solution. In Figure 3, the results for ‘case 3’ with a homogeneous mixture (x_H2O = 20 vol.%, x_N2 = 80 vol.%) and gas temperature (T = 1000 K) are given. It is noted that the temperature and flow field is frozen and only RTE equation is solved. It is evident that the volumetric heat sources and wall heat fluxes predicted by the new Monte Carlo solver agree reasonably well with the reference results. The visible fluctuations (<±3%) are to the Monte Carlo method and the only noticeable deviation is a slight underprediction in the wall heat flux.

The limitation of the MC method is the statistical error, which can be reduced either by increasing the number of photon histories or improving the algorithm. The latter gains importance in containment flows because the temperature difference is relatively small compared e.g., to combustion applications.

In containmentFOAM, the algorithm is improved in two ways: First, by importance sampling [27], the photons are distributed according to the cell temperature, using the inherent bias that the hotter cells contribute more to the radiative source term. Therefore, it allows reducing the number of photons without losing accuracy. Second, the shifted forward MC method, proposed in [28], or the emission reciprocity method [29] are applied. These advanced algorithms significantly improve the efficiency in quasi-isothermal media without increasing the number of photon histories. Tests revealed that the number of photon histories and thus the computation time can be reduced by 1 to 2 orders of magnitude using the advanced algorithms.

Ongoing work aims to validate the model for accidental conditions, i.e., low temperatures (<600 K), high vapor fraction (>70 vol.%), and moderate pressures (<10 bar) based on new experiments, conducted in the OECD/NEA HYMERES-2 project or the German THAI project. To further improve the MC numerical efficiency, a coarsening algorithm for the radiation grid is considered, considering that the gradients in the radiative intensity are less pronounced than in the fluid flow.

2.2.5. Aerosol Transport and Decay Heat Distribution

Besides noble gases, aerosols are the main carriers of the fission products in the containment atmosphere. Thus, their transport and deposition determine the decay heat distribution in the containment atmosphere. System code analysis revealed a mean volumetric source in the order of 100 W/m³ gas. It is expected that the airborne and deposited heat sources will affect the buoyancy driven flows and transport processes and thus need to be considered in the prediction of the containment atmosphere mixing processes.

As an initial step a Lagrangian approach was adopted in containmentFOAM. It models aerosols as dispersed Lagrangian particle clouds using a loose coupling with the Eulerian continuous gas phase ([3,30]). The model considers all relevant forces, such as drag, thermophoretic force, and gravity, while buoyancy, pressure gradient, virtual mass, and particle-particle interactions are neglected. Due to the large variation of aerosol particle size from 0.01 to 10 μm [31], drag and lift forces on the particles are modeled by incorporating the Cunningham correction factor, which is important for particles below 1 μm. Turbulent dispersion is considered using the Continuous Random Walk (CRW) model based on normalized Langevin equation [32]. It was implemented in containmentFOAM to overcome deficiencies in the stochastic Discrete Random Walk (DRW) model, available in OpenFOAM-6.

A first distinct validation of the CRW implementation is conducted against experimental data of particle deposition in turbulent pipe flow by Liu and Agrawal [33]. All particles touching the wall are considered as deposited in this simulation. As shown in Figure 4, the dimensionless particle deposition velocity is predicted with reasonable accuracy in forced convection flows. For the particles in the diffusion-impaction regime, the model slightly under-predicted the deposition, which is currently under investigation.

Besides the Liu-Agarwal experiments, an anisothermal turbulent particle-laden pipe flow by Dumaz et al. [34] was reproduced with good accuracy. The decay heat of the particles, modeled as a volumetric heat source, was revealed to have a significant effect on the fluid flow, and the particle deposition in a small scale differentially heated cavity (Rayleigh number Ra ~ 10⁹). However, the analysis also revealed that the Euler-Lagrangian approach is too expensive for technical scale applications.

Consequently, a Eulerian treatment of the aerosol phase with the mixture model [35] is currently under investigation to provide a computationally more efficient approach. In addition to advection due to the bulk motion of the gas phase, Brownian diffusion is the dominant transport mechanism in the sub-micrometer aerosols, whereas particle inertia results in a relative velocity with the gas-phase for the larger aerosols (>1 μm) and makes them drift from the flow streamlines. Subsequently, both particle drift and diffusion are accounted for in aerosol transport based on the work of [36]. In the late phase of a severe accident, stagnant zones may form (e.g., stable atmospheric stratification). In these regions, aerosol deposition is only influenced by gravity. Aerosols of different sizes settle under gravity at different rates; therefore, predicting the evolution of aerosol size distribution becomes necessary. The size distribution is affected by processes such as agglomeration, de-agglomeration, and phase-change phenomena (evaporation or condensation) on aerosols. The evolution of the particle size distribution is computed by solving the population balance equation (PBE) with the sectional modeling approach [37]. Some of the substances (CsI, CsOH) present in the radioactive aerosols are hygroscopic and they grow by absorbing from water in the humid containment atmosphere, thereby affecting the size distribution. This steam condensation on aerosols involves both heat and mass transfer and hence affects gas distribution. Since a typical aerosol present in a containment atmosphere is a mixture of multiple species, the hygroscopicity for a multicomponent material is calculated using the works of [38].

2.2.6. Safety Systems

Several possibilities exist to model the operational behavior of safety systems, such as porous media or code coupling. To model the efficiency of a PAR, which is determined by its heat balance, driving the buoyant chimney flow and its interaction with the surface chemistry at the catalysts, the mechanistic PAR model REKODIREKT (RD) [39], developed at Forschungszentrum Jülich, has been coupled to containmentFOAM. RD models the entire PAR based on a detailed representation of a sub-channel between two catalyst elements and a chimney model to predict the buoyancy driven flow through the PAR. It has been developed based on the detailed REKO database and enables the prediction of the H₂ and CO conversion rate under all relevant conditions.

Based on previous experience [40], a domain decomposition approach is used. The coupling occurs on two boundaries: At the ‘PARinlet’ interface, the area averaged gas composition, and temperature as well as the system pressure is evaluated in containmentFOAM and transferred to RD (see also Figure 5). The latter computes the conversion rate and delivers gas composition and temperature at the ‘PARoutlet’ interface as well as the buoyant mass flow rate through the PAR. The explicit transient coupling is realized by means of a Python wrapper. Besides being the interface between C++ (OpenFOAM-6) and Fortran 90 (RD), the python wrapper organizes a consistent initialization of RD, and takes care of the data logistics; such as intermediate data storage, data conversion (fields <> scalars), and backup and restart files. In principle, an arbitrary number and type of PARs can be considered. Due to the fast runtime of RD, both codes are executed serially. In parallel runs, RD is only executed on the master exchanges data with the different partitions. The necessary user input for PAR modeling is the definition of the coupling interface, the PAR geometry, and logical criteria that define the activation of the PAR, such as a time point (typical in validation experiments) or a minimum reactants concentration.

Besides PARs, there are other safety relevant technical components to be considered in a containment-scale application, e.g., the opening of pressure relief flaps and doors is modeled on basis of an arbitrary mesh interface (AMI), which opens depending on a differential pressure condition. It enables to consider opening of a new flow path during the accident sequence. The modeling of heat exchangers used e.g., within the containment cooling system, is in progress.

3. Verification and Validation

3.1. Overview

Over the last ten years, significant efforts have been taken to develop and validate a baseline set of models necessary to predict accidental flows in a nuclear reactor containment and connected buildings as well as system behavior e.g., passive autocatalytic recombiners (PAR) on basis of the commercial CFD package ANSYS CFX (up to v17.2) [1]. The defined set of baseline models were validated systematically on experiments conducted in different facilities and at different scales. Technical scale experiments on the establishment of an atmospheric stratification, its stability, and remobilization under different conditions, which were conducted within the OECD/NEA-SETH and SETH-2 [41] projects in the MISTRA (CEA, France) and PANDA (PSI, Switzerland) facilities, were employed for validation. A similar phenomenology was addressed at the well instrumented small-scale MiniPanda facility (ETHZ, Switzerland [42]) and the OECD/NEA PANDA benchmark [43]. Within the recently concluded OECD/NEA-HYMERES and ongoing HYMERES-2 projects, this database was extended by test configurations, where a diffuse flow behind a flow obstruction is created and leads to erosion of the stratification [44,45]. Natural convection flows driven by wall heat transfer were investigated in the MISTRA and THAI (Becker Technologies, Germany) facilities [46]. Recently, the effect of thermal radiation heat transfer on these mixing processes was elaborated [19,21]. Detailed models were developed to analyze PAR performance [39,40]. An excerpt of the validation matrix is given in

Table 2. Validation matrix—comparative assessment of experiments at different scales (excerpt).

Phenomenon	Separate Effect Test	Coupled and Integral Effect Test
Basic Transport Phenomena:
• Stratification built-up	Garage	THAI: HM2, TH21 PANDA: SETH, PE0 TOSQAN: ISP-47 MISTRA: MERCO-0, MM3
• Mixed/free convection	SETCOM / CONAN
• Wall /bulk condensation	SETCOM / CONAN
• Gas radiation heat transfer	Goutiere et al., Cassol et al.
• Buoyancy turbulence	FLAME
Stratification re-mobilization:		THAI: TH22, TH24, TH26 PANDA: SETH-2, HYMERES MISTRA: NATHCO, MERCO-2
• Jet/momentum driven	MiniPanda
• Plume driven
• Diffuse flows
• Wall bounded natural circulation
• Thermal plumes
Safety systems: • PAR operation	REKO	THAI: HR
Aerosol transport: • Drag, lift and thermophoretic force • Turbulent dispersion	Liu-Agrawal STORM, DIANA, TUBA	THAI: AW, HD Phebus KAEVER
Scenario (scaled):		THAI: ISP-47 MISTRA: MERCO PANDA: PE BMC: VANAM M3/M4

.

To build on this comprehensive experience, the baseline models were transferred into containmentFOAM, rather than utilizing similar ones available in OpenFOAM^®. Likewise, validation runs can also consider a code-to-code comparison with existing results. The validation follows the established Best Practice Guidelines (BPG) [47] to minimize the numerical errors.

Both the containmentFOAM solver and model library have already undergone a first fundamental validation (e.g., [3,5,6,16]), is continuously further validated in the frame of ongoing benchmarks such as within the national THAI or the OECD/NEA HYMERES-2 projects and existing application-oriented validation cases. Along with the systematic and continuous validation process, a methodology for uncertainty assessment is developed and implemented [48,49]. In the following section, an application-oriented validation case, focusing at H₂ mixing and mitigation by means of PARs, is summarized.

3.2. Validation Case: PAR Operation under Accidental Conditions

A comprehensive test series on the performance of passive auto-catalytic recombiners (PAR) has been performed within the OECD/NEA THAI projects for different PAR types. For this validation run, the test THAI HR12, conducted in the frame of the OECD/NEA THAI-1 project [50] is considered. Its test configuration is illustrated in Figure 5 (left). A Framatome PAR FR380T is mounted at the outer side of the inner cylinder. The vessel is preconditioned at accident relevant conditions i.e., an initial gas and structure temperature of T₀ ~ 150 °C at P₀ ~ 3 bar and an initial steam concentration of X_H2O,0 ~ 60 vol.%. H₂ is injected with a mass flow rate of up to 0.45 g/s via a ring feed line below the PAR in two phases. The modeling challenge of this application-oriented test case is the strong effect of gas radiation heat transfer and wall condensation on the vessel pressure and thus also recombination rate. Besides, there is an initially low oxygen concentration of X_O2,0 ~ 8 vol.%. After the first injection phase, the PAR has consumed most of the available oxygen and runs into oxygen starvation at the beginning of the second H₂ injection phase, i.e., the conversion rate and efficiency is impaired.

A fully block-structured hexahedral mesh is employed (cf. Figure 5. THAI HR12 test: Scenario (left), fully hexahedral block-structured grid (right) Figure 5 right), and no symmetries are considered.

The grid is refined in the free shear layers of the H₂ plume, close to the walls and in the vicinity of the PAR. The flow in the vessel is solely driven by the buoyancy of the injected H₂ and the hot exhaust gas plume of the PAR and is quite unsteady. Therefore, the remaining volume of the facility is meshed with almost homogeneous cell size. The block topology and grid spacing is based on previous validation experience [40,41] and was proven to allow for a grid-independent prediction of the gas transport and mixing. The grid quality was optimized according to the BPG [49], in terms of the minimal face angle and the maximum volume and aspect ratios. Since the THAI inner cylinder and vessel walls represent a significant heat capacity, and there are visible heat losses that lead to a stagnant zone in the vessel sump, a conjugate heat transfer approach is employed. Fluid and vessel domain are initialized by the measured vertical temperature profile, obtained by circumferential averaging of the available sensors. The gas composition in the test was nearly homogeneous and is initialized by the average measured concentrations. The individual nozzles of the ring feed line are not resolved but modelled as an inlet boundary patch with transient definition of the measured gas mass flow rate, composition, and temperature. This modeling decision is expected to be applicable as the jet to plume transition length is in the range of a few centimeters and thus the conservation of the momentum of the individual jets is not important for further transport and mixing. The turbulent kinetic energy and eddy frequency are calculated assuming a low turbulent intensity of I = 1% and an eddy length scale of L = mm (diameter of the injection line holes). The injection line is assumed to be adiabatic. All walls are defined as no-slip walls, using the blended wall function approach (compare Section 2.2.3). The emissivity of the vessel walls is defined as ε = 0.6 (oxidized steel). The vessel oil heating system and the isolated parts are not resolved and simply modelled by an effective heat transfer coefficient and secondary side temperature (either the environment T_env = 25 °C or average oil temperature). The coupling and data exchange with RD occurs on two boundaries (see Figure 5 right): At the ‘PARinlet’ boundary, the area averaged gas composition, and temperature as well as the system pressure is evaluated in containmentFOAM and transferred to RD, which computes the gas composition and temperature at the ‘PARoutlet’ as well as the buoyant mass flow rate through the PAR. The model and numerical settings correspond to the baseline approach defined for the validation of containmentFOAM (refer to Table 1 for the models and [3] for the numerical methods). The validation case was run at time step of Δt = 0.1 s, which corresponds to a CFL number less than seven, and took ~6 days (4531 CPUh) on 32 Intel^® Xeon^® E5-2637v4 CPU @ 3.50GHz. In this, the runtime of RD accounted for less than 1 h.

In this experiment, the vessel atmosphere is well mixed for most of the transient. Consequently, for the sake of brevity of this paper, primarily the thermo-fluid dynamic boundary conditions for PAR operation (gas composition, temperature, and vessel pressure) as well as the REKODIREKT, output values are compared against the experimental data. The reactants (H₂ and O₂) concentration at the PAR inlet is predicted to be of good consistency for the first injection phase, while the PAR inlet concentration of H₂ is visibly under predicted in the second phase (see Figure 6 left), but still reveals a comparable reduction between inlet and outlet concentration (ΔX_H2 = X_H2,in–X_H2,out), as indicated by the arrows. Considering the PAR inlet and exhaust gas temperature (Figure 6, middle), a good agreement with a slight overestimation of the second peak temperature can be found for the complete transient. This implies that the reaction/heat release rate as well as the thermal inertia of the PAR (indicated by means of the temperature of the PAR housing T_box, where no experimental measurement was available) is still modelled properly, despite the too low inlet H₂ concentration in the second phase. The resulting density differences create a buoyancy driven flow through the chimney, which is also predicted in reasonable consistency with the experiment (see Figure 6, right), indicating that the pressure losses in the PAR are modelled correctly. It should be noted that the uncertainty of the flow rate measurement (approx. 10%) is of high impact and propagates further into the calculation of the reaction rate.

On basis of the consistently predicted boundary conditions, the major PAR performance parameter, the H₂ recombination rate is evaluated and validated against the measurements. (Figure 7, left). The simulation reveals a well comparable trend, both qualitatively and quantitatively. To highlight the effect of oxygen starvation qualitatively, the recombination efficiency (X_H2,out/X_H2,in) is compared in Figure 7 (middle). One can observe the drop due to increasing lack of oxygen in the second injection phase, which is predicted consistently with the experimental observation.

The observed small quantitative deviation between experiment and simulation, comprises of various sources (e.g., differences between predicted PAR flow rate, pressure, temperature etc.). It is certainly as acceptable because the recombination rate is the relevant performance criterion and in correlation-based models, the efficiency is typically a user defined parameter.

The global heat and mass balance especially governed by gas radiation and condensation heat transfer is well predicted as depicted in Figure 7 (right) on the basis of the vessel pressure.

Concluding, the application-oriented validation of containmentFOAM regarding H₂ mitigation yielded no systematic lacks. Future work will address the parallel recombination of H₂ and CO, which has recently been investigated at the THAI facility.

4. Dissemination

It is intended to release the baseline version of containmentFOAM under the GNU Public License (GPL) v3 for beta testing (upon reasonable request) once fundamental quality assurance is completed in the near term. New model improvements and extensions shall follow as soon as the corresponding research has been published. The focus on open-source software development aims at enabling open collaborations as well as the free dissemination to and application of the research outcomes by a broad community.

In this spirit, containmentFOAM is further developed to become an open basis to predict containment flows, pressurization, H₂ and aerosol behavior, and associated safety measures. In the mid-term, containmentFOAM shall become a part of the multiphysics simulation platform for education and research in nuclear applications, developed within the open-source initiative launched under the aegis of the IAEA [51]. On the other side, it shall establish an open link to related non-nuclear research fields e.g., safety assessment of industrial processes or the applications of hydrogen as an energy carrier.

In contrast to commercial multi-purpose CFD packages, containmentFOAM is rather an expert tool that demands a skilled user. Application of CFD to containment flows requires the combination of multiple models to represent the underlying physics and safety systems. The single models are typically developed within individual PhD research projects and comprise of comprehensive expert know-how (e.g., in the definition of numerical parameters or the combination of available sub-models). Their definition is conducted primarily within individual text files, so-called dictionaries (see Figure 8, left), but may also involve links to other dictionaries. It is obvious that there is a certain potential for user dependent errors, resulting from a false or incomplete model definition or inconsistencies among the dictionaries, which must be avoided to ensure trustworthy and reproducible results. Classically, this can be done e.g., by comprehensive documentation (document or Wiki), commented dictionaries, and tutorial cases. Nevertheless, such documentation cannot consider all details and dependencies and is often lagging the code development. For these reasons, the containmentFOAM package relies on a guided case setup process (see Figure 8, right).

It consists of a Java GUI, which guides the user through the case setup process. It is based on comprehensive JSON-based templates, and rules that allow on the one hand a flexible extension and on the other hand to define logical dependencies (e.g., if there is no steam, the wall condensation model is not available) to ensure consistent model definition in the various dictionaries. Furthermore, tooltips allow direct access to explanations and avoid searching in a comprehensive documentation. The definition of the JSON files includes the ‘best practices’ of the original model developer, is based on his/her experience, and integrated along with the model development and V&V. The numerical methods and default settings specified for validation runs in [3] are the basis for the ‘simulation control’ template.

By this means, but also a continuous documentation of the modeling basis and application experience, the comparability of results obtained by different users and correct application of the models shall be maintained, even if there is no direct contact between developers and users. Furthermore, a standardization is the only mean to ensure a certain quality of the simulation results and enable identification of model deficiencies.

Currently, a flexible monitor is developed for live analysis of simulation progress in terms of convergence, performance, and monitor variable histories. Different filters, such as floating average or FFT are available to analyze the comprehensive solver output during run time and identify potential issues at an early stage.

5. Summary and Outlook

The dedicated CFD simulation tool containmentFOAM for the analysis of transport and mixing phenomena inside the containment and connected buildings of a nuclear reactor is developed based on the open source CFD toolbox OpenFOAM^® (currently version 6). It comprises of a multi-species solver integrated in a conjugate heat transfer framework as well as a comprehensive model library to represent the containment phenomenology. The latter consists of the baseline model for containment analysis, developed and validated at Forschungszentrum Jülich since more than 10 years ago, in addition to ongoing new developments e.g., regarding thermal radiation, or aerosol transport in the containment atmosphere, or a wall treatment for mixed convection.

containmentFOAM has undergone a fundamental validation against separate effect tests and is now continuously validated against application-oriented coupled and integral effect tests. In this context, a framework for quantification of the propagation of input uncertainties was developed and integrated in collaboration with the Universität der Bundeswehr Munich. An application relevant validation case in the THAI facility was discussed and it demonstrated the numerical efficiency and predictiveness for an application to H₂ mixing and mitigation inside the containment. In comparison with commercial solvers used before, containmentFOAM yields a comparable and partly even better performance. Besides the systematic validation, prototypic full-scale applications will be conducted to investigate the robustness of solver, numerical methods, and models on complex geometries and unstructured meshes.

Publication of containmentFOAM is intended in the near term as a beta version and in the long term as part of the open source multiphysics platform, developed under the aegis of the IAEA. As a major measure to minimize potential user errors, a guided user interface is developed. It enables a guided and standardized case setup under consideration of developers’ best practices and established CFD Best Practice Guidelines.