Gas mixing caused by interacting heat sources. Part ii: Modelling

Abstract The analyses of hypothetical accident scenario in various containment designs include the evaluation of the effectiveness of Passive Autocatalytic Recombiners (PARs), which are installed to mitigate the hydrogen risk, or are considered by utilities for retro-fitting the plants to implement revised regulatory requirements motivated by the Fukushima accident. The assessment requires validated computational tools. To build confidence in the codes, their capability must also be assessed against separate effect tests addressing specific phenomena under conditions that are representative of those expected in an accident scenario in a real plant. Under the conditions of a severe accident, the large amount of hydrogen release requires the installation of several units, and therefore the study of their interaction is necessary for the optimal layout. Within the OECD/NEA HYMERES project, four experiments have been performed in the large-scale PANDA facility, addressing the thermal effect of one or two PARs, which was simulated by means of heaters (Part I). The tests were of increasing complexity, to permit a step-by-step validation of the codes. These experiments have been analyzed with the GOTHIC code. In general, the GOTHIC code in conjunction with a coarse mesh could predict the mixing process reasonably well. The only notable discrepancy with the experiments was the overprediction of the mixing rate when the lower heater was activated. The analyses revealed that this was due to the incorrect representation of the interaction of the flow out of the lower heater with the obstruction, and better results could be obtained by refining the mesh around the heaters. It was concluded that a coarse mesh is adequate for representing the flow produced by multiple heat sources in an open geometry, but detailed models are generally required if the flows could interact with structures.


Introduction
The release of hydrogen in the containment that can occur in a Light Water Reactor (LWR) during a severe accident is a serious safety issue, as recently demonstrated by the course of the accident in three reactors at the Fukushima site (OECD/NEA, 2016). The hydrogen risk has been studied for a long time, and it was recognized that the risk to the integrity of the containment would be enhanced by stratification, due to the possibility of formation of explosive mixtures in specific compartments or zones. In order to mitigate such a dangerous hydrogen accumulation, in various reactor designs the containment is equipped with Passive Autocatalytic Recombiners (PARs), which produce a significant amount of heat. This heat, transferred to the gas, induces buoyancy driven flows through the PARs and in the containment. For the evaluation of the effectiveness of these devices, appropriate computational tools are needed, which account for the 3-D distribution of gases in large volumes (Royl et al., 2000;Xiong, 2009;Reinecke et al., 2010;Sonnenkalb et al., 2015;Visser et al., 2015;Lopez-Alonso et al., 2017;Saghafi et al., 2017;Papini et al., 2018;Choi et al., 2018). Several experimental programmes addressed the behaviour and the effect of single PARs, either using hydrogen and real units (Gupta, 2015) or their thermal effects on mixing using helium and heaters. For instance, in the SETH-2 and ERCOSAM projects (NEA/CSNI, 2012;Paladino et al., 2016) experiments have been conducted where a heat source was activated after the stratification (of both helium and steam) had been already established in the upper part of one vessel. Under the conditions of a severe accident, however, the large amount of hydrogen release requires the installation of several units, and therefore the study of effect of the position of the PARs on their effectiveness (Meynet et al., 2008;Park and Bae, 2014;Hoyes and Ivings, 2016;Lee et al., 2016). Moreover, the evaluation of their interaction is necessary for the optimal layout (Deng and Cao, 2008;Huang et al., 2011;Bachellerie et al., 2003;Serrano et al., 2016;Papini et al., 2018). Although an extensive database exists for the assessment of the overall modeling capability of the codes, very few tests at large scale have been carried out (e.g., Kanzleiter, 1997) with real PARs. Similarly, only few data exist for the basic validation of advanced containment and CFD codes in relation to specific phenomena, such as the interaction between units (e.g., in the MISTRA facility, NEA/ CSNI, 2018), These tests featured mostly units at the same elevation or slightly staggered with respect to vertical elevation, and activated at the same time. Therefore, within the HYMERES project (NEA/CSNI, 2018), one of the aims of the new tests was the investigation of the effect on mixing of the interaction between two heaters located at notably different elevations, and activated at different times to reflect the delay in the start of PARs located in regions where the light gas concentration increase is slower. In the new tests , a few additional effects are addressed, which are playing a role in real containment analysis, notably the effect of wall condensation. Indeed, previous experiments have been conducted where the condensation at the walls was either prevented (NEA/CSNI, 2012) or was very small, as a result of the heat-up of the metallic structures during the earlier period of the transient (Paladino et al., 2016). Under those conditions, designed for simplifying the basic validation of codes, independently of the existence of nearly stagnant conditions or the use of a continuous release below the stratification front, the mixing above the heater inlet was complete, but the helium remained confined in that region, with the stratification front moving very slowly downwards. It was concluded that a heat source is not capable to produce global mixing. This finding is in accordance with previous analyses showing that a PAR (heat source) at high elevation can only produce strong convective currents above the inlet (Dabbene, 2011). It was also recognized that the limited mobility of helium in these tests was mainly due to the absence of condensation , which is not representative of the actual conditions in Containment. Therefore, the new experiments aimed to provide a database for studying the capability of the codes to simulate the effect of heat sources under conditions closer to the prototypical ones. These tests thus featured large condensation rates, comparable with those produced by condensation on the containment walls, and the activation of the heat sources at the time the helium reaches the entrance of the heater box, to represent the condition that triggers the start of the chemical reaction in a PAR. The scoping calculations  permitted to identify a configuration for the tests that could reproduce the main features of the thermal interaction of two PARs' in a real plant, although the consumption of the light gas (hydrogen in a real PAR) cannot be represented. With respect to the condenser tubes, all the 3-D effects of condensation over the entire containment wall could not be fully represented by the vertical condenser located in one sector of the vessel. However, the condensation process produced by the three tubes was adequate to reproduce extent and time scale of the downward transport of the light gas to the entrance of the lower heater. Moreover, additional planning calculations (not reported in Andreani et al. (2018)) with heat exchange area distributed over four sectors of the vessel periphery showed that no substantial difference in the gas distribution would be produced by condensation on one side of the vessel in comparison with a more homogeneous partition.
Since the tests were aimed to produce data for the validation of the codes, the obstruction plate close to the exit of the lower heater (not considered in the planning calculations) was also installed because the analyses (especially pre-test) of previous tests (Andreani et al., 2019) had shown that codes had difficulties to simulate experiments with this configuration. The obstruction plate, as well as the injection pipe, are not necessarily intended to represent a specific object in the containment, but, since they are close to the lower heater outlet, represent an additional challenge to the codes. In principle, the crane would constitute a similar obstruction for certain layouts of the PARs, as it can be inferred from the set-up of the analyses of Lopez-Alonso et al. (2017) and Papini et al. (2018), and the injection pipe could represent an obstacle similar to the steam generator external wall.
The paper addresses the validation of the GOTHIC code (EPRI, 2016) using the results of the HYMERES experiments. The code, associated with the use of a coarse-mesh, was successfully validated against test with heater in previous projects (Andreani and Kelm, 2012;Andreani et al., 2015), and in this paper a very similar modelling approach is used for analysing the more challenging conditions of the HYMERES tests. In particular, as discussed above, the two main effects addressed in the tests, namely the effect of condensation on the mobilization of light gas and the interaction between the two heat sources are looked at in detail. This will permit to draw some conclusions on the capability of the code to simulate the phenomena associated with the operation of two heaters in presence of condensation driven natural convection.

