Multiscale modelling of the phénix dissymmetric test benchmark

Nuclear Engineering and Design xxx (xxxx) xxxx

Contents lists available at ScienceDirect

Nuclear Engineering and Design journal homepage: www.elsevier.com/locate/nucengdes

Multiscale modelling of the phénix dissymmetric test benchmark ⁎

H.J. Uitslag-Doolaarda, , F. Alcaroa, F. Roelofsa, X. Wangb, A. Krausc, A. Brunettc, J. Thomasc, C. Geffrayd, A. Gerschenfeldd a

Nuclear Research and Consultancy Group (NRG), Petten, the Netherlands Karlsruhe Institute of Technology (KIT), Karlsruhe, Germany Argonne National Laboratory (ANL), Argonne, USA d DEN-Service de Thermohydraulique et de Mécanique des Fluides (STMF), CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, France b c

A B S T R A C T

Validated thermal-hydraulic calculation methods are essential for safety assessments of nuclear reactors, since experiments are at least expensive and mostly impossible or undesirable. In the framework of the validation of thermal-hydraulic codes for liquid metal-cooled pool-type reactor applications, the end-of-life dissymmetric test performed in the Phénix sodium reactor was selected to carry out a benchmark activity. This paper covers the achievements of the multiscale models on the simulation of this benchmark. The Phénix sodium-cooled fast reactor contains two pools in which characteristic 3D effects like thermal stratification and jets play a significant role. One of the two active heat exchanger sets fails in the dissymmetric test, leading to changes in the thermal stratification and asymmetric conditions in the reactor. Therefore 3D modelling using Computational Fluid Dynamics codes (CFD) is required to model this transient in sufficient detail. However, several aspects of the reactor, like the core and the secondary side of the heat exchanger, can be modelled more efficiently by System Thermal-Hydraulics codes (STH). Therefore multiscale models are developed by coupling CFD and STH codes, exploiting the strengths of both codes in terms of modelling options and computational costs. The results provide insight in the local flow patterns and the temperature distribution during the dissymmetric transient. The results and modelling choices of the participants are compared in this paper and the main findings are discussed in the light of lessons learned from the blind phase of the benchmark and providing recommendations for future modelling of pool-type liquid metal-cooled reactors.

The Phénix reactor is an experimental liquid sodium-cooled fast reactor operated by CEA in France from 1973 till 2009. Several end-oflife transients have then been performed for producing measurement data for better understanding of accident scenarios in sodium-cooled fast reactors. Within the H2020 SESAME project, an effort is made on further development of the thermal-hydraulic safety analyses for metal-cooled pool-type reactors. Usually such safety analyses are performed with 1D system thermal-hydraulics codes (STH). STH codes were developed and validated in the past decades for computing the thermal-hydraulics in complete reactor systems for nominal and accident scenarios. In the past decade, the 3D computational fluid dynamics codes (CFD) have been more and more applied for reactor safety analyses, since these provide detailed insight in local flow patterns and thermal stratification. However, CFD calculations are much more computationally expensive than STH calculations. Therefore, multiscale simulations, combining the strong parts of both types of codes, have been further developed in the SESAME project, where different STH and CFD codes are coupled by multiple parties. Several benchmarks based on measurement data in experimental facilities or from actual reactors have been defined in the SESAME project for validation of these multiscale ⁎

models. One of these benchmarks is dedicated to the end-of-life dissymmetric test performed in the Phénix reactor. The dissymmetric test starts from conditions similar to nominal conditions, at a power level of about 340 MWth, with 3 primary pumps (PPs) and 4 intermediate heat exchangers (IHXs) fully operating. The IHXs that belonged to the second PP were replaced by so-called DOTEs (fake IHXs), which are inactive components that do not let the sodium between the hot and cold pool. This transient starts with the trip of one of the pumps in a secondary loop belonging to two IHXs. This causes an asymmetric temperature distribution in the cold pool due to different outlet temperatures at the primary side of the IHXs of the different loops. Subsequently, the reactor automatically shuts down, amongst others leading to a pump speed reduction in the secondary loop of the remaining two IHXs and a reduction of the core power to decay heat. After the scram is completed, the measurements were continued till the total transient time of 1800 s. One of the main challenges in this benchmark is to model the temperature distribution, having stratification and asymmetry, in the pools and the IHXs. To properly model this, accurately resolving the flow patterns determined by a balance of forced flow on one hand and

Corresponding author. E-mail address: [email protected] (H.J. Uitslag-Doolaard).

https://doi.org/10.1016/j.nucengdes.2019.110375 Received 28 May 2019; Received in revised form 23 September 2019; Accepted 3 October 2019 0029-5493/ © 2019 Elsevier B.V. All rights reserved.

Please cite this article as: H.J. Uitslag-Doolaard, et al., Nuclear Engineering and Design, https://doi.org/10.1016/j.nucengdes.2019.110375

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

buoyancy forces on the other hand is of critical importance. This balance is the most delicate in the cold pool near the IHX outlet windows and the PP skirt inlets. This paper proceeds with a summary of the main components of the Phénix reactor, followed by a description of the multiscale models developed by each benchmark participant. This includes highlights of improvements based on the lessons learned from the blind phase. Subsequently, the results that were obtained in the open phase are presented, compared among the participants and compared to the measurement data. The discussion of these results leads to the conclusions and recommendations for further development of the Phénix model and the multiscale modelling of pool-type liquid metal-cooled reactors in general. 1. The phénix reactor The Phénix reactor was a liquid sodium cooled fast reactor operated by CEA in France from 1973 till 2009. Over this time, it was used both for experimental purposes and for power generation. It helped acquiring knowledge about the sodium cooled reactor technology in numerous fields (neutronics, materials and thermal-hydraulics). It also allowed for the validation of technological choices for this kind of reactor designs and demonstrated that good breeding performances could be achieved, thus proving that fast reactors could be used to better use the natural resources and reduce the radioactive waste. More details about the Phénix reactor and the experience acquired during its operation time can be found in (Guidez, 2013). The primary sodium of the Phénix reactor, schematically shown in Fig. 1 (from IAEA, 2013) in side-view, is contained in the hot pool (red), the cold pool (green) and the vessel cooling bypass system (green outer part). The internals of the reactor consist amongst others of three PPs

