16 O 1.952e-02 23 Na 6.352e-03 56 Fe 1.861e-02 235 U 1.542e-05 238 U 7.599e-03 238 Pu 5.833e-05 239 Pu 1.238e-03 240 Pu 5.773e-04 241 Pu 1.617e-04 242 Pu 1.743e-04 241 Am 2.71e-05 10 B 6.388e-03 11 B 2.587e-02 12 C 8.065e-03 23 Na 1.094e-02 56 Fe 1.256e

The purpose of this paper is to study and compare different neutronic approaches used to calculate ULOF transients in sodium cooled fast reactors. A first objective is to compare two codes used for spatial kinetics calculations. Two neutronic models are compared. The first one is based on a deterministic (discrete ordinate S N ) approach, using tabulated self-shielded cross sections, where the core reactivity and the power shape distribution are evaluated at each time step of the transient calculation. The second model relies on the Transient Fission Matrix (TFM) approach, condensing the response of a Monte Carlo neutronic code in time dependent Green functions characterizing the local transport in the reactor. This second approach allows a fast estimation of the reactivity and of the flux redistribution in the system during the transient with a precision closed to that of the Monte Carlo code. Both models have been coupled to the thermalhydraulics and applied on an ASTRID representative assembly. This application case is supposed to be sensitive to power redistributions. A second comparison between spatial kinetics and point kinetics calculations has been led to study this point. Finally we obtain a good agreement between spatial and point kinetics on the ULOF calculation, while some discrepancies are observed between the TFM and the S N approaches on the power level stabilization, due to difference on the feedback estimation in both models.

I. INTRODUCTION

Studying reactor behavior during various transient situations requires a coupling between neutronics and other physics such as thermalhydraulics. This coupling acts through the feedback effects linking the temperature (Doppler), the density or the geometry deformation to the neutron behavior in the core. These feedback effects induce an evolution of the total power value and on its distribution in the core. Various approaches exist to model the neutron kinetics depending on the required precision, on the studied system complexity and on the computational resources accessible. In order to enhance the precision of the neutronics modeling, many developments concern the consideration of spatial flux redistributions during the transient. In the frame of the study of sodium fast reactors, the capability to predict the spatial decoupling is an important point due to the geometrical heterogeneities of core concepts such as the low void concept ASTRID (Ref. 1). This kind of reactor represented in Fig. 1 is composed of both axial and radial heterogeneities. In this frame, a comparison of neutronics models and of calculation tools is performed between the French Alternative Energies and Atomic Energy Commission (CEA) and the Israel Atomic Energy Commission (IAEC). These studies aim to define the validity domain of point kinetics models for system code applications.

Fig. 1. Axial cut of the reference core of ASTRID, with the fuel in red, the fertile matter in green, and the sodium plenum in yellow.

The one dimensional representative assembly of the ASTRID reactor used for this benchmark is described in the second section. The third section is devoted to the presentation of the calculation codes used in the study. The two neutronics approaches based on Monte Carlo and deterministic calculations are detailed, together with their spatial and point kinetics applications and their coupling to thermalhydraulics. Finally, the results obtained on a simplified ULOF transient are presented in the fourth Proceedings of ICAPP 2017 April 24-28, 2017 -Fukui and Kyoto (Japan) section of the paper together with a comparison to static Monte Carlo calculations on different snapshots of the fuel temperature and sodium density distribution taken at representative time steps of the transient.

II. CASE PRESENTATION

The application case is a 1D assembly representative of the ASTRID sodium cooled reactor at beginning of life, characterized by a negative sodium void effect. The radial distribution is neglected in this 1D study, assuming that there is no sodium boiling in any assembly. Such a reactor is very sensitive to power redistributions with the two fissile areas separated by a fertile area. A sodium plenum at the top of the geometry amplifies this effect by increasing the neutron leakages from the superior fissile area when the sodium density is reduced.

II.A. Case geometry

Figure 2 presents the case geometry. This configuration is very heterogeneous with fertile areas between fissile areas. The sodium plenum is optimized to ensure a negative sodium void effect (Ref. 2). If the sodium density decreases, the neutron absorption in the B 4 C increases so that this negative feedback leads to a reduction of the power.