The PANDA facility and the experiments
The PANDA facility and the experiments are described in Part I of the paper , and therefore only a few main features are given here. The PANDA facility and the configuration of the tests are shown in Fig. 1, where the condenser (three tubes) and the heaters are also included in the schematic. The heating elements are plates, installed in the lower part of the housing, which has an opening at the bottom (inlet), and one opening (outlet) on the upper part of the side facing the axis of the vessel. A horizontal obstruction plate is installed at 1 m above the exit of the injection pipe.
Four experiments have been conducted in the two upper vessels, with the same initial conditions (a mixture of air and steam -54% and 46%, respectively-, with the thermal structures in approximately thermal equilibrium with the fluid at around 115 • C). In all tests, superheated steam at about 150 • C (about 30 g/s) was injected for about 2000 s, helium (about 3.2 g/s) was injected in the time span between 500 s and 1000 s, and cold water flowed in the tubes during the entire duration of the tests. Only the time of activation of the heaters was varied, with the "zero" test being conducted without operating the heaters to provide data for validating the modelling of the condenser and its effect on the downwards transport of helium. Additional information is given by Paranjape et al. (2020). For the sake of reading convenience, the test matrix is also summarised in Table 1, and the boundary conditions (inlet mass flows and heaters powers) for test HP3_1_5 (with both heaters in operation) are also shown in Fig. 2. The figure also shows that the heat rates actually transferred to the fluid (calculated by the code) are lagging behind the electric powers due to the thermal response of the heater plates. For the other three tests, the two mass flows were practically the same, whereas the power input was absent (Test HP2_0_1), given to the upper heater only (Test HP2_1_2), or to the lower heater only (Test HP2_2_2).

The GOTHIC code
The thermal-hydraulics module of GOTHIC (EPRI, 2016) is based on a two-phase, multi-fluid formulation, and solves separate conservation equations for mass, momentum and energy for three fields: a multicomponent gas mixture, a continuous liquid, and droplets. In addition, species balances are solved for each component of the gas mixture. GOTHIC includes a full treatment of the momentum transport terms in multi-dimensional models, with optional models for turbulent shear, and for turbulent mass and energy diffusion. The options for turbulence are the mixing-length model and several variants of the k-ε model. In this work, the standard k-ε model has been used. The hydraulic model of GOTHIC is based on a network of computational volumes (one, two or three-dimensional) connected by flow paths. In contrast to standard CFD packages, in GOTHIC the subdivision of a volume into a multidimensional grid is based on orthogonal co-ordinates. The actual geometry of volumes with curved surfaces (e.g., a cylindrical vessel) can be represented, however, by blocking groups of cells. This feature was used to represent the PANDA vessel (see the figure displaying the nodalisation presented in the next section, where the blocked cells painted in dark grey appear in the subdivision of the volumes). The numerical solution of the transport equations is based on a semi-implicit method.
The method is first-order in time, whereas for the space-discretisation of the advection term both a first-order upwind and a bounded secondorder method are available. In this work, the bounded second-order method was used, which proved in past investigations to greatly reduce numerical diffusion on a coarse mesh in simulations addressing buoyancy-dominated flows. The thermal effect of structures is fully represented: 1-D or 2-D heat conduction equation is solved in the solid and built-in models for all heat transfer modes are available for the thermal coupling between structures and fluid. In particular, very advanced condensation model are included in the heat transfer package. For these analyses, the default model (DLM-FM) based on diffusion theory is adopted. Simplified wall-to-wall and wall-to-gas (in the cell adjacent to the wall) radiation models are also available, but no model (radiation in a participating medium) exists for radiation within the gas and from the bulk of the gas to the wall. In these analyses, no radiation is considered. The code has been extensively used for 3-D analyses of accident scenarios in Containments (e.g. Bocanegra et al., 2016). The code has built-in models for the simulation of recombiners, and has been used for optimizing the layout of the PARs in containment (Lopez-Alonso et al., 2017;Papini et al., 2018). Since helium was used instead of hydrogen, these models are not relevant for the tests. The code also has a model for heaters, but this does not account for the thermal inertia of the heating elements and the heat source cannot be distributed among several cells. Therefore, the choice was made to model the heater plates explicitly, using thermal conductors with prescribed, time-varying volumetric heat sources. Since the chimney effect produced by the tall case in real PARs cannot be produced by ad-hoc models (as it occurs with the PAR model implemented in GOTHIC), the heater height must be subdivided in a stack of cells in order to generate the appropriate gravity head, which, in turn, drives the natural convection through the heater.
The version of the code used in the present analysis is GOTHIC 8.2a (D170121), which is essentially the same as GOTHIC 8.2 (QA) (EPRI, 2016) with the addition of a few patches, which were afterwards implemented in the later versions of the code.