Fig. 2. Schematic overview of the Phénix reactor as compiled top-view.

(see Fig. 2). Two of them are surrounded by a couple of IHXs and the third one (PP2) is surrounded by two DOTEs, which are hollow cylinders replacing the space of two IHXs. Following the flow path (see Fig. 3 from IAEA, 2013) of the primary sodium starting from the pumps, the sodium enters the diagrid, which is located below the core, where about

Fig. 1. Schematic overview of the Phénix reactor as compiled side-view. 2

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Fig. 3. Sodium flow in the primary circuit.

Fig. 4. Top view of the core. 3

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Fig. 5. top-view of the CFD model focussed on the core, showing the thermo-couple tree ‘COLTEMP TC pole’ as red dots.

10% of the sodium flows downwards into the vessel cooling bypass system and 90% upwards through the core. The sodium enters the hot pool through the fuel assembly (FA) outlets, partly impinging at the core plug and the fuel handling system. The core (see Fig. 4) is surrounded by so-called “logs”, being neutron shielding steel and graphite cylinders submerged in the hot pool. The sodium in the hot pool proceeds to the cold pool via the IHXs, subsequently re-entering the PPs. The sodium that had entered the vessel cooling bypass system enters the cold pool after cooling the vessel. The Phénix reactor was equipped with different temperature measurement devices:

2. Multiscale model NRG The multiscale model developed at NRG (Petten, the Netherlands) (Uitslag-Doolaard, 2018) is a code coupling between the STH code SPECTRA and the CFD code ANSYS CFX. SPECTRA is a thermal-hydraulic system code applicable for transient analyses in Light Water Reactors (LWR), High Temperature Reactors (HTR) and Liquid Metal Fast Reactors (LMFR) (Stempniewicz, 2009). It is selected because of experience in LMFR modelling and because NRG has access to the source code. ANSYS CFX (ANSYS CFX, 2013) is a CFD software tool that delivers reliable and accurate solutions quickly and robustly across a wide range of CFD and multi-physics applications. CFX is designed for detailed three-dimensional analyses of fluid flow and (conjugate) heat transfer in the fluid and solid structures in both steady state and transient.

- At the PP inlet (measured at the top of the pump skirt); - Several thermo-couples in the hot pool above the shielding of the core (red dots showing the thermo-couple tree called ‘TC pole COLTEMP’ in Fig. 5); - At the IHXs inlet on the primary side; - Thermo-couple tree in the cold pool next to IHX M, also covering the temperature at the outlet window (line of brown dots called ‘IHX M TC pole’ in Fig. 6); - Thermo-couple tree in the cold pool below a DOTE (line of blue-grey dots called ‘TC pole DOTE’ in Fig. 7); - At the inlet and outlet of the secondary side.

2.1. Spectra model In the SPECTRA model, both the cold pool and the hot pool are divided in three azimuthal sections of each 120°. A section of the SPECTRA model of the pool is provided in Fig. 8 (left). The pools are confined by the vessel cooling bypass region. The external vessel that surrounds the bypass cooling region is included in the open phase model and is designed to a heat loss of ~520 kW to the environment under nominal conditions. The heat exchange between the pools and bypass through the vessel walls is also included in the open phase, though computed by CFX and provided to SPECTRA. Heat structures are present in different portions of the pool to account for the reactor internals. Although the geometry details are not provided, the structure masses are accounted for in order to preserve the heat capacities. For

Additionally, the core power and the primary and secondary pumps rotational speeds were measured. The coming sections describe the multiscale models of the primary and secondary system of the Phénix reactor developed by the benchmark participants.

4

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Fig. 6. side-view of the CFD mdoel showing the thermo-couple tree ‘IHX M TC pole’ as brown dots.

2.2. Cfx model

the open phase of the benchmark, the so-called LOGS elements are implemented assuming a steel mass of 130 t and 56 t of graphite, since it followed from the analyses of the blind phase that its heat capacity has a significant influence on the IHX inlet temperature. These elements surround the core and receive no flow but are rather standing in contact with the hot pool sodium. Heat loss from the hot pool through the roof is taken into account using a specific heat loss correlation provided as boundary condition for the benchmark. The core model consists of eight types of sub-assemblies (SAs): inner core with free or linked fuel, outer core with free or linked fuel, fertile blanket, reflector/shielding, control elements and spent fuel elements. All SA types are uniformly subdivided in the three azimuthal sectors and connected to the corresponding diagrid sectors, at the inlet, and hot pool sectors, at the outlet. The nodalization scheme adopted for the fuel SAs is provided on the left side of Fig. 9. The SPECTRA model of the IHXs is provided at the right side of Fig. 9. The primary sodium flows downwards in the shell side exchanging heat to the secondary sodium which flows upwards in the tubes, coming from the central downcomer. The downcomer and the IHX outer shell walls are insulated, therefore the heat exchange occurs only through the HX tube walls. The red parts in Fig. 9 represent the locations where the boundary conditions are imposed for the secondary side. The blind phase contained two sets of duplicated IHXs, which is expanded to four separate IHXs in order to better catch the potential asymmetry. The PP characteristic is derived from the homologous curve data provided in the benchmark input data file. The sodium properties and the corresponding heat transfer correlation are implemented in SPECTRA from Fink et al. (1995), Lyon (1951), Churchill (1977) and Mikityuk (2009).

