MULTIPHYSICS MODELISATION OF AN UNPROTECTED LOSS OF FLOW TRANSIENT IN A SODIUM COOLED FAST REACTORS USING A NEUTRONIC-THERMAL-HYDRAULIC COUPLING SCHEME

Sodium cooled fast neutron reactors (SFR) are one of the selected reactor concepts in the framework of the Generation IV International Forum. In this concept, unprotected loss of cooling flow transients (ULOF), for which the non-triggering of backup systems is postulated, are regarded as potential initiators of core melting accidents. During an ULOF transient, spatial distributions of fuel, structure and sodium temperatures are affected by the core cooling flow decrease, which will modify the spatial and energy distribution of neutron in the core due to the spatial competition of neutron feedback effects. As no backup systems are triggered, sodium may reach its boiling temperature at some point, leading to local sodium density variations and making the transient fluctuate in a two-phase flow physics where thermalhydraulics and neutronics may interact with each other. The transient phenomenology involves several physic disciplines at different time and spatial scales, such as core neutronics, coolant thermal-hydraulics and fuel thermo-mechanics. This paper presents the results of thermalhydraulic/neutronic coupled simulations of an ULOF transient on the SFR project ASTRID. These coupled calculations are based on the supervisor platform SALOME to link the neutron code APOLLO3® to the system thermal-hydraulic code CATHARE3. The physical approach used by the coupling to describe the neutron kinetic is a quasi-static adiabatic one, updating the normalized spatial power distribution periodically by performing static neutron calculations, while a point kinetic model associated to a neutron feedback model calculates the power amplitude variations.


INTRODUCTION
The Gen-IV design approach for future reactor is based on improving safety, managing natural resources and optimizing economic performance. With respect to safety of sodium cooled fast reactor, the gen-IV French SFR prototype ASTRID-600MWe proposes advanced features such as an in-vessel core catcher, an innovative nitrogen Brayton power conversion cycle and a low sodium void effect core design [1]. The core geometry is specifically designed to reduce the overall sodium void reactivity effect by increasing neutron leakages in the upper part of the core. The sodium feedback effect is a key parameter to improve the core safety behavior during transients involving sodium concentration variations. During an unprotected loss of flow transient, the core cooling flow decreases in few seconds without any reactor scram. As no core backup-systems are triggered by assumption, the reactor behavior during the transient is entirely driven by the spatial competition of its neutron feedback effects. Sodium may reach at some point its boiling temperature, paving the way to the establishment of a two-phase flow dynamics with large sodium vapor slugs interacting with neutron physics through the sodium feedback effect. By itself two-phase flow dynamics in heating tubes can generate, under certain circumstances, dynamic instabilities like chugging phenomena which cause periodic burst of vapor followed by condensation phenomena [2]. The sodium void neutron feedback can amplify or mitigate these phenomena following the situation. In order to study the dynamic aspects of the coupling between neutron physics and two-phase flow during an ULOF transient, an adiabatic kinetic model [3] is proposed completed by a coupling scheme based on the neutron code APOLLO3® [4], used to update periodically the spatial power distribution, and the system thermalhydraulics code CATHARE3 [5]. The results obtained by the use of the coupling platform will be compared to those of decoupled calculations with constant power profile.

DESCRIPTION OF THE ASTRID REACTOR
ASTRID is a SFR project supported by the French Alternative Energies and Atomic Energy Commission (CEA) and its industrial partners in the framework of the Generation IV international forum [1]. It is a pooltype reactor with a thermal power of 1500 MWth (600 MWe). The reactor core is located inside a primary pool and is continuously cooled by a liquid sodium flow produced by three mechanical pumps (cf. Figure  1). The energy produced by nuclear fissions is transferred to the secondary loops through intermediate heat exchangers (IHX) immersed inside the primary vessel and then to an innovative nitrogen Brayton power conversion cycle through sodium/gaz heat exchangers (SGHE). The core consists in a lattice of 288 hexagonal sub-assemblies containing 217 fuel pins filled with a mix of uranium and plutonium oxide (MOX). Reactor control is ensured by the mean of two diversified control rod systems containing boron carbide rods with various enrichment in 10 B (nominal and emergency system). Reflector and absorbing media wrap the driver core and help to decrease neutron leakage and to limit irradiation damages on main primary components. The inherent low sodium void effect is achieved by maximizing neutron leakage in case of sodium draining in the upper part of the core. This is done by a combination of technological solutions such as the addition of an upper sodium plenum located directly above fissile pins, the insertion of a fertile zone inside inner core and finally a taller outer core than inner one (see Figure 2).