The model used for the analyses
The schematic of the model is presented in Fig. 3. The two interconnected vessels are represented with three subdivided volumes (1 s, 2  s and 3 s), connected by 3-D flow connectors. The two heaters are represented with two subdivided volumes (4 s and 5 s), connected to Vessel 1 by two pairs of 3-D connectors, representing inlet and outlets, respectively. Differently from previous analyses (Andreani and Kelm, 2012), the heaters are represented by separate volumes rather than being embedded in the mesh. This choice was motivated by convenience in the preparation of the input, due to the need to use for each heater only one blockage rather than a combination of several volume and/or face blockages. Moreover, due to the need to ensure that no spurious flow paths between the fluid inside the heater and the surrounding vessel result from the blocking of groups of cells with thin elements, the use of separate volumes eliminates the need for a tedious verification that this is actually achieved. Finally, the possibility to subdivide the heaters independently of the mesh of the surrounding volumes permit easier performance of mesh sensitivity studies. The user convenience of this approach (at least in the authors' opinion) in the pre-processing is not decisive if the model includes only one heater, but increases with the number of components. In Section 8, sensitivity studies will be presented with the use of the embedded mesh approach. The condenser (three tubes, as described in the Part I of the paper, where the experimental setup is illustrated) and the inlet section (three injection pipes) are also represented by subdivided volumes (7 s and 6 s, respectively), connected by flow paths. The upper collector and the vertical hose carrying the water outside the vessel are represented with two LP volumes (9 and 10), connected to the condensers through flow paths (one for each of the tubes). Finally, another LP volume (10) is used to store the water drained from the trays installed to collect the condensed water at the bottom of the tubes. This volume is connected to cells (blocked for liquid flows in vertical as well as transversal directions) in Vessel 1 where the water is extracted through three flow paths. The water flowing into the condenser tubes is fed to the individual injection pipes using three separate flow boundary conditions, and a pressure boundary condition closes the water circuit of the primary side of the condenser. The steam and helium injections are represented by two flow paths, which connect the gas sources to the bottom of the straight section of the injection pipe (see nodalisation in Fig. 4). Various thermal structures represent the vessel walls. The thermal effect of the structures is mostly represented using conductors assigned to the blockages: however, in the critical areas, namely the bottom and top headers, additional conductors are used to ensure the correct distribution of heat transfer areas in those regions. The heat transfer from the vessel to the water in the tubes is simulated using three separate spanned conductors, one for each tube (Fig. 3). The nodalisations on both sides of the tubes (Fig. 5) are matched in all directions, to ensure that the heat transfer areas of the subconductors on the primary side are correctly connected with the corresponding cells in the vessel. Additional conductors are transferring heat from the vessel to the water injection pipes, the upper water collector and the hose from the upper plenum to the vessel exit. These additional heat sinks must be represented because it has been observed  that the spurious condensation on inlet and outlet tubes is not negligible Fig. 4 shows the nodalisation of Vessel 1, as well as the subdivision of the heaters. The actual geometry of the vessels is approximately represented by blocking groups of cells, which are painted in dark grey. The nodalisation of Vessel 1 is similar to the one that was successfully used for previous analyses, which featured a rather coarse mesh (Andreani and Kelm, 2012). A few modifications were necessary to account for the size, number and position of the two heaters. Therefore, starting from a rather homogeneous mesh, with typical cell size of 0.2 m (the inner diameter of the injection pipe), the grid was "stretched" in the x-direction on the fully closed side of the heaters and "compressed" in the region between the heaters and the axis of the vessel. As for the vertical direction, the original mesh already included smaller meshes (10 cm) in the zone of the interconnecting pipe, due to the need to adequately resolve the counter-current flow between the two vessels. The global mesh was adapted for a smooth transition in the vertical direction between the exit of the lower heater and the fine mesh at the elevation of the interconnecting pipe. Finally, at the bottom of the vessel, a very large mesh (0.8 m in height) was used, considering that the modeling of this region would not have any influence on the simulation results. As regards the global quality of the mesh for Vessel 1, the cell size varies between 13 and 25 cm in the X-direction and between 20 and 26 cm in the Y-direction, respectively. In both cases, the maximum ratio between the sizes of adjacent cell was 1.5 and 1.2 in the X-direction and Y- direction, respectively. As for the subdivision in the Z-direction, if one disregards the lowest layer of cells, the minimum and maximum cell size were 10 cm, in correspondence to the lower heater exit, and 32 cm for the layer below the base of the condenser. The largest cell size ratio for adjacent cells was about 1.7. The subdivision of Vessel 2 was the same in the Z-direction to avoid that a poor quality of the mesh could negatively affect the gravity driven flows prevailing between the two vessels. In the X-direction, the cells increase gradually from 20 to 50 cm moving from the IP to the opposite wall, whereas in the Y-direction the cells increase from 20 to 30 cm moving away from the IP. The heaters appear in the nodalisation as blockages, since the heater housings are represented by separate volumes. The same applies to the condenser tubes. The obstruction plate is represented by blocking one cell face. The mesh of Volume 1 s consists of around 13'400 cells. Volumes 2 s and 3 s are represented with a coarser mesh, including 750 cells and around 7000 cells, respectively.
The mesh adjacent to the exit of the upper heater consists of one cell, because the details of the flow through and around this volume were shown to have a small effect on mixing (Andreani and Kelm, 2012). The exit zone of the lower heater instead is represented by two vertical subdivisions to permit a gross representation of the possible countercurrent flow at the exit, and its effect on the rising plume and its  interaction with the upper plume. The inlets to the heaters are represented in Vessel 1 by only one cell. For both heaters (cases), the crosssectional subdivision consists of 30 cells (6 × 5), whereas the vertical nodalisation is slightly finer for the lower heater than for the upper heater (8 instead of 7 vertical subdivisions, the additional subdivision being in the exit section, as shown in Fig. 4). The outlet area is subdivided in 4 × 3 cells for the upper heater, and 4 × 4 cells for the lower heater. Therefore, for the reference simulations, the 3-D connectors (3 and 5 for the inlets, and 4 and 6 for the outlets), which permit the flows between the vessel and the heater cases, hydraulically connect unmatched grids on the two sides. The cells in the lowest level (grey zone in Fig. 4) are partially blocked (57% porosity) to represent the reduction of flow area due to the heating elements (plates). A distributed thermal conductor with internal heat generation spans over all cells in this zone. In the reference calculations, the heat transfer from/to the housing was not considered, due to the difficulty to properly account for radiation heat transfer (Andreani and Kelm, 2012). The heat losses have only been represented in a sensitivity study carried out for tests HP2_2_2 and HP3-1-5, which addressed the effect of the exit gas temperature on the plume trajectory (Section 7).

Performance of the condenser
Test HP2_0_1 served to characterise the performance of the condenser and to observe its effect on the distribution of helium produced by the natural convection flow generated by condensation on the tubes. Fig. 6 shows the time history of pressure and power extracted by the condenser, showing the power of the tubes only and the total (which also accounts for the contribution of condensation on the inlet pipes and the outlet conduits). The power extracted is calculated from enthalpy balances on the water flow inside the tubes. It can be observed that the agreement for the pressure is excellent (the difference between calculated and experimental curve are barely visible) for the entire duration of the test. The power of the condenser is slightly overpredicted during the initial period of the steam injection phase (before helium injection), notably overpredicted during the helium injection phase, and again slightly overpredicted in the last phase of the experiment, the peak power being reasonably well captured. The overprediction in the early phase of the test is possibly due to the small disagreement of the steam vertical distribution (see below) at the very beginning of the test. It can also be remarked that the helium injection has a moderate effect on the pressure increase, but a larger effect on cooler power due to heat transfer degradation produced by the small increase of the non-condensable gases concentrations in the vessel. This degradation is underpredicted in the simulation, possibly due to a combination of effects (differences in the fine vertical distribution of helium, some reduction of accuracy of the condensation model in presence of a light gas, etc.). Similar to the experiment, the calculated powers extracted by the three tubes are approximately equal.
Sensitivity studies have been conducted on the effect of the cooler, and the contributions of the various heat transfer surfaces to the total power extracted. If the cooler were not activated, the pressure would have had a larger increase during the steam injection phase, and would have remained nearly constant after its termination. This sensitivity study shows that the effect of the condenser is already substantial after 1000 s, and becomes very large after termination of steam injection. Fig. 6 also shows that the power removed by the tubes, as well as the total power extracted by the cooler, which includes the spurious condensation on the inlet and outlet tubes and collector and accounts for up to 20% of the power extracted , Both powers are calculated fairly well. However the representation of the spurious condensation is not critical, because in absence of this additional heat sink the thermal load of the tubes would increase and the resulting pressure would be only slightly higher (see dashed lines in the figure) than for the case where all condensing surfaces are considered. 1 Another important aspect of the operation of the condenser is its capability to produce the transport of helium towards the bottom of the vessel and thus to the inlet of the lower heater. Indeed, the cooler was actually intended to reproduce the effect of prolonged condensation on the concrete walls of a plant on the transport of hydrogen to PARs located below the gas source. Fig. 7 shows the comparison between the calculated and experimental helium concentrations at various locations in the lower part of Vessel 1 for the first 4000 s. In the figure, also the results obtained for the hypothetical case without condenser are presented. It can be noted that the condenser has a substantial effect on the transport of helium to the lower part of the vessel, with the calculated timing of the concentration increase being close to that measured in the experiments. In the calculations (and to some extent in the experiments), the faster helium concentration increase in the lower header (Pos. T_20) than at higher position N_26 indicates that the accumulation in the lower part of the vessel partly occurs due to filling from below. This interpretation is supported by the observation of the flow field in Vessel 1 (Fig. 8, left), which shows that the condenser actually produces a downward oriented natural convection, with helium-rich mixture being displaced to the bottom of the vessel. The interpretation of the results is complicated by the counter-current flow in the IP, with return flow from Vessel 2 (Fig. 8, right), which produces at the elevation of the IP a strongly stratified layer, which "filters" the mixture composition in the downwards flow produced by the cooler. As already pointed out in the planning calculations , the effect of the condenser on the transport of helium would be somewhat stronger if no connection existed between the two vessels. It has to be considered, however, that due to the steam stratification resulting from the injection at 4 m, the condensation rate (and associated natural convection) below the injection level is much smaller than in the region above (about a factor 2 according to calculations), and therefore the return flow from the second vessel has a minor impact on the evolution of the thermal-hydraulic conditions in Vessel 1. Fig. 7 (right) presents the results for the helium concentration at the lower heater inlet for the test without heaters (HP2_0_1) and for the two tests (HP2_1_2 and HP3_1_5) where the upper heater was activated. It is shown that the energy input from the upper heater has practically no effect on the helium transport, as the density decrease due to the higher temperatures of the thermal plume is negligible with respect to the large density differences due to concentration gradients produced by the injections and the condensation on the cooler tube surfaces.

Mixing produced by one heater in tests HP2_1_2 and HP2_2_2
The upper heater is activated in Test HP2_1_2 (as well as in test HP3_1_5, which is discussed in the next section) immediately after the start of the helium injection, and produces a thermal plume that interacts with a rather homogeneous ambient above the IP, resulting from the injections and the operation of the condenser. Fig. 9 shows the 1 The "adaptivity" of the power that can be extracted by the cooler to modifications of the total area of the condensing surfaces results from the nonlinearity of the driving force controlling heat/mass transfer, which depends on local values of single-phase water temperature, tube wall temperature and steam partial pressure in the vessel. This complex dependence of global heat transfer on distributions of variables on both sides of the tubes results in total condensation rate not linearly proportional to the total heat transfer area. development of the steam and helium vertical distributions (along the vessel axis) in Vessel 1 at various times, comparing the results with those obtained in test HP2_0_1 (without heater). In general, the agreement between calculations and experimental results is good, with a small deviation in the steam distribution at the very beginning of the steam injection and the only large deviation in the helium distribution occurring during the helium injection (for instance at 700 s in Fig. 9, right), where the calculations seem to indicate a rise of concentration with elevation at the highest positions. This is due to the bending of the plume above the obstruction plate (not shown), which originates from the interaction of the jet with the obstacle and with the flow produced by the condenser. This behavior was not observed in previous analyses of tests with obstructions (Andreani et al., 2020), but in the cases analysed earlier the density difference between injections and ambient was much smaller and no condenser was installed.
As for the effect of the operation of the heater, the thermal plume produces a little faster mixing of the fluid, decreasing to the elevation of the heater inlet (~5 m) the height of stratification front. This behavior   confirms the previous findings for tests in PANDA (Andreani and Kelm, 2012). Another beneficial effect is to reduce the persistence of the peak helium concentration reached at the end of the injection and prevent the long-term rise of helium concentration produced by condensation (as shown by the lower helium concentrations above 6 m at 9000 s in test HP2_1_2).
In test HP2_2_2, the lower heater was activated 1000 s after termination of the steam injection. This implies that the thermal plume produced by the heater had to rise in an ambient that at the time of the heater activation was strongly stratified. The evolution of the flow produced by the plume is shown in Fig. 10. Since density distribution depends more on mixture composition than on temperature for the range of conditions investigated in these experiments , the thermal plume is initially negatively buoyant with respect to the upper region of the vessel, and remains confined in a region below the inlet of the upper heater for a long time after the heater has reached at 4000 s its maximum power (Fig. 10, left). Later, after the stratification erosion induced by the plume has gradually reduced the strength of the stratification, the plume can finally rise to the top of the vessel (Fig. 10, right), and rapidly mixes the fluid in the entire vessel. Fig. 11 shows the evolutions of the vertical distributions (along the vessel axis) of steam and helium concentration (along the axis above the pipe exit and at 1.43 m from the axis below that level) and the gas temperatures (at 0.65 m from the axis, between the axis and the heater). The calculated results are generally in good agreement with data, except for the timing of the erosion of the stratification under the effect of the thermal plume produced by the lower heater. Indeed, the temperatures at 4000 s are close to the experimental values, which confirms that the plume could not rise during the first period after the activation of the heater. Later, the calculation shows a too fast progression of the leading edge of the plume, which, as shown in Fig. 10, at 5700 s reaches the top of the vessel. In the experiment, however, the temperatures at the highest elevations rise only at 6200 s. The time histories of steam and helium concentrations also show that the calculation overpredicts by 500 to 700 s the rate of upwards displacement of the density interface (stratification front). These results seem to contradict previous results (Andreani and Kelm, 2012), where the calculations showed a very good agreement for the time of the mixing, starting from similar conditions. The reason for the overprediction of the mixing rate is the slight underprediction of the inclination of the thermal plume exiting the heater, which causes the flow to avoid nearly completely the interaction with the obstruction plate. In the experiment (Fig. 12), in the PIV Fieldof-View, the velocity field appears to be strongly affected by the presence of the plate, which slows down the upwards propagation of the flow structure. In the calculation, however, only the edge of the plume interacts with the obstruction, and therefore the flow can rise nearly unhindered in the vessel. Since the experimental information on velocities at the inlet and outlet of the heater is not available, it is not possible to provide a complete explanation for the different inclination. Some considerations, however, can be made.
First of all, it can be inferred from the comparison of the temperatures at three positions close to the heater exit (Fig. 13, left) that the average temperature is overpredicted, since the value at the lowest position is strongly overpredicted. Therefore, it is likely that the calculated gas outflow is hotter than that in the test, possibly due to the neglect of heat losses from the housing. On the other hand, the measured lower temperature at the bottom of the exit section indicates the possibility of reverse flow, which probably produces high exit velocity at the top of the exit section, and therefore a more horizontal flow. The details of the flow at the exit section cannot be captured with a coarse Fig. 9. Gases vertical distributions for tests HP2_0_1 and HP2_1_ 2: steam (left) and helium (right).    mesh. Another reason for the discrepancy could be the too small temperature drop outside the heater. In fact, the calculated temperatures at level L (0.6 m above the heater) are higher than the measured ones (Fig. 13, right). This can be due to underprediction of mixing, but also to the neglect of radiation heat transfer, which in reality tends to reduce the gas temperatures. Both effects contribute in the calculation to the too strong buoyancy of the flow in vicinity to the heater exit and thus cause the nearly horizontal jet exiting the housing to evolve into a nearly vertical plume at too short distance from the heater.
A sensitivity study with a model including heat losses and a finer mesh in the heater exit region has then be performed to investigate the effect of the plume trajectory on the global mixing. This study is presented in the next section, where the test with two heaters is discussed in more detail.

Interaction of the two heaters
In test HP3_1_5 both heaters were activated, the upper one at 500 s, the lower at 3000 s. Fig. 14 presents the results of the velocity and temperature fields at various times. The rise of the thermal plume produced by the lower heater is hindered by the stratification at the start of its operation, similarly to test HP2_2_2, and the operation of the upper heater affects only marginally the rate at which the stratification is eroded. In principle, this result was not obvious, because the heat release from the upper heater would tend to decrease the density above the heater inlet, and thus increase the resistance against the mixing produced by the lower plume. On the other hand, the mixing produced by the upper heater reduced the stratification within the upper region of the vessel, which reduced this resistance. Due to the overwhelming effect of the mixture composition on density stratification, the changes due to the moderate temperature increase was practically negligible, and therefore the mixing in test HP3_1_5 was close to that in the test without upper heater. The comparison of the calculated evolution of the gas temperature vertical distribution on the axis (Fig. 14, right) above the pipe injection exit with measured data shows that simulation correctly represents the balance between the various effects. However, as shown by the general trends in the distributions at 5700 s and 6200 s the lower plume in the simulation merges with the upper plume some hundreds seconds earlier than in the test: in fact, at 5700 s, the experimental profile (blue triangles in Fig. 14, right) shows a decreasing trend between the lower heater exit and about 6.3 m, and an increasing trend above the upper heater outlet. This shows that at 5700 s the lower plume was still confined below the upper hot layer produced by the upper heater, or, in other words, the two plumes had not yet merged. In the simulation, however, already at 5700 s the profile is monotonic, which shows that the two plumes had already merged. The faster rise of the lower plume, similarly to test HP2_2_2, causes the calculated mixing process to be faster than in the experiment (this will be more evident in the time history of the helium concentrations shown below). For test HP3_1_5, no PIV data are available, but it can be argued that, similarly to tests HP2_2_2, the calculated plume is less inclined than in the test, and this results in less interaction with the obstruction plate, and thus in faster upwards propagation of the plume. Indeed, the comparison of the horizontal gas temperature distribution at 0.6 m above the heater exit (not shown) shows a smaller lateral distance of the calculated peak value from the heater than in the test (similarly to Fig. 13, right). In order to verify this conjecture, a simulation with a slightly modified model for Vessel 1 was performed for both tests HP2_2_2 and HP3_1_5. In this refined nodalisation (Fig. 15), the density of the cells in the regions of the two heaters was increased to match the grid inside the volumes representing the heaters, and heat losses from the heater housings were considered. The subdivisions of the two heater cases were not modified. The total number of cells for Vessel 1 in the refined mesh thus increased to around 25 ′ 600 cells. The computational time for simulations with this mesh increased to about 5 days. It has to be considered that the refined mesh is far from being of "CFD quality", because, especially in the Ydirection, the very large cell size change near the heaters is expected to lead to deterioration in the accuracy of the results. However, due to the very large computational times (several weeks) that a more regular grid would have implied, a denser mesh respecting best practice rules could not be afforded. Since the interest here was to explore whether we could get the indication that a different mesh could produce a more accurate flow pattern outside the lower heater, and this, in turn, could result in a more accurate calculation of the time of complete mixing, the results are presented, although they are not completely trustworthy. Fig. 16 (left) shows that the resulting thermal plume was more inclined, and strongly interacting with the obstruction plate. Separate sensitivity studies revealed that this change in the results was mostly due to the refined nodalisation, and only marginally affected by the consideration of the heat losses. The modification of the flow structure is also illustrated by the horizontal distribution of gas temperatures (Fig. 16, centre), where the position of the peak value is practically coincident with that inferred from the experimental measurements. It can be observed that the calculated temperatures, at some distance from the heater outlet, are overpredicted, although the temperatures at the heater outlet (Fig. 16, right) are underpredicted (except for the lowest position BG27). The underprediction of the heater exit temperatures is possibly due to the overprediction of the heat losses (but it could also due to effect of the Fig. 14. Test HP3_1_5: gas temperature and velocity fields at three times (left) and evolution of temperature vertical distribution above pipe exit elevation (right). direct radiation from the housing and the hot gas on the reading of the thermocouples). The too small temperature drop between the heater exit and the level L could be due to one or both of two effects: 1) underprediction of mixing of the bent flow; 2) neglect of radiation heat transfer from the hot plume to the surroundings over the large length of nearly vertical flow.
The latter effect is likely to be more important because for the upper heater, where the available horizontal profile is at an elevation closer to the top of the heater, the temperatures are still underpredicted (not shown). This observation would suggest that the disagreement does not develop in the immediate vicinity of the heater, where the flow is bent, but at some distance, where the flow is approximately straight and mixing is usually overpredicted using a coarse mesh, rather than being underpredicted. In this far-field zone, the neglect of radiation heat transfer is probably causing the overprediction of temperatures.
Finally, Fig. 17 shows the time histories of helium concentrations at various positions, calculated with the reference and the refined mesh, and, for the sake of comparison, the results for test HP2_2_2. It is noted that the agreement using the reference mesh is slightly better for test HP3_1_5, where the mixing was a little faster, because the upper heater produced a less strong stratification at the time of the start of operation of the lower heater. This effect is captured by the code using the reference, coarse mesh. Better agreement for both tests is obtained with the refined mesh, since the lower thermal plume is more inclined and is hindered in its upwards propagation by the interaction with the obstruction plate (Fig. 16). In all simulations, the largest disagreement between calculated and experimental results was observed for position L_26 during the steam injection. Since this location is at height of the top  of the interconnecting pipe, the mixture composition is strongly affected by the complex flow in the pipe and thus in the middle of strong gradients. Since the focus of this paper is on the effect of heaters, and this is not influenced by the precise position of the interface within the interconnecting pipe height, this specific aspect of the simulations has not been discussed.
Since the main interest in this work (as well as in previous works over the last two decades at PSI) was to test a coarse mesh approach for a transient controlled by stratification and mixing, no comprehensive mesh sensitivity study has been performed. On the one hand, a much coarser mesh of Vessel 1 would not be adequate, because it has been observed in all past studies that a mesh with cells much larger than the size of the injection pipes produces excessive numerical diffusion during the stratification build-up. Moreoever, due to the need to resolve the counter-current flow between the two vessels, the subdivision of the Vessel 1 in the Z-direction close to the interconnecting pipe must be sufficiently fine. A much coarser subdivision for the remaining part of the vessel would thus generate a strongly irregular mesh with unknown consequences for the trustworthiness of the results. On the other hand, a homogenously finer mesh (e.g. halving the size of the mesh in all direction everywhere in the fluid domain) than that used in the this work would result in very large computational effort, which could not be afforded and, in authors' opinion, would be in contrast with a proficient use of GOTHIC, which is best suited for reasonably accurate, but fast running simulations. The goal of mesh independent results is practically not achievable with GOTHIC for the present complex situation, and therefore a comprehensive mesh sensitivity study is not performed. The last section is instead dedicated to a partial sensitivity study, which addresses the modelling of the heaters, with the intention to derive some suggestions for plant analyses.