The geometry for the CFD model is provided by CEA and is split into the three azimuthal sections. The CFD model contains the two pools; the IHXs, the DOTEs, the core plug, the fuel handling system, the vessel cooling bypass system and the skirts of the PPs. The pools were hydraulically separated by a cut in the IHXs, but the blind phase results showed that complete IHXs in the CFD code are essential to properly compute the flow through and the stratification in the IHXs. Conjugate heat transfer is included in all vessels, in the PP skirts and in the caps at the top of each IHX. The core outlet and the vessel cooling bypass system are the inlets of the CFD model and the PP inlets form the outlets. The CFX mesh consists of structured hexahedral cells where possible, e.g. the PP skirts, the IHXs and the major part of the hot and cold pools, supplemented with tetrahedral cells, comprising of 1.17 million cells in total. The core outlet is modelled as five concentric rings, getting the weighted average temperature and cumulative mass flow rate of the 1) inner core and control SAs; 2) outer core SAs; 3) blanket SAs; 4) spent fuel SAs and 5) shielding SAs. The LOGs are taken into account as surface surrounding the core, acting as a heat source or sink in the hot pool according to the input from SPECTRA. The IHXs are divided in the same axial sections as in the SPECTRA model, such that the heat exchanged to the secondary side in SPECTRA is applied as volumetric heat sink in CFX. Summarizing the computational setup, the SST k-ω turbulence model is applied with buoyancy turbulence production, the same temperature-dependent properties for sodium are applied as in SPECTRA, a constant time-step of 0.05 s is applied, the turbulent Prandtl number is set to 2 and the pool levels are modelled as fixed free slip walls. 5

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Fig. 7. side-view of the CFD model showing the thermo-couple tree ‘DOTE TC pole’ in the cold pool below the DOTE as grey-blue dots.

these results are provided as boundary conditions at the coupling interfaces and the initial conditions in CFX stand-alone are computed. Eventually, SPECTRA and CFX converged to a coupled steady state, which constitutes the starting point for the transient.

2.3. Coupling setup To obtain the initial conditions for the dissymmetric transient, first a stand-alone SPECTRA steady state result is computed. Subsequently

Fig. 8. SPECTRA model nodalization developed by NRG of one azimuthal section of the pool (left); CFD model developed by NRG (right). 6

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Fig. 9. Nodalization of a fuel assembly (left) and an IHX (right) in NRG’s SPECTRA model.

pumped back to the diagrid. About 10% of the coolant from the diagrid flows into the vessel cooling loop, and enters the cold pool from the surface region. The hot pool, cold pool and vessel cooling bypass system are connected by the argon layer with the cross connections. The argon gas works like a pressurizer to stabilize the system pressure. The most important components in the primary loop are the 3 primary pumps, connecting the cold pool and the diagrid. Each pump is discretized into several parts: the pump inlet skirt is simulated as pipe component; the turning point is modelled as a branch which is connected to the pump model; the pump outlet is also modelled as a pipe which is connected to the diagrid. In ATHLET, the differential pressure control model is used for the pump component, and the pump speed is given as a constant value of 540 rpm. The core is divided into 6 zones, including inner core, outer core, blanket zone, steel shielding zone, control rod and safe rod zone, and incontainment storage fuels. In ATHLET model, the flow path of fuel assembly in each zone is modelled. In each core zone, one fuel assembly is modelled and multiplied by a factor of the amount of fuel assemblies in this zone. The cross section of the flow path inside fuel assembly varies in the flow direction. In the axial direction, the fuel assembly is divided into bottom nozzle, expansion plenum, blanket, fission fuel, expansion plenum, upper axial blanket, and shielding. The sodium flow path is thermal connected with the fuel rod model. The fuel rod is modelled as several layers, the inner part of rod is the MOX fuel pallet, and the outer part is the cladding material austenitic steel. The internal heat source is given as a user input.

An explicit domain overlapping coupling method is used. The full system is modelled in SPECTRA as described above. The SPECTRA model contains an accurate representation of the heat transfer phenomena in the core and the IHXs. Furthermore, the flow distribution in the PPs is well computed at system level, though the temperature of the incoming sodium is predicted more accurately by CFX. The results in the pools, IHXs and the temperature distribution in the skirts of the PPs are computed in CFX and overwrite the SPECTRA results in order to resolve the thermal stratification and complicated three-dimensional flow patterns. The mass flow rate, temperature, heat flux and heat source or sink are exchanged between CFX and SPECTRA after each time-step via external data exchange routines. The data-exchange timeinterval equals the CFX time-step, resulting in a RMS courant number of 0.22 during the complete transient simulation, since the PPs are operating at constant speed along the whole transient. 3. Multiscale model kit 3.1. Athlet model As shown in Fig. 10, the pool type reactor is filled with working fluid of sodium, and the top of the system is covered by the argon gas. The sodium coolant from the diagrid goes through the reactor core, and is heated by the fuel rods. The thermal power in 6 different core zones is given as an input. After heating by the core, the coolant flows into the hot pool, which is modelled as several big branches and pipes. Then the coolant from hot pool is split into 4 intermediate heat exchangers, and the heat is released through thermal conduction to the secondary side. In the secondary side, the time dependent inlet and outlet boundary are given as input for two secondary loops. Each secondary loop goes into two heat exchangers. After cooling by the 4 heat exchangers, the sodium coolant enters the cold pool, then mixes with the coolant in the bottom region, and flows into the pump inlet skirt, and finally it is

3.2. OpenFOAM model The hot pool and cold pool are isolated in the CAD model of Phénix. The core, heat exchangers and pumps are simulated by the STH code. The CAD models provided by CEA for the separated hot pool and cold pool are used to create the CFD models. As shown in Fig. 11, the inlet 7

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Fig. 10. ATHLET model for Phénix.