Theoretical Aspects of Neutron Kinetics
In order to resolve the neutron population evolution in the reactor during the transient, we have to solve the set of coupled Boltzmann (Eq. 1) and kinetic precursor equations (Eq. 2). The first equation describes the spatial, temporal and energetic evolution of the neutron flux inside the reactor (Φ(‫ݎ‬ ⃗, Ω ሬሬ⃗ , ‫,ܧ‬ ‫,))ݐ‬ whereas the second one describes the evolution of delayed neutron families. The operators ‫,ܮ‬ ‫,ܪ‬ ‫ܨ‬ , ܳ takes respectively into account neutron losses by leakage and absorption (Eq. 3), neutron arrivals due to scattering process (Eq. 4) and finally prompt (Eq. 5) and delayed neutrons arrivals (Eq. 6) induced by fission. Note that Σ ௧ , Σ ௦ , Σ correspond to the total, the scattering and the fission macroscopic cross sections, ߚ, ߥ, Χ to the total delayed neutron fraction, the neutron yield and the prompt neutron spectrum and finally Χ ௗ , ߣ , ߚ , ‫ܥ‬ to the delayed neutron spectrums, the decay constants, the delayed neutron fractions and the concentrations of the delayed neutron families.
The direct resolution of equations (1) and (2) can be done and refers to the spatial kinetic method. The quasi-static methodology corresponds to a rewriting of these equations introducing a factorization of the time dependent neutron flux into the product of an amplitude and a shape function (Eq. 7). A set of four equations, which are equivalent to the (1) and (2), can then be obtained. The first one deals with the amplitude function (Eq. 8 and 9) and is also known as point kinetic equations. The second one deals with the shape function (Eq. 10 and 11). New effective parameters (ߚ ෨ , Λ ෩ , ‫ܥ‬ ෩ ) completed by the dynamic reactivity of the system(ߩ) are introduced. They formally depend of the shape function. In order to solve these sets of differential equations, various level of assumptions can be introduced with respect to shape function processing (see Table 1 and ref. [6]). The legacy Quasi-Static method (QS) solves the kinetic problem by assuming that the shape function varies slowly compared to the amplitude function ( ). The Improved Quasi static Method (IQM) improves the time processing of the shape function and tries to solve the kinetic problem with the fewest approximation. These two methods require to perform inhomogeneous source calculations to solve the equation (10). On the contrary the adiabatic method is a simpler model in which all the temporal derivate terms in equations for the shape function (Eq. 10 and 11) are neglected. It implies that the multiplication factor should remain close to 1.0 during the transient, that the precursor distribution is instantaneously in equilibrium with respect to the shape and finally that the impact of the temporal inertia terms on the shape is negligible. Shape evaluations along the transient are then reduced to EPJ Web of Conferences 247, 07001 (2021) PHYSOR2020 https://doi.org/10.1051/epjconf/202124707001 a suite of stationary flux calculations without inhomogeneous source calculations. The adiabatic method was chosen as a first step for the kinetic method in the coupling scheme to treat ULOF transient due to its pretty simple integration in the calculation scheme. This model do not suit well for fast transients (~100$/s) in thermal reactors as stated in the reference [6] due to the elimination of both flux and precursor spatial inertia terms. However studied ULOF transient involves weak reactivity ramp (reactivity ramp < 0.3 $/s) which guaranties a multiplication factor close to one. Moreover, as the order of magnitude of the averaged neutron velocity (v) in a SFR is a thousand time upper than in thermal reactor, the flux inertia assumption ( Nevertheless, more robust validations based on comparisons with improved quasi static or even spatial kinetic models, have to be performed to assess the validity of the adiabatic assumptions in our case [3] [7].

Feedback Coefficient Model
Thermal feedback effects in a reactor comes from the dependencies of the nuclear reaction rate probabilities to thermal parameters such as material temperatures and densities, control rod relative insertion and irradiation state of the fuel. The dynamic reactivity of the kinetic system can be computed at each time step using the feedback coefficient approach based on the following major feedback effects in SFRs (Eq. 12). Note that most of the feedbacks have a spatial dependence which is taken into account using spatial distributions of feedback coefficients obtained using the standard neutron perturbation theory. ߜߩ = ߜߩ ே + ߜߩ + ߜߩ ி௨ + ߜߩ ௗ + ߜߩ ௐ + ߜߩ ோ x The sodium expansion effect is a combination of two opposing effects: a negative reactivity insertion due to an increase of the neutron boundary leakage and a positive one due to a neutron spectrum hardening and a decrease of neutron capture rates by the sodium. The result in terms of reactivity depends on the localization of the sodium expansion in the core. x The Doppler effect opposes any fuel temperature increase by strengthening 238 U sterile capture.
x The fuel axial expansion opposes any temperature increase by decreasing fuel concentrations in the central part of the core. The driving temperature is the fuel temperature itself for fresh fuel or the clad one for irradiated fuel, for which fuel pellets are closely tied to the clad by thermal and irradiation effects.
x The clad and wrapper axial expansion reinforces any temperature increase by decreasing amount of steel sterile neutron capture in the central part of the core.
x The Core-Vessel-Control Rod differential expansion (CVCR) opposes any reactor temperature increase by triggering control rod insertion in the core due to the competition of relative expansion