Sensitivity studies
In this section, a series of parametric studies for Test HP3_1_ 5, with both heaters in operation, will be presented, which address the sensitivity of the main results to the details of the modelling of the heaters. In addition to the two meshes considered in the previous section (the "Reference" model using unmatched grid at the boundary between heaters and vessel, and the "refined" mesh using matched grids), the use of two other meshes with matched grids for the volumes representing the heaters will be investigated. These feature a single-cell or 2x2 cells representation of the cross-sectional areas of the heaters, respectively. Since the 2x2 subdivision was adopted for the SETH project (Andreani and Kelm, 2012), this model will be referred to as the "SETH" mesh. The sensitivity study is completed by the "basic" model, with only one cell representing the cross-sectional area. The three models ("basic", "SETH", and "refined") with matched grids will also be used in conjunction with a different approach, where, instead of using separate volumes for the heaters, connected to the vessel by 3-D connectors, the heaters boxes are embedded in the mesh of the vessel, with the box walls represented by partly blocking some cells (Andreani and Kelm, 2012). Such an approach is illustrated in Fig. 18 for the case ("basic" mesh) with the coarsest subdivision of the heater boxes (only on cell representing the inlet). It can be recognized that the fluid inside the box is surrounded by blockages of small thickness, which prevent the lateral and top flows between heater boxes and vessel, and force the flow through the inlet and the outlet of the heaters. Again, the light grey area represents the cell where the heater plates (heat sources) are located.
The sensitivity study thus includes the seven simulations listed in Table 2 (two already discussed in the previous sections). In principle, the two approaches (with the same subdivisions) should lead to the same results, because the numerical method for solving the conservation equations on the 3-D connectors is fully consistent (EPRI, 2016) with the one-domain mesh set-up. However, probably also due to small differences in the results of the pre-processing for hydraulic diameter of the cells and the slightly different intrinsic pressure losses in the two approaches, slight differences may appear, which, in long transients, grow with time. Indeed, tests with the Basic and SETH meshes applied to a single heater showed, for nearly steady-state conditions, that the solutions for both temperature and velocity fields were not fully identical for the separate volume approach and the embedded mesh representation of the heater. For both meshes, the velocity at the exit of the housing is slightly smaller in the separate volume approach. Obviously, neither of the two solutions can be considered accurate, because the complex flow field inside and around the heater cases cannot be simulated on a very coarse mesh.
The results of the sensitivity studies are shown in Fig. 19. In the first (Fig. 19 a)), the results obtained with the reference mesh are compared with the Basic and E-Basic mesh, because all these models feature a single node representation of the heater inlet and two-cell representation of the heater outlet on the vessel side. It is observed that the basic mesh performs as well as the Reference mesh, which indicates a finer subdivision of the volume representing the heater is actually not required for obtaining globally successful results. The figure also show the somewhat surprising result that the embedded mesh permits to reproduce perfectly the mixing at the highest elevation considered. The difference in the results using the Basic mesh and those obtained with the E-Basic mesh can only be explained with small differences in the numerical solution, which grow with time (see above). Fig. 19 b) compares the results of the two approaches for the SETH mesh. Again, the results using the embedded mesh are slightly better than those using separate volumes for the heaters. Finally, Fig. 19 c) shows the predictions with the Refined mesh: in this case, the simulations predict practically identical results. Fig. 20 shows the velocity field predicted with the three embedded meshes: in all cases, the plume produced by the lower heater fully impacts the obstruction plate, which confirms that the key for the accurate prediction of the mixing is a correct representation of the flow interaction with that obstacle.
In summary, the sensitivity study shows that, independently of the model chosen for representing the heaters, the global mixing process can be successfully predicted using a coarse mesh. Even a very coarse modelling of the housings with a few cells (one for the inlet, two vertical subdivision for the outlet, a vertical stack of a few cells to represent the chimney effect of the PAR case) seems to be adequate for simulating the mixing produced by the hot plume generated by the heaters. As for the choice of the specific approach (separate volumes or embedded mesh), the present study only shows that for the specific conditions simulated, the models with embedded mesh resulted in less inclined heaters exit plumes, and nearly mesh-independent and more accurate results. This could not be the case under other conditions, since the two approaches are in principle equivalent. It can be concluded, however, that the use of both approaches can be recommended for parametric studies, especially in cases where the trajectories of the plumes could affect the accuracy of the predictions (like in the case of a possible interaction with a structure located at short distance from the heater exit).
The relevance of these analyses for full-scale containment applications are discussed in the conclusions.