results are used as the constant boundary conditions for the CFD steady state simulation. The core outlet mass flow rate and temperature are used as the inlet boundary of the hot pool, while the IHX inlet pressure is used as the outlet boundary of the hot pool. For the cold pool, the IHX outlet mass flow rate and temperature are given as the cold pool inlet boundary, while the primary pump inlet pressure is used as the outlet boundary of the cold pool. After finishing the steady state calculation of CFD and STH, the results of the STH simulation are taken as the initial field of the transient simulation. The CFD part of the coupled simulation is initialized with the flow field from this steady state. However, the initial temperature is assumed to be uniform in each pool. Firstly, the ATHLET standalone transient simulation is carried out with the time dependent user input boundary, including core power, secondary mass flow rate and temperature. Secondly, the STH transient simulation result is imposed as the time dependent boundary of the CFD model for the hot plenum. The hot pool outlet locates at the same position as the IHX inlet. Thirdly, the mass flow rate and temperature are provided as the inlet boundary of the open loop ATHLET transient case, while the hot pool inlet pressure is used as the pressure boundary of the

patch for the hot pool is divided into 3 rings for inner core, outer core and blanket zone. The hot pool has 4 outlets which are the inlet windows of the 4 operating heat exchangers. As Fig. 12 shows, the 4 inlet boundaries are placed just above the outlet windows of the heat exchangers. The positions of the inlets are specified to get a reasonable velocity field at the outlet windows of the heat exchangers. The inlet temperature of the CFD model is obtained from the ATHLET IHX outlet window. The outlet boundaries of the cold pool are placed at the turning positions of the pumps. The CFD mesh is generated with OpenFOAM snappyHexMesh, which is a mesh generator of OpenFOAM for complicated geometry. The meshes close to the surfaces are refined and the boundary layers are added to the wall boundaries.

3.3. Coupling setup The coupling starts from the ATHLET standalone simulation. In the standalone simulation, both the hot pool and the cold pool are included in the closed loop. The ATHLET standalone steady state simulation

Fig. 11. The inlet panel (left) and outlet panel (right) of the hot pool. 8

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Fig. 12. The inlet panel (left) and outlet panel (right) of the cold pool.

Fig. 13. Nodalization of the primary (left) and secondary (right) systems in Argonne’s SAS4A/SASSYS-1 model.

number of SFR-related analyses at Argonne, and a coupling framework for SAS4A/SASSYS-1/STAR-CCM + analyses has also been demonstrated previously in multiscale SFR transient analyses, including Thomas et al. (2012). These efforts provided a base capability which was expanded over the course of the project to better simulate the Dissymmetric Test conditions. Many modelling decisions were driven by the extent of the coupling framework. All previous work at Argonne had employed single-volume coupling, mostly for evaluation of thermal stratification in the hot pool. This was true also in the blind phase of the Dissymmetric Test, which involved coupling only in the cold pool in order to capture the flow asymmetry there. In the open phase, the coupling framework was expanded to accommodate coupling multiple volumes, allowing for more thorough representation of multi-dimensional flow. However, this still includes only “compressible volumes,” i.e. pools and other large wellmixed volumes, not components such as the IHX.

STH simulation. The difference between the open loop ATHLET model and the ATHLET standalone model is that all the components in the region of hot pool are removed from the ATHLET model for Phénix, and the inlet and the outlet boundaries are added to the two ends of the open loop model. Lastly, the transient simulation result of the ATHLET open loop case is used as the boundary of the CFD model for the cold pool, which ends one complete cycle of the coupling simulation. The bounding surface of the CFD-STH model which is connected to the environment is treated as an adiabatic boundary condition. So the heat loss from the system to the environment is not considered in this simulation. 4. Multiscale model ANL The multiscale model used by Argonne National Laboratory in this work was a coupling between the liquid-metal-cooled reactor safety analysis code system SAS4A/SASSYS-1, Fanning et al. (2017), and the commercial CFD code STAR-CCM+, Siemens (2017). This approach was chosen to leverage a number of prior efforts. Argonne participated previously in the Phénix End of Life Natural Convection Test Benchmark, Thomas et al. (2012) and IAEA (2013), and obtained good agreement with experiment through SAS4A/SASSYS-1 simulation. The SAS4A/SASSYS-1 model from that benchmark was used as the starting point for the model employed here. STAR-CCM + has been used in a

4.1. Sas4A/SASSYS-1 model For the STH model, separate channel models were developed for each of the five flow zones in the core. This was similar to the approach used in the Natural Circulation Test, but with the flow rate and heat generation increased accordingly. The nodalization of the primary and secondary sodium coolant systems are provided in Fig. 13. All volumes 9

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

4.3. Coupling setup

in SAS4A/SASSYS-1 utilize the perfect mixing model for this test. Unlike in prior work, the four IHXs and three PPs are all explicitly represented, instead of being simulated with a multiplicity condition. This allows for the application of the different secondary sodium loop conditions, and simplifies the coupling process. The simplified secondary system model has been developed in order to accommodate the IHX tube-side inlet temperature data that CEA provided for the boundary condition of the secondary system. The IHX inlet temperature is specified by applying time-dependent values for the steam generator outlet, which is immediately upstream of the IHX. The remainder of the secondary system model consists of a pump that draws sodium from an expansion volume and simple connecting pipes with no bends, constrictions, or heat losses. Some differences were present between the STH model used in the stand-alone calculations and the STH model used in the coupled simulations. Notably, in the coupled case, only one volume is used to represent the cold pool. This was done with the intent of simplifying and stabilizing the coupling. In all cases, the core, diagrid, pumps, IHXs, and vessel cooling network were all modelled only with SAS4A/ SASSYS-1.