II.B. Material characteristics

The material temperatures and the isotopic reference compositions of the different areas are given in Tab. I. As mentioned, the calculations presented here are performed at beginning of life. In order to compute the feedback coefficients, a material temperature increase of +300 K (Doppler effect) and a sodium density reduction of 1% are considered.

The nuclear database used in this study is JEFF 3.1 (Ref. 3).

II.C. Thermalhydraulics characteristics

The sodium is injected in the assembly at nominal condition at a flow rate of 31 kg/s, and at 400°C corresponding to a density of 0.85514 g/cm 3 . The sodium heat capacity is 1265.4 J/kg/°C and its expansion factor is 0.000315 K -1 .

The fuel thermal conductivity and heat capacity depend on the temperature and the burnup. The correlations used are accessible in Ref. [START_REF] Inoue | Fuel-tocladding gap evolution and its impact on thermal performance of high burnup fast reactor type uranium-plutonium oxide fuel pins[END_REF] and represented in Figs. 3 and4 at the beginning of the irradiation. Note the important variation of this parameter in the fuel temperature range.

III. CALCULATION CODE PRESENTATION

Two calculation codes are used in this study. Both are composed of a neutronics -thermalhydraulics coupling. This study aims to compare the neutronic modeling used: the deterministic S N (section III.A.2) and the TFM Monte Carlo (section III.B.2) models. Both are Proceedings of ICAPP 2017 April 24-28, 2017 -Fukui and Kyoto (Japan) capable to perform spatial kinetics and point kinetics calculations.

III.A1. Modeling used at IAEC

III.A.1. Thermalhydraulics

The thermalhydraulic part of the code describes the generation and removal of heat from an average fuel pin inside the core. The heat source, associated to the fission rate, is provided using the axial power density calculated by the neutronic part of the code throughout the transient. The assembly thermalhydraulics is modeled as a 1D flow with an imposed flow rate at the assembly inlet. At each time step, the local velocity is calculated using the imposed flow rate boundary condition and the sodium density variation. The heat removal from the cooling sodium channel is described by axial convection, and the heat removal from the fuel pin to the coolant by radial diffusion along the fuel pin and clad, assuming that the axial diffusion is negligible. The time dependent flow rate is an external parameter assumed to be known and thus not calculated by the code.

As mentioned, fuel and coolant temperatures provided by the thermalhydraulic part of the code are used to calculate either the space dependent cross-sections if one is interested in solving the spatial neutron transport equation, or to calculate the reactivity if one is interested to the describe the neutron evolution under the point kinetics (PK) assumption.

III.A.2. Neutronics

The deterministic calculations are based on the traditional two level lattice/core scheme. First, the selfshielded cross sections are computed at CEA by the ECCO (Ref 5) code cell, using the fundamental mode assumption for each kind of materials of the 1D core description. For fissile material, a buckling search algorithm is used to obtain the critical flux for the cross section collapsing to a 33 energy group mesh. For the subcritical materials such as fertile or structural parts of the 1D subassembly, the process is based on source calculations using the spectrum coming from previous fissile calculations. The isotopic cross sections are computed once for all at the different temperatures and considered as constant over time in the following calculations. Their dependency with the temperature is assumed to be linear with the sodium density variation and logarithmic with the temperature variation. The fundamental mode hypothesis used here does not seem to be the most realistic assumption regarding the 1D model of the core, but these calculation schemes inherited from PHENIX and SUPERPHENIX studies are the only ones available in ECCO package.

The neutronic part of the code then solves the time dependent transport equations, written for a one dimensional slab geometry (Ref. 6):