Conclusions
The experiments performed within the HYMERES project addressing the effect of heat sources on mixing in presence of condensation and a flow obstruction, as well as the interaction of two sources at different elevations, have been analysed using the GOTHIC code and a coarse mesh. In general, good agreement has been obtained for all cases. In particular, the effect of condensation on the downwards transport of light gas has been properly represented, although this occurred slightly faster in the calculation than in the tests. In addition, the small effect of the operation of the upper heater on the distribution of gases was well reproduced. In this case, the mixing effect of the heater, in agreement with previous findings, is limited to the upper part of the vessel above the inlet of the heater housing. The global mixing produced by the operation of the lower heater alone could also be calculated fairly well, and the small acceleration in the mixing at the highest elevations resulting from the simultaneous operation of the two stacked heaters could be captured.
The only notable discrepancy that has been observed was the faster  upwards propagation of the thermal plume generated by the lower heater in both tests where this was activated. The comparison of the velocity field with PIV measurements in one case and horizontal profiles of gas temperatures for both tests revealed the unexpected influence of the flow obstruction located above the injection. In fact, the interaction of the flow issuing from the lower heater exit was not correctly represented with the reference model, because the thermal plume was too buoyant and barely touching the plate. A refined model using a separate control volume, and all meshes where the heater boxes were integrated into the global computational grid permitted to reproduce the trajectory of the plume, and this resulted in better results. These results confirm previous analyses, i.e. that a coarse mesh is capable to represent the main phenomena resulting from the operation of heat sources, and the new analyses show that this is true also in presence of distributed condensation. However, for cases where the flow produced by the heat sources could interact with structures (in this case was a simple horizontal plate), the mixing is affected by the orientation of the jets/ plumes, and the effect of the flow pattern should be investigated by means of parametric studies using different meshes and/or resorting to CFD codes. The analyses presented in this work are expected to be of some relevance for plant analyses, although the typical modelling approach with GOTHIC is to use the PAR component (EPRI, 2016;Papini et al., 2018), where the metallic housing is not explicitly represented, and the chimney effect produced by the tall case is modelled by an ad-hoc treatment of the gravity head and use of adequate design parameters. In the case, however, that the housing is explicitly modelled (for some of the components in critical areas where a combustible mixture can buildup and/or for sensitivity studies), the present results provide some useful indications. For plant analyses, the present studies thus suggest that, depending on user convenience (in authors' opinion the use of separate volumes simplifies the pre-processing), either of the two approaches (using separate volumes or embedding the housings in the global mesh) can be chosen, without concern about a possible degradation of the global results. In both cases, the subdivision of the heaters in a few cells is sufficient to obtain reasonably accurate results. Moreover, if close to the exit of some PARs there are structures with which the hot plumes can interact, parametric studies using both approaches should be considered, because they could lead to slightly different flow patterns and different mixing times. If such discrepancies appear, and these are reckoned to be of some importance, the user is advised to resort to more detailed modelling, e.g. using CFD codes.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.