A file-based coupling scheme is used, with SAS4A/SASSYS-1 being the primary driver. The system is first solved in SAS4A/SASSYS-1, which provides boundary conditions for the CFD modelling, which in turn provides its own output values back to SAS4A/SASSYS-1 as a correction to the standard “well-mixed” volume approach. The complexity of this transient motivated the expansion of the coupling framework to accommodate multiple volumes in the open phase, which proved beneficial in modelling the system and provided improvements over the blind phase results. 5. Multiscale model cea 5.1. Cathare model The CATHARE model, Garcin (2016), of the Phénix reactor comprises the entire primary sodium circuit but only the IHXs with their inlet and outlet pipes for the secondary sodium system. The main reactor vessel components are the core, the hot and cold pools, the primary pumps, the reactor vessel lower plenum, the strongback, the upper-core structure, the diagrid and the shell side of the IHX. The transient IHX secondary inlet flow rates and temperatures measured during the dissymmetric test are used as boundary conditions for the secondary sodium system. In the context of the work reported here, only the primary pumps, the diagrid, the free surfaces and the entire secondary side are actually used as modelled in CATHARE. The rest is overlapped by either TrioMC or TrioCFD.

4.2. Star-CCM + Model For the CFD model, geometry was generated from CAD provided by CEA with the benchmark specifications. Initial blind phase simulations represented only the cold pool in CFD, while open phase simulations also modelled the hot pool. Cold pool coupling was performed at the four IHX outlets, the three PP inlets, and the vessel cooling path outlet. Hot pool coupling was performed at each of the five core outlet channels and the four IHX inlets. A mesh of roughly 0.7 million polyhedral cells was constructed in the open phase, with refinement in the vicinity of the IHX outlets and pump inlets. A finer mesh was also tested that provided little impact on large-scale system behaviour. An extrusion was used at the PP inlet faces (i.e. the domain outlets) to prevent reversed boundary flow. A prismatic cell layer was used at all solid walls. Thermal properties were the same as used in the STH. A URANS approach was used for turbulence, with the Realizable Two-Layer k-ε model of Shih et al. (1994). A segregated solver was used to solve the equations for the system. A slip wall was used to treat the upper surfaces of the pools. All walls were treated as adiabatic in the CFD, with heat losses to the environment modelled in SAS4A/SASSYS-1. To obtain the best initial conditions for the CFD, in spite of the coupling limitations, a steady-state CFD model was developed that included the hot and cold pools as well as the vessel cooling flow, and the heat transfer through the redan and the vessel cooling wall. This provides an initial temperature field that has hotter sodium at the top of the cold pool, which differs from the uniform initial pool temperature used in the blind phase. This may be expected to impact the buoyant mixing in the cold pool during the transient. Given that the IHX internal flow was not modelled, an estimate of the flow pattern at the IHX outlet is necessary. This was performed by including a small upstream-extruded internal IHX volume as the inlet for the steady-state case, which provides the initial conditions for the transient. The flow direction vectors and turbulence quantities on the IHX outlet face were extracted from this steady simulation. The directions of these vectors were applied as a fixed boundary condition at the IHX outlet window in the transient simulation, and held constant throughout the duration, while the velocity magnitudes were provided by SAS4A/SASSYS-1 from the coupling. Although this may impact the results due to lack of stratification in the IHX outlet window, this was deemed a reasonable approximation in this particular transient given that the PPs continue operating as normal, and the primary forced flow was expected to be relatively high and roughly constant over the duration of the transient.

5.2. TrioMC model The core model includes the 456 fuel assemblies of the core. The fuel pins (fissile and fertile), the control rods and the hexagonal cans are modelled in the sub-channel code TrioMC. The code computes the thermal-hydraulics of the sodium flow within a fuel assembly and accounts for the heat conduction within the fuel pins and hexagonal cans. Thanks to MATHYS, which integrates and couples TrioCFD, TrioMC and CATHARE, the heat exchange between the inter-wrapper flow and the hexagonal cans is also taken into account. Correlations are used to account for the pressure drops due to the wire wrap (Cheng and Todreas, 1986) and to compute the heat transfer coefficients. The nodalization comprises 17 meshes in the axial direction. 5.3. TrioCFD model The TrioCFD model includes the hot pool, the inter-wrapper volume, the inter-neutron shielding volume, the neutron shielding elements themselves, the upper core structure, the primary side of the IHXs (porous medium), the cold pool, a portion of the vessel cooling system and a portion of the pump skirts (inlet). The final mesh comprises more than 4 million elements. This mesh size results from a trade-off between the numerical accuracy and the computational cost. A view of the final mesh is offered in Fig. 14. No turbulence model is applied. The Boussinesq approximation is used to account for buoyancy effects. The Gunter and Shaw (1945) correlation is used to predict the pressure losses in the radial direction in the IHXs. Singular pressure loss coefficients are applied at the connection of the vessel cooling with the cold pool to account for the reduction of flow area at this location. Conjugate heat transfer is accounted for: - Between the inner and - Between the inner hot vessel; - Between the inner and - Between the inner and 10

outer sides of the vessel cooling; and outer cold sides of the conical primary outer sides of the primary pumps; outer sides of the hexagonal cans.

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Fig. 14. TrioCFD Mesh of the Phénix Primary Side. Table 1 Summary of the main properties of the multiscale models. Participant

NRG

KIT

Argonne

CEA

STH code CFD code Coupling method

SPECTRA CFX Explicit; Domain overlap Two-way

ATHLET OpenFOAM Explicit; Domain overlap One-way

SAS4A/SASSYS-1 STAR-CCM+ Explicit; Domain overlap Two-way

# nodes STH # cells CFD (million) Turbulence modelling

407 1.17, hex and tets SST k-ω Constant Prt = 2.0 Automatic wall function Porous medium with heat sink Yes, walls and roof in STH Included in CFD

450 4.1, hex-dominated SST k-ω