1 ݒ ߲߶ሺ,ݔ ,ܧ ߗ, ݐሻ ݐ߲ = ,ݔሺܮ ,ܧ ߗ, ݐሻ߶ሺ,ݔ ,ܧ ߗ, ݐሻ + 1 4ߨ ߯ ሺܧሻߣ ܥ ሺ,ݔ ݐሻ + ܵሺ,ݔ ,ܧ ߗ, ݐሻ, ߯ 4ߨ ܥ߲ ሺ,ݔ ݐሻ ݐ߲ = -߯ 4ߨ ߣ ܥ ሺ,ݔ ݐሻ + ܨ ௗ ሺ,ݔ ,ܧ ߗ, ݐሻ߶ሺ,ݔ ,ܧ ߗ, ݐሻ for ݅ = 1, … , ,ܫ (1)
where ܫ is the number of delayed neutron precursor families. The symbol ,ݔሺܮ ,ܧ ߗ, ݐሻ denotes the timedependent transport operator including streaming, absorption, scattering and prompt fission, and ܨ ௗ ሺ,ݔ ,ܧ ߗ, ݐሻ߶ሺ,ݔ ,ܧ ߗ, ݐሻ is the delayed neutron precursor group production rate.

The code solves the steady state equations for the neutron flux and its adjoint version for the adjoint flux by the ܵ ே method of discretization for the neutron angle ߗ (Ref. 7) with N = 16 and with 33 energy groups corresponding to the self-shielded cross sections prepared as described above. The boundary conditions are a flux leakage on the axial boundaries, and a volumetric leakage modeling the radial leakages.

The time evolution of the neutron flux is described following the assumption of the Predictor-Corrector version of the Improved Quasi-Static (IQS) method (Refs. 7 and 8): the time-dependent neutron flux is split into the product of two functions called amplitude and shape:

߶ሺ,ݔ ,ܧ ߗ, ݐሻ = ܲሺݐሻ߮ሺ,ݔ ,ܧ ߗ, ݐሻ (2)
The amplitude function ܲሺݐሻ depends only on time and provides the bulk information about the power change, while the shape function ߮ሺ,ݔ ,ܧ ߗ, ݐሻ depends on all the variables describing the (time-dependent) power profile deformation. The Predictor-Corrector version of the IQS method evaluates a predicted angular flux using a macro time scale for discretization, typically the same as the one chosen for the thermalhydraulic part of the code. The method assumes that the error resulting from the macro time discretization is related to the amplitude function alone. The predicted flux is used to estimate the shape function and to the calculated effective parameters for the point kinetic equations. Then one solves the point kinetic equations on a micro time scale, for a corrected estimated amplitude function. Once new amplitude is known, one gets a corrected value for the angular flux and the precursor concentrations. During the transient, local cross sections are given using a pre-calculated polynomial dependence on the coolant density and fuel temperature, Proceedings of ICAPP 2017 April 24-28, 2017 -Fukui and Kyoto (Japan) provided by thermalhydraulic part of the code. For the purpose of comparison, the point kinetics version of the solution is trivially achieved by assuming that the shape function remains constant throughout the transient. For this case, the code calculates at the beginning of the transient the local reactivity coefficients, which drives the point kinetic equations and then calculates at each time step the reactivity using the local values for the coolant density and fuel temperature provided by the thermalhydraulic part of the code.

III.B. Modeling used at CEA

III.B.1. Thermalhydraulics

The simplified thermalhydraulics calculations are performed with a 1D motion using the fluid mechanics calculation code OpenFOAM (Ref. 9). Despite reference tools for thermalhydraulics applied to fuel bundle exist, we choose a fast calculation tool developed in previous work (Ref . 10). The velocity field is represented in Fig. 5.

The section of the geometry simulated in OpenFOAM represents the average sodium axial section that changes in the different regions. Since the fluid dilatation is not considered in this thermalhydraulics computation, the fluid velocity only depends on the sodium fraction with a constant flow rate. Then, if the sodium section is smaller, the fluid velocity is larger as we can see in Fig. 5 with the smaller sodium section in the pin area (gas, fertile and fissile areas). Both the sodium and the fuel average temperatures are calculated during the transient. The heat exchange coefficient between the two media is computed on the fly at each axial position and time step. Using the previous average temperature estimation, the pin temperature radial distribution is reconstructed together with the sodium-fuel heat exchange. The pin temperature field is presented in Fig. 6. The difference on the mean temperature between the reference calculation (10000 bins in the fuel and 1000 bins in the steel clad) and the calculation used in this study (only 16 bins in the fuel and 4 in the steel clad) is limited to 0.3 K which is well negligible compared to the temperature variation in the pin. After a sensitivity study, the CEA calculations performed in the following use 16 bins in the fuel and 4 in the steel clad

III.B.2. Neutronics