Description of the Calculations Models
The coupling algorithm is based on a direct transfer of relevant physics parameters between the neutron code, APOLLO3® [4], and the system thermal-hydraulic code, CATHARE 3 [5]. The coupling platform, SALOME [8], supervise the adiabatic kinetic time loop, data transfers and mesh projection operations. Punctual kinetic equations are resolved by CATHARE3 using the feedback coefficients model to compute the global dynamic reactivity of the system at each time step whereas the normalized power shapes are periodically updated by resolving the stationary version of the Boltzmann equation with APOLLO3®. The reactivity variation due to sodium density variation and Doppler Effect, computed by APOLLO3® between two updates, is moreover used to correct the sodium and Doppler reactivity contributions estimated by the feedback model. As neutron flux calculations are costly in terms of calculation time, the supervising algorithm tries to limit as possible the update needs. Update criteria are implemented to fulfill this objective, such as the global reactor power variation (>5%) and the local sodium density variation (>0.02 g.cm -3 ). The values of these criteria come from specific sensitivity studies.

Neutronic model of the core
ASTRID nominal state parameters, such as isotopic compositions, feedback coefficients, kinetic data and initial power distributions, were estimated using the diffusion solver available in the ECCO/ERANOS neutron code [9] with the JEFF 3.1 nuclear library. These initial data were provided as constants to CATHARE3 for its point kinetic model. An additional set of cell calculations were performed with ECCO to create a parametric cross section database with respect to fuel temperature and sodium density. APOLLO3® [4] was used during the coupling calculation for the periodical update steps. During such a step, sodium density and fuel temperature distributions calculated by CATHARE3 are transferred to neutronic solvers using mesh projection operations. APOLLO3® evaluates then the updated cross sections by interpolation in the parametric cross section databases and performs 3D flux calculations using either the MINARET SN transport solver [4] with a S4 angular scheme and linear discontinuous Galerkin spatial finite elements or the MINOS solver [4] with diffusion approximation and linear Raviart-Thomas finite element. MINARET calculations capitalize on the one third core symmetry to reduce the computation time.
Normalized power maps and reactivity variations are then transferred to CATHARE3 to update power maps and correct sodium and Doppler reactivity contributions estimated by CATHARE3 with its feedback model.

Thermal-hydraulic model of the core
The core was modeled by a set of parallel hydraulic pipes (one per sub-assembly) connected to the diagrid and the hot plenum. The one third core symmetry allowed to reduce the number of pipes to 96. The whole primary loop was modeled (hot and cold plenums, IHX, primary pumps, diagrid and strongback). Secondary loops were described using static temperature and pressure boundary conditions. The core thermal-hydraulic was solved by CATHARE3 using a 6 equations models associated with a specific revision of the sodium two phase flow closure laws (revision 3, [12]]), which improves the predictive capability of the code in the field of sodium two phase flow dynamics. Fuel temperatures were evaluated for each sub-assembly by resolving the time-space dependent heat equation for an averaged fuel pin whose power had a time dependence according to the neutron kinetic of the transient. The correlation laws used in CATHARE3 for the fuel thermal conductivity and heat capacity are respectively the Philipponneau's law, implemented in the reference code for SFR thermomechanic under irradiation studies GERMINAL V2.2 [10], and the Carbajo's one [11]. The thermal conductance of the gap between clad and pellets depends strongly on the fuel irradiation state. GERMINAL V2.2 was used to compute the gap conductance values which were assumed to be constant during the transient. Finally as the calculation model was not defined to treat with material melting and fuel displacements, the calculation stopped as soon as clad material reached its melting point.