699 0.7, polyhedral Two-Layer Realizable k-ε Constant Prt = 2.0 Two-Layer STH

CATHARE TRUST (TrioCFD) Explicit Domain overlap (STH-CFD), domain decomposition (sub-channel – CFD) Two-way 2281 4.02 tets only None None Porous medium coupled with STH secondary side Yes (external vessel in CFD and roof in STH)

Wall treatment IHX model CFD Heat loss to environment Conjugate heat transfer through vessels

Wall function Simulated with STH, Not included in CFD model No, adiabatic wall boundary Not included in CFD

Yes, walls and roof in STH CHT only through redan in CFD

Included in CFD

Thermal-hydraulic coupling is applied at:

No conjugate heat transfer was applied at the primary side of the IHXs at the outer shell. Heat losses were accounted for at the outermost shell of the vessel. After the blind benchmark exercise, issues were identified at the outlet of the core: during this phase, the core was modelled by CATHARE and the sodium exiting the fertile zone was “stuck” at the core outlet because the inter-wrapper volume was not accounted for. This motivated the addition of TrioMC in the current model which significantly improved the model in this regard. Several improvements were also made in the modelling of heat losses and conjugated heat transfer.

- The pump inlets (a pressure boundary condition is specified at the inlet of pump 2 – the one located between the DOTE dummy IHXs); - The upper interface of the hot pool (in order to avoid high computational costs, no interface tracking is performed in the CFD calculation: instead, the CFD domain ends below the lowest level of the free surface. The STH code is used to compute the elevation of the free surface). The same approach is applied for the vessel cooling and the cold pool upper interface; - The core outlet (CATHARE provide the core inlet mass flow rate and enthalpy to TrioMC and TrioMC the core outlet mass flow rate and temperature to TrioCFD); - The core vessel cooling inlet.

5.4. Coupling setup The calculation is first started by running CATHARE in stand-alone mode so that it can reach a stable steady-state. Then, TrioCFD and TrioMC are started. At first, only the thermal feedback at the IHXs is activated: CATHARE provides the boundary conditions to TrioMC and TrioCFD but only the energy transferred at the IHXs is provided by TrioCFD to CATHARE. Finally, the full coupling is started and the calculation runs normally. A summary of the coupling phases is given below:

Thermal coupling is applied at the IHXs to account for the heat exchange between the primary and the secondary side. The current model requires more than 3 days to run the transient using 280 CPUs. 6. Multiscale models The choices that are made by each benchmark participant for developing the multiscale model of the Phénix reactor are summarized in Table 1.

- 17–7000 s: CATHARE runs in stand-alone mode; - 7000 to 7100 s: TrioCFD and TrioMC are started using the boundary conditions from CATHARE. The thermal feedback from TrioCFD is active; - From 7100 s on, CATHARE, TrioMC and TrioCFD run and exchange data at each time step; - The transient is initiated at t = 10 000 s.

7. The dissymmetric test The events in the dissymmetric test are summarized in Table 2: trip of one of the secondary pumps, reactor shutdown followed by a scram. 11

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Table 2 Events occurring in the dissymmetric test. Time (s)

Event

0 5

Secondary pump trip (on one loop): speed reduced from 700 to 100 rpm in about 13 s. The 100 rpm correspond to the rotational speed of the auxiliary pump rotor Automatic reactor shutdown: insertion of the control rods (1.4 mm/s) in 45 s Turbine trip Secondary pump speed reduced (on the other loop) from 700 to 110 rpm in about 60 s Scram End of benchmark test

48 1800

the IHX outflow is going to be redirected toward the top of the cold pool as illustrated in Fig. 16. This is due to the density difference of the hot sodium from the IHX entering the relatively colder cold pool. The timing of the redirection of the flow is one of the challenges of the benchmark. Besides, the flow rate in the second secondary loop decreases at a slower rate (within about one minute). Therefore, the hot shock at the outlet of the corresponding IHXs is delayed and of a lesser amplitude. This triggers the dissymmetric patterns in the reactor. The power extracted by the IHXs is dissymmetric during the first phase of the transient. In the second part of the benchmark, after the events have occurred, a new equilibrium will be reached. In this part, the models are challenged on the mixing of the sodium in the cold pool, with two sets of active IHXs and one set of inactive DOTEs. Furthermore, the modelling of the heat losses through structures and towards the environment are key parameters in this long term behaviour of the Phénix reactor. 8. Results In the dissymmetric transient, the main events occur in the first 50 s. The majority of the consequences of these events occurs in the first 200 s of this benchmark. Therefore, the discussion of the results focusses on the first 200 s and ends with a brief description of the modelling of the remaining 1600 s of the transient. The end-of-life dissymmetry test data mainly consist of temperature measurements at several locations in the pools and the PPs, so the presented results are focused on comparing computed temperatures amongst the participants with the experimental results. The assumed added value of multiscale simulations is to compute the three-dimensional flow pattern and temperature distribution in the relevant regions with higher accuracy than one-dimensional codes. Due to the steady primary forced circulation, a direct jet flow of sodium from the core outlet to the IHXs inlet windows is observed. Several thermocouples are installed at the inlet windows and outlet windows of the IHXs to measure average inlet and outlet temperatures. The mass flow averaged temperatures as computed by the multiscale calculations are compared to these data in Figs. 17 and 18. The averaged inlet temperature in Fig. 17 shows a reasonably accurate prediction for all participants. The results computed by NRG and CEA resolve more temperature fluctuations. This might be due to the modelling choice of including the complete primary side of the IHXs in the CFD domain, resulting in a more realistic flow path towards and inside the IHXs and more direct feedback between the secondary side in the STH code and the primary side in the CFD code. ANL starts with an accurate estimate of the temperature but begins to underestimate after 50 s. From 130 s onwards, NRG’s and KIT’s results are in good agreement with the measurement data. The fast flow rate decrease in the secondary loop 1 leads to a hot shock at the outlet of the primary side of the IHXs K and M which propagates to the cold pool. This is seen in Fig. 18, where the average outlet temperature of IHX M shows a sharp temperature rise at a few seconds after the start of the transient. This rise is reasonably computed by all participants, though all results show a faster rise than the measurements. The temperature peak is underestimated by CEA, well