The neutronics approach is based on a perturbative version of the Transient Fission Matrix (TFM) approach. The TFM approach described in Refs 10, 11 and 12 is based on the utilization of fission matrices and average time transport matrices to compute neutron kinetics. The matrices are estimated using a Monte Carlo calculation once per core configuration, and prior to the transient calculation. Different matrices ܩ , depending on the neutron spectrum (prompt ߯p or delayed ߯d) and the neutron multiplicity (prompt ߥp or delayed ߥd), are estimated during the same calculation. Those matrices contain the neutron propagation spatial and temporal behavior of the system, so that no more Monte Carlo calculation are required during the transient calculation thanks to an interpolation of the matrices performed on the fly to follow the system evolution. An interpolation model is thus implemented in the TFM approach.

Due to the large influence of the crossed volume during the neutron transport between its creation position and its fission position in an heterogeneous reactor, a correlated sampling (CS) technics (Ref. [START_REF] Bernnat | A monte carlo technique for local perturbations in multiplying systems[END_REF]) has been applied to study sodium cooled reactors such as ASTRID.

The CS technics provides a perturbed estimation of the fission matrices using a neutron weight perturbation that depends on the effect of a crossed section modification due to a sodium density or fuel temperature (Doppler) modification. Considering ܩ ఞ ೣ ఔ ೣ , any of the matrices of the TFM approach, and applying a CS weighting for a sodium density variation of -1% or a local temperature of +300 K, the effect of a local perturbation in the subvolume ݇ on ܩ is written ܩ ෨ ఞ ೣ ఔ ೣ ௗ for the density effect and ܩ ෨ ఞ ೣ ఔ ೣ ௗ for the Doppler effect. These matrices are called perturbed matrices. Finally, using ܶሺ݇ሻ and ∆ߩ ே ሺ݇ሻ the fuel temperature and the sodium density variation in volume ݇, the matrix ܩ is interpolated using a Proceedings of ICAPP 2017 April 24-28, 2017 -Fukui and Kyoto (Japan) linear interpolation for the density effect and a logarithmic interpolation for the Doppler effect as:

ܩ ఞ ೣ ఔ ೣ ሺ߂ߩ Na ሺ݇ሻ, ܶሺ݇ሻሻ = ܩ ఞ ೣ ఔ ೣ -∑ ܩ ෨ ఞ ೣ ఔ ೣ ௗ .߂ߩ Na ሺ݇ሻ + ∑ ܩ ෨ ఞ ೣ ఔ ೣ ௗ . logሺ்ሺሻ ் ref ⁄ ሺሻሻ logሺሺ் ref ሺሻା300ሻ ் ref ⁄ ሺሻሻ
(3) During transient calculations, the matrices of the TFM approach are then updated using this interpolation model in order to estimate the fission neutron redistribution and the reactivity variation calculated with the temporal integration presented in Ref. 11. The final kinetics equations solved for the prompt neutrons (N p ) and precursors of delayed neutrons of each family i (P i ) are the following: the fraction of delayed neutrons of family i. The matrix-vector multiplications correspond to the source terms, each source corresponding to a specific matrix. These matrices are updated at each time step using the interpolation model with the perturbed matrices and the sodium density and fuel temperature modification.