ULOF CALCULATION RESULTS
Before starting the ULOF transient, a careful care was done to define the initial reactor state. A convergence step between APOLLO3® and CATHARE3 calculations was done in order to insure a coherent and stabilized initial maps of sodium and fuel temperature. Figure 3 presents the converged initial reactor power map. Moreover to validate the transfer to CATHARE3 of the thermomechanical properties computed by GERMINAL V2.2, a specific comparison of averaged fuel temperatures has been done (cf. figure 4). The unprotected loss of flow event corresponds to the loss of all the primary pumps without triggering of the reactor emergency system. During such an event, the active core cooling flow decrease is driven by the mechanical inertia of the primary pumps (halving time of 11.0s) and fades away in few seconds (see Figure  5). The sodium temperature at the core outlet quickly increases whereas the inlet core temperature is maintain constant (400°C) due to the operation of IHXs. The sodium overheating, mainly located in the upper part of the core, induces sodium expansion which drives negative reactivity contributions due to firstly the negative sodium void effect, secondly the Core/Vessel/Control Rod (CVCR) differential expansion and finally to fuel axial expansion induced by clad expansion (cf. Figure 6). These negative feedbacks drive a slow global power decrease whereas the sodium quickly overheats until it reaches its saturation temperature (~900°C) which triggers the start of a new phase of the accident with a two-phase flow dynamics interacting with neutron physics. Two-phase flow dynamics in a heating tube can leads to the development of dynamic instabilities like chugging phenomena which cause periodic burst of vapor followed by condensation phenomena [2]. These phenomena are mainly driven by pressure drops dynamics on the two-phase flow featuring large sodium vapor slugs, and result from the competition between floatability friction and acceleration forces. The sodium void effect can enhanced or mitigate these instabilities following the situation. Sodium boiling in the upper part of the core induces a negative reactivity insertion, which will decrease power and hasten vapor condensation, which will then induce positive reactivity insertion and so on (see Figure 7). The figure 7 compares the total core reactivity estimated along the transient by point and adiabatic kinetic models. As expected, the boiling onset is reached earlier with a neutron transport solver (MINARET) than with neutron diffusion solvers (MINOS / ERANOS) as these latter overestimate sodium void leakage effect. The reactor power at the boiling onset is consequently affected with a deviation of 7% (see Figure 5). After the boiling onset, adiabatic kinetic models brings more reactivity variation than the point kinetic model. This can be explained by looking at the figure 8 which presents the power deviation map between boiling and initial states for the adiabatic kinetic model with MINARET. Local power depressions around voided assembly head can be observed (~- 10%) which leads to stronger local sodium temperature decreases which will enhance reactivity oscillations. Note that this effect is more pronounced with MINOS as a consequence of the sodium void effect overestimation. These results militate for the use of a specific kinetic treatment in order to correctly account for spatial flux deformation along ULOF transient when boiling occurs as it seems to enhance reactivity oscillations. Finally the boiling front evolution is of major importance with respect to the core degradation, as its penetration in the fissile part of the core will trigger clad degradations by loss of cooling (see Figure  9). The adiabatic MINOS calculation, which overestimates the sodium void effect, is the only simulation which avoid clad melting during the first 150 seconds of the transient. The results, presented here, are sensitive not only to the feedback and kinetic models, as previously seen, but also to thermo-hydraulics and thermomechanical chosen models. No major conclusions concerning the core degradation itself can be inferred apart from the necessity to conduct robust uncertainty and numerical bias analysis.

CONCLUSIONS
During an ULOF in a SFR the core cooling flow decreases in few seconds without any reactor scram leaving the reactor in a power state mode with limited cooling capacities. The reactor behavior is then driven by the spatial competition between neutron feedback effects, especially the sodium void neutron feedback. The induced sodium temperature increase leads in a few seconds to the start of sodium boiling phenomena paving the way to two-phase flow dynamics in interaction with neutron physics, especially in ASTRID, whose core was designed to enhance negative sodium feedback effects in the upper part of the core. The establishment of stable boiling front in the upper part of the core is of major importance with respect to the prevention of core degradation. In order to study the dynamic aspects of the coupling between neutron physics and two-phase flow, an adiabatic kinetic model was proposed, completed by a coupling scheme based on the neutron code APOLLO3®, used to update periodically the spatial power distribution, and the system thermal hydraulics code CATHARE3. A detailed sub-assembly by sub-assembly core model was used in both codes as long as a realistic description of the isotopic compositions and thermal properties of each fuel sub-assemblies. Both point kinetic and adiabatic simulations, with transport (MINARET) and diffusion (MINOS) solver updates, were performed. Diffusion solvers overestimate the sodium void effect in the plenum which explain a reduced power at the boiling onset and stronger reactivity oscillations afterwards. During the boiling phase, adiabatic models enhances reactivity oscillations with regard to the point kinetic model thanks to the local power deformations which strengthen chugging effects, which can be observed in both adiabatic and point kinetic simulations. These results militate for the use of a specific kinetic treatment in order to correctly account for spatial flux deformation along ULOF transient when boiling occurs. Finally, the adiabatic kinetic model was used as a first step but further validation studies have to be performed with more robust kinetic models (IQS model, Spatial Kinetic Model).