Fig. 15. Thermal flow path (nominal).

The following data describing these events are provided as a function of time as input to the participants: flow rate and inlet temperature of the secondary loops; core power. Looking at the global phenomena in the dissymmetric test, the events lead to a loss of cooling in two IHXs, later followed by the other two IHXs and a reduction of the core power. The loss of cooling performance of the IHXs influences the flow patterns in the cold pool as sketched in Figs. 15 and 16. As can be seen in Fig. 15, in normal operation conditions, the sodium from the IHXs flows down to the bottom of the PP skirt. This is still true at the beginning of the transient. Slowly,

Fig. 16. Thermal flow path (transient). 12

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Fig. 17. Average inlet temperature of the primary sodium in IHX J.

though they all overestimate the temperature rise. The temperature increase for ANL is slightly steeper and occurs sooner, while KIT predicts a steep temperature rise but it is delayed. The difference between the timing and steepness of the slope, which is most accurately computed by CEA and NRG, might have the same origin as the difference in the predicted fluctuations in the IHX inlet and outlet temperatures, namely that CEA and NRG included the complete IHXs in their CFD model, while KIT and ANL did not. This could lead to a difference in outflow direction and temperature distribution in the redundant sodium. Since the outflowing jets from the IHXs are partly directed towards the PP skirts, this directly influences the temperature computed at the top of the PP. On the one hand, hot sodium leaving the IHX outlet window could directly enter the PP skirt. On the other hand, the jet could impinge on the PP skirt, heating the sodium in the skirt by conduction. Whether the jet impinges on or enters into the PP skirt is determined by the balance between the momentum of the jet and the buoyancy forces between the cold sodium present in the cold pool and the hot sodium leaving the IHX. Furthermore, the thermal initial condition in the CFD model of KIT assuming that the pools are at uniform temperature result in more cold sodium in the hot pool, which contributes to the later temperature rise at the top of the PP skirt. PP2 is surrounded by 2 DOTE’s which are passive components. Therefore the flow near PP2 is less disturbed and the temperature in its surroundings is more uniform. The temperature at the top of PP2 (Fig. 20) shows a single smooth slope which is nicely followed by NRG,

predicted by NRG and ANL and overestimated by KIT. After the peak, KIT is closest to the experimental results, while ANL, CEA and NRG underestimate the temperature. Concluding the modelling of the IHXs for the first 200 s, all participants predict the trend quite well and predict the temperature within 5–30 °C of the measurements. From the IHX outlet, the sodium flows into the cold pool and into the skirts of the PPs. The initial temperature at the top of the skirt of PP1 is well modelled by all participants (Fig. 19). The hot sodium from IHXs K and M, due to the triping of the secondary pump, first reaches PP1 after about 25 s (see Fig. 19). The temperature then drops shortly at t = 65 s and subsequently increases again. This phenomenon is not observed by PP2 (see Fig. 20) and PP3 (see Fig. 22) which only see a temperature increase. Since PP2 is located further away from the 4 IHXs and it is surrounded by two passive components (DOTEs), it can only suck hot sodium after it mixed in the cold pool. This explains the delayed temperature increase of this pump. The initial small temperature rise observed at the inlet of PP1 it is not observed at the inlet of PP3 either because the secondary flow rate in the adjacent IHXs is decreased more slowly. The sodium at the pump skirt level has enough time to get mixed with sodium from the cold pool which leads to a steady temperature increase. Coming back to Fig. 19, a similar first temperature rise in PP1 is resolved by all participants (Fig. 19), though the timing and peak temperature differ. Between 60 s and 150 s a smooth temperature increase is measured, more or less resolved by ANL, CEA and NRG,

Fig. 18. Average outlet temperature of the primary sodium in IHX M. 13

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Fig. 19. Average temperature at the top of the skirt of PP1 during the first 200 s of the dissymmetric transient.

Fig. 20. Average temperature at the top of the skirt of PP2 during the first 200 s of the dissymmetric transient.

Fig. 21. Average temperature at the inlet if IHX J during the complete dissymmetric transient.

Proceeding to the second part of the transient, it should be noted that the two secondary cooling loops equally contribute to the reactor cooling in the long term, as soon as the secondary pumps reach an equal mass flow rate. From there on, symmetry is established again in the

ANL and CEA, though the results from CEA and ANL show an additional peak near 100 s. The results obtained by KIT show a delay of the temperature rise by about 35 s, which is a similar delay as observed in Fig. 19 for PP1. 14

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Fig. 22. Average temperature at the top of PP3 during the complete dissymmetric transient.

entire reactor. The hot and cold pool temperature then steadily decrease at a rate of 90 K/h (see Figs. 21 and 22). The modelling achievements for this long term behaviour of the Phénix reactor are assessed looking at the IHX inlet temperature (Fig. 21) and the temperature at the top of PP3 (Fig. 22). The fact that the temperature decreases at both locations is predicted by all partners. The results from KIT show a very good agreement in the long term. While the temperature gradient at the IHX inlet window is already in good agreement with the measurement results, the temperature at the top of PP3 almost overlaps the experimental data. This might however be a fortunate result due to the different initial enthalpy present in the system caused by a different assumption for the initial temperature. ANL, CEA and NRG underestimate the cooling rate. The temperature computed by ANL is lower than the temperature from NRG and CEA, being closer to the measurements, especially at the top of PP3. The cooling rate predicted by ANL and NRG is very similar, while CEA computes a slightly lower cooling rate.