݀ ݐ݀ = ܩ ఞ ఔ 1 ݈ + ܩ ఞ ఔ ߣ - 1 ݈ ݀ ݐ݀ = ߚ ߚ ൭ܩ ఞ ఔ 1 ݈ + ܩ ఞ ఔ ߣ ൱ -ߣ (
As mentioned, the TFM approach requires to compute these specific matrices once, prior to the transient calculation. All the matrices are calculated with a modified version of the Serpent 2.1.21 code (Ref. [START_REF] Leppänen | The Serpent Monte Carlo code: Status, development and applications in 2013[END_REF]. One billion of neutrons are simulated and the system boundary is a neutron leakage.

Note that the reactor radial dimension has been adjusted to reach criticality, the value obtained with Serpent is k eff = 0.99980 ± 0.00002.

Concerning the point kinetic resolution, the power shape is constant and provided by the spatial kinetic calculation at nominal power. Instead of interpolating the fission matrices during the transient, the reactivity weight of each perturbed matrix is performed prior to the point kinetic calculation to produce the spatial feedback coefficient distribution in the reactor used during the transient. The distributions obtained for the ASTRID representative assembly considered here are presented in Fig. 7 and8.

We can see a good global agreement between the CEA-TFM (red curves) and the IAEC-S N (blue curves) approaches. A difference can be noticed in the gas plenum between 200 and 210 cm on the density feedback. This difference is due to the deterministic modeling approximation: this area is thin and the neutron spectrum is very different between the neutrons going from the fuel to the sodium plenum and the neutrons that are reflected from the sodium plenum. This effect is visible on the global feedback effect presented in Tab. II, where the reference value corresponds to a direct Monte Carlo calculation. The point kinetic values corresponds to the sum of the local contributions of Fig. 7 and8. The last column of Tab.II corresponds to a local sodium density decrease of -10 % in the upper gas plenum. Note that the value of PK-CEA corresponds to a sum of local contributions. There is a bias since the crossed effects between these local perturbations are not taken into account. If we generate the perturbed matrix associated to a global perturbation on the whole reactor, the reactivity variations are respectively of -20.2 pcm and -171 pcm for the density and the Doppler effects. These values are consistent with the direct Monte Carlo calculation (see Tab. II). The difference between CEA (Monte Carlo with local correlated sampling approach) and IAEC (S N perturbations using the direct and adjoint flux) comes from the neutron spectrum variation and anisotropy at the interface between the different areas. We can observe a larger effect on the density feedback in the Proceedings of ICAPP 2017 April 24-28, 2017 -Fukui and Kyoto (Japan) gas plenum in Fig. 8, also quantified in the right column of Tab. II. It may be explained either by the difference in the self-shielded cross sections between the reflected neutrons and the leaking ones due to the spectrum difference and the large steel density variation; or by the neutron leakage model used during the cross section selfshielding process.

IV. ULOF transient calculation

IV.A. Introduction

The Unprotected Loss Of Flow accident consists in a sodium flow rate reduction due to a pump failure. We consider as benchmark an imposed exponential reduction of the flow rate down to the minimal flow rate due to the natural convection of 7 % of the nominal value, with a decrease rate of ܶభ మ = 20 ݏ (see Fig. 9). Fig. 9. Flow rate evolution imposed during the ULOF transient.

As detailed in section II, two calculation codes from CEA and IAEC are used in this benchmark and both of them can use spatial kinetics (SK) and point kinetics (PK). They will be referred to as CEA-SK, CEA-PK, IAEC-SK and IAEC-PK. Just a sentence to remind that the CEA neutronics module uses the TFM approach based on a Monte Carlo pre-calculation, and the IAEC neutronics module uses a S N deterministic approach. One additional calculation called CEA-PK IAEC has been performed, using the IAEC PK parameters in the CEA numerical tool to check that the results obtained are the same with equivalent neutron kinetic inputs.

IV.B. Transient analysis

This section presents the results of the ULOF calculation using the CEA-SK and CEA-PK calculation codes. The SK approach provides the power redistribution, and the PK approach the local feedbacks. Figure 10 shows the assembly at the beginning (top) and at the end (bottom) of the transient calculated with CEA-SK. Due to the loss of flow, the sodium temperature is increasing during the transient, due to the negative feedback effect, the power level decreases. According to the power decrease, the fuel temperature also decreases to adjust the temperature difference between the fuel and the sodium. Finally, the power level stabilizes once the positive feedback of the fuel compensates the negative feedback of the sodium. The final power level here is around 20 % of the nominal power. This value results of this balance effect between the fuel and the sodium, so that it is sensitive to the feedback coefficients. The power redistribution in the assembly is around -5 % at the fuel -sodium plenum interface and +2 % at the bottom. Due to the sodium dilatation, the neutron leakage from the fuel to the B 4 C increases, resulting in a strong local reduction of the power shape. This effect results in a redistribution of the power in the lower fissile area. Fig. 10. Overview of the reactor at the beginning (top) and at the end (bottom) of the transient with a representation from left to right of: the power, the solid (fuel and structures) temperature, the heat transfer coefficient, the sodium temperature and the power redistribution.

Figure 11 presents the evolution of the axial power in the assembly (ordinate) as a function of time (abscissa) during the transient (left), together with the power redistribution in the assembly (right) calculated with CEA-SK. The redistribution is the normalized power variation compared to the initial value. Note that the horizontal dashed lines separate the different areas of the geometry (fuel, fertile matter, gas plenum…). The initial and final values correspond to the results presented in Fig. 10. During the transient the power level evolves, and we can see that the power redistribution is maximum at Proceedings of ICAPP 2017 April 24-28, 2017 -Fukui and Kyoto (Japan) around 100 s. As expected, the power level in the fertile area is low because of the small amount of fissile matter. Figures 12 and13 present the evolution of respectively the fuel temperature and the sodium temperature (left), together with their impact on the reactivity (right) due to the sodium dilatation and the Doppler effect using CEA-PK. The power production areas directly correspond to the positions where the fuel temperature is higher, and where the sodium temperature increases. Note that the sodium temperature stays below the boiling temperature (1155 K). We can see on the right the effect on the reactivity of the temperatures variations. The power reduction is due to the strong negative sodium density feedback effects in the sodium plenum. The contribution of the fuel area is slightly positive, but thanks to the progressive increase of temperature this component is smaller than the negative one. The fuel temperature reduction, induces a positive reactivity increase with a large component in the middle and the bottom of the upper fissile area. We can note a skin effect at the interface between the fuel and the gas plenum. Due to the sodium temperature increase, even if the power level is smaller, the middle fertile area temperature slightly increases, inducing a small positive reactivity component.

IV.C. Influence of the neutronics model

Proceedings of ICAPP 2017

April 24-28, 2017 -Fukui and Kyoto (Japan)

The system evolution corresponds to the description in section IV.B. Due to the feedback effect, the power level decreases during the transient, resulting in a reduction of the temperature difference between the sodium and the fuel. The reactivity reaches a minimum value during the transient, around -70 pcm, and slowly comes back to 0 pcm with a time constant associated to the precursor equilibration. All the calculation codes and neutronic models provide a similar behavior during the transient.

Prior to the precise comparison of the different neutronics models, one has to notice (Fig. 14, 15, 16) that the ULOF CEA-PK IAEC (green dashed line) calculation with the IAEC point kinetic parameters shows a very good agreement with the IAEC-PK calculation. This confirms that the differences in the thermalhydraulics modeling do not impact the results and the differences described in the next section are due to the neutronics models.

IV.C.1. Spatial kinetics versus point kinetics models

The PK and SK modeling shows a very good agreement respectively between CEA-PK and CEA-SK, and between IAEC-PK and IAEC-SK results. The power stabilization is the same and the temperature variations are superimposed. The difference on the reactivity maximum variation is limited to 3 pcm. Using point kinetics calculations is correct in such a slow transient, even for a heterogeneous geometry as the case considered here.

IV.C.2. TFM versus S N approaches

The larger difference is obtained between the TFM (CEA-PK/SK) and the S N (IAEC-PK/SK) neutronics models. The reactivity difference during the transient increases to a maximum value of 10 pcm. Due to the difference on the feedback effects, the final temperature equilibration is not the same. The final power level results of this coupling between the sodium density decrease with a negative reactivity contribution and the fuel temperature with a positive contribution. Finally, the propagation of the feedback difference between the two approaches implies a power level difference of 5% (18% for IAEC and 23% for CEA).

IV.D. Comparison with direct Monte Carlo calculations

In order to obtain the reference reactivity and power shape redistribution in the core during the transient, comparisons have been performed at specific time steps of the transient. From the coupled SK tools, the sodium density and the fuel temperature are extracted at t = 0 (initial state), t = 20 s and t = 400 s, and direct Monte Carlo calculations have been performed on these configurations. Note that the initial distributions (t=0) are different between CEA and IAEC since the equilibrium power shape are not exactly the same as it will be discussed (Fig. 17).

The reactivity and the power redistribution are calculated using the coupled SK results. In addition, Monte Carlo calculations are also performed using only density or only Doppler effect in order to compare the individual contributions. These individual reactivities are compared to the PK results to decompose the different contributions.

Tables III and IV respectively present the different reactivity variations for the IAEC and CEA results. The two middle columns are the decomposed reactivity variations, and the column on the right is the total variation. Note that since the Doppler and the sodium density effects have opposite contributions with the same order of magnitude, the total is small and then the relative difference may be very large between the reference and the calculated one. We can see that the larger differences on the individual contributions concern the density effect. Concerning the IAEC results it confirms the effect of the upper gas plenum in Fig. 8: the reactivity variation is overestimated in this area. A difference of 7±3 % at 20 s and 12±1 % at 400 s is also observed for TFM, the latter Proceedings of ICAPP 2017 April 24-28, 2017 -Fukui and Kyoto (Japan) being significant compared to the statistical error, and possibly illustrating the limit of the linear interpolation for the matrices interpolation model. Note that at t=0 s, the density variation between the bottom and the top of the geometry is around -4.3 %, while at 20 s it is around -7.1 % up to -14 % at 400 s.

Figs. 17 and 18 respectively present the flux distributions and redistributions in the core. The different time steps are presented in red (0 s), green (20 s) and blue (400 s). The reference Monte Carlo calculation is in plain line, and the results from IAEC and CEA are in dashed line. Even if a small difference due to the choice of a deterministic or Monte Carlo approach can be observed on the power distribution (Fig. 17), we can see on their variation (Fig. 18) that the three resolutions (S N , TFM and MC) are in very good agreement. The small difference that appears at 400 s is linked to the limit of the linear interpolation on the density effect. Finally, the capability of the IPK model to predict the flux redistribution without new spatial kinetic calculation (only sum of individual contributions) shows that the decoupling between the upper and lower fissile areas is very limited.

V. CONCLUSIONS

This study presents the impact of the neutron modeling on a representative ASTRID assembly during an ULOF calculation. Different neutronics models are tested. The first one, the TFM approach, is based on a Monte Carlo precalculation of the system response associated to an interpolation model, and the second one is a deterministic discrete ordinate S N approach solving the neutron equation in time, position, energy, and angle. These two approaches are used to perform both spatial and point kinetics calculations through the generation of local feedback coefficients.

In order to check the agreement of the thermalhydraulics and coupling tools, the ULOF transient has been calculated using the same point kinetics parameters with the two calculation codes and a very good agreement is obtained. The axial power shape redistribution during the considered transient is limited to +2/-4 %, the sodium density decrease at the top of the reactor inducing a power shift in the lower fissile area. However, the global system evolution is very similar for each code between its spatial kinetic and its point kinetic resolutions. The difference on the reactivity variation is limited to 3 pcm. The main discrepancy observed concerns the neutronics approaches TFM and S N used. The difference on the reactivity is around 10 pcm and, due to the discrepancy of the global feedback coefficients, the difference between the final power stabilization levels is of around 5%. Comparisons with static Monte Carlo calculations on different snapshots of the transient show that both approaches have a very good capability to evaluate the power redistribution in the core.

Proceedings of ICAPP 2017

April 24-28, 2017 -Fukui and Kyoto (Japan)

The study of this ULOF scenario highlights that, even if the spatial kinetics provides an information on the power redistribution, the neutronics modeling approach has a larger impact on the transient evolution. If the sodium remains liquid, point kinetics and spatial kinetics are in good agreement. Complementary studies will be required on other scenarios that may induce a spatial decoupling, with a larger perturbation amplitude such as sodium boiling or full core calculations with radially non uniform perturbations.

Fig. 2 .

 2 Fig. 2. Case geometry.

Fig. 3 .

 3 Fig. 3. Fuel thermal conductivity as a function of the fuel temperature.

Fig. 4 .

 4 Fig. 4. Fuel heat capacity as a function of the fuel temperature.

Fig. 5 .

 5 Fig. 5. Velocity field.

Fig. 6 .

 6 Fig.6. Radial temperature distribution associated to an imposed pin-power of 50 kW/m and a sodium temperature of 700 K. For the study presented in this paper, the calculation is performed with 16 bins in the fuel and 4 in the steel for the CEA calculation, and respectively 20 and 5 for the IAEC.

  [START_REF] Inoue | Fuel-tocladding gap evolution and its impact on thermal performance of high burnup fast reactor type uranium-plutonium oxide fuel pins[END_REF] where ݈ is the effective prompt lifetime calculated with the time matrix, andఉ ఉ బ

Fig. 7 .

 7 Fig. 7. Axial Doppler feedback coefficient.

Fig. 8 .

 8 Fig. 8. Axial density feedback coefficient.

Fig. 11 .

 11 Fig. 11. Power evolution map.

Fig. 12 .

 12 Fig. 12. Fuel temperature evolution map (left) and impact on the reactivity (right).

Fig. 13 .

 13 Fig. 13. Sodium temperature evolution map (left) and impact on the reactivity (right).

Figures 14 ,Fig. 14 .

 1414 Figures 14, 15 and 16 present the temperature, reactivity and power evolutions during the transient using the different codes and SK/PK models.

Fig. 15 .

 15 Fig. 15. Reactivity evolution.

Fig. 17 .

 17 Fig. 17. Power distribution at t=0 s (red), 20 s (green) and 400 s (blue) with IAEC (top) and CEA (bottom) results, reference in plain line and results from the coupling in dashed line.

Fig. 18 .

 18 Fig. 18. Power relative redistribution at t=0 s (red), 20 s (green) and 400 s (blue) with IAEC (top) and CEA (bottom) results, reference in plain line and results from the coupling in dashed line.

Transient coupled neutronics-thermalhydraulics study of ULOF accidents in sodium fast reactors using spatial kinetics: comparison of the TFM Monte Carlo and S N approaches A

  . Laureau 1 , Y. Lederer 2 , A. Krakovich 3 , L. Buiron 1 , B. Fontaine 1

	1: CEA, DEN, DER, Cadarache 13 108 Saint-Paul Les Durance Cedex, France
	axel.laureau@cea.fr
	2: Dept. of Physics, Nuclear Research Center-Negev, Beer-Sheva, Israel
	3: Dept. of Nuclear Physics & Engineering, Soreq Nuclear Research Center, Yavne, Israel

TABLE I .

 I Material temperature and composition -10 24 atoms per cm3

	Fert -1000 K	Fiss -1500 K	B 4 C -740 K
	16 O 1.952e-02
	23 Na 6.352e-03
	56 Fe 1.861e-02
	235 U 1.977e-05
	238 U 9.742e-03

TABLE II .

 II Global reactivity variation in pcm

	Case	Density	Doppler	-10% density
	-1% ߂ߩ ெିோ -21.2 ±1.6 ߂ߩ ିா -20.1 ߂ߩ ିூா -30.7	+300K -171.8±1.6 -182 -159	in 200-210 cm -51.2±2.4 -49.1 -111

TABLE III

 III

		. Reactivity variation -IAEC
	Case:	Doppler	Density	Doppler and
	20 s ߂ߩ ெିோ ߂ߩ ିூா difference Case:	47±2 pcm 48 pcm 2±4 % Doppler	-66±2 pcm -87 pcm 32±4 % Density	density -19±2 pcm -39 124±25 % Doppler and
	400 s ߂ߩ ெିோ 268±2 pcm -221±2 pcm ߂ߩ ିூா 266 pcm -269 pcm difference -1±1 % 21±1 %	density 48±2 pcm -3 pcm -106±8 %
	TABLE IV. Reactivity variation -CEA
	Case:	Doppler	Density	Doppler and
	20.3 s ߂ߩ ெିோ ߂ߩ ି்ிெ difference Case:	39±2 pcm 41 pcm 5±5 % Doppler	-77±2 pcm -71 pcm -7±3 % Density	density -42±2 pcm -32 pcm -24±6 % Doppler and
	400 s ߂ߩ ெିோ 231±2 pcm -278±2 pcm ߂ߩ ି்ிெ 246 pcm -246 pcm difference 6±1 % -12±1 %	density 45±2 pcm 0 pcm -100±9 %

ACKNOWLEDGMENTS

The authors wish to thank the IN2P3 department of the CNRS (National Center for Scientific Research) for its support during the development of the TFM approach, the IAEC (Israel Atomic Energy Commission) for its support during this research and to Dr. Erez Gilad for his contribution to the code development on the IAEC side.