• • • •

Acknowledgements This work is part of the SESAME European project which has received funding from the Euratom research and training program 20142018 under grant agreements No. 654935. All the work of NRG described in this paper was funded by the Dutch Ministry of Economic Affairs. Argonne work was performed under an International Nuclear Energy Research Initiative (INERI) collaboration. Argonne acknowledges support from the U.S. Department of Energy Office of Nuclear Energy, Office of Advanced Reactor Technology. Furthermore, the authors gratefully acknowledge the contributions of all colleagues involved in this project.

9. Conclusions and recommendations Because of the complexity of the phenomena involved and of the size of the physical domain, the modelling of the Phénix reactor proved to be a challenging task. One has to constantly seek the best compromise between the accuracy and the computational cost which is prone to quickly escalate. The results show that the challenge actually consists of two challenges, namely correctly computing the thermal-hydraulics of the first three minutes of this transient and finding the correct parameters to accurately compute the remaining 27 min. For the short term, it is found that the IHXs should be included in the CFD model in order to correctly compute the momentum and stratification of the sodium leaving the IHXs. For the long term, three of the four participants underestimate the cooling rate, so deeper investigation of the heat losses from and the thermal inertia in the Phénix reactor. Several recommendations have been formulated for potential improvements of the Phénix multiscale models:

References ANSYS CFX 15.0 user guide, 2013. Cheng, Shih-Kuei, Todreas, Neil E., 1986. Hydrodynamic models and correlations for bare and wire-wrapped hexagonal rod bundles—bundle friction factors, subchannel friction factors and mixing parameters. Nucl. Eng. Des. 92 (2), 227–251. Churchill, S. W., Bernstein, M., 1977, Journal of heat Transfer, 99 p. 300, 1977 (see F.P. Incropera, D.P. DeWitt, “Fundamentals of Heat and Mass Transfer”, 4th Edition, 1996, Section 7.4 eq. 7.57). Fanning T.H., Brunett A.J., Sumner T., eds., 2017. “The SAS4A/SASSYS-1 Safety Analysis Code System: User's Guide,“ Argonne National Laboratory, ANL/NE-16/19. J.K. Fink L. Leibowitz Thermodynamic and Transport Properties of Sodium Liquid and Vapor, Argonne National Laboratory 1995 ANL/RE-95/2. Garcin, J. CATHARE 2 V25_3mod6.1 code: General description, DEN/DANS/DM2S/ STMF/LMES/NT/16-033/A, 2016. Guidez, J., 2013. Phénix: the experience feedback. EDP Sciences. Gunter, A. Y., and W. A. Shaw., A general correlation of friction factors for various types of surfaces in cross flow, Trans. ASME 67.8 (1945): pp. 643–660. IAEA (2013), IAEA TECDOC No. 1703: Benchmark Analyses on the Natural Circulation Test Performed during the PHENIX End-of-life Experiments. Gunter, A.Y., Shaw, W.A., 1945. A general correlation of friction factors for various types of surfaces in cross flow. Trans. ASME 67.8 643–660 In press.

• A CFD study of the IHX comprising the primary and secondary side • •

through the walls between the various flow domains or including more structures to account for thermal inertia; The incorporation of a more detailed partitioning of the pool geometries in the STH; Further analysis of the flow distribution between the core and the vessel cooling, which has a direct impact on the flow through the IHXs; In the long term, it could be interesting to model the diagrid in CFD. In this way, the dissymmetric flow due to the different pump inlet temperatures could be propagated through the core; The initial temperature distribution in the Phénix reactor at the start of the transient should be assessed in more detail to understand the good results from KIT for the long term behaviour.

might also be helpful to better understand the behaviour of the latter at different flow conditions and to check that the heat transfer is correctly predicted by the coupled codes.; Assessment of mesh refinement in the CFD domain, though CEA observed limited effect on the overall behaviour; The long term cooling rate might be affected by more heat loss to the environment, other heat loss mechanisms, such as higher IHX power removal, more directly accounting for conjugate heat transfer 15

Nuclear Engineering and Design xxx (xxxx) xxxx

H.J. Uitslag-Doolaard, et al.

Description, Volume 2 – User’s Guide, Volume 3 – Subroutine Description, Volume 4 Verification and Validation, NRG report K5024/10.101640, Arnhem, April 24, 2009. Thomas J.W., et al., 2012. “Validation of the Integration of CFD and SAS4A/SASSYS-1: Analysis of the EBR-II Shutdown Heat Removal Test,” ICAPP '12, Chicago, USA. Thomas J.W., Fanning T.H., Dunn F.E., et al., 2012. “Analysis of the Phénix End-of-Life Natural Convection Test with SAS4A/SASSYS-1,” ICAPP '12, Chicago, USA. Uitslag-Doolaard, H.J., Alcaro, F., Roelofs, F., Zwijsen, K., 2018, System thermal-hydraulic and multiscale simulations of the dissymmetric transient in the Phénix reactor, ICAPP 2018.

Lyon, R.N., 1951. Liquid metal heat transfer coefficients. Chem. Engng. Progr. 47 (2), 75–79. Mikityuk, K., 2009. Heat transfer to liquid metal: review of data and correlations for tube bundles. Nucl. Eng. Des. 239, 680–687. Shih T.-H., et al., 1994. “A New k-ε Eddy Viscosity Model for High Reynolds Number Turbulent Flows – Model Development and Validation”, NASA TM 106721. STAR-CCM+, 2017, (C) CD-adapco, LTD (A Siemens Company). Available: https://mdx. plm.automation.siemens.com/star-ccm-plus. Stempniewicz, M.M., 2009, SPECTRA Sophisticated Plant Evaluation Code for ThermalHydraulic Response Assessment, Version 3.60, August 2009, Volume 1 – Program

16