A Laureau

Y Lederer

A Krakovich

L Buiron

B Fontaine

CEA/IAEC benchmark on neutronic approaches for ULOF and UTOP transient calculations in SFR

Keywords: ASTRID, benchmark, Spatial & point kinetics, S N & TFM neutronics models, ULOF, UTOP 16

The main purpose of this benchmark paper is to study and compare point and spatial neutronic approaches used to calculate ULOF and UTOP transients in sodium cooled fast reactors. A second objective is to compare deterministic and Monte Carlo calculations with two different calculation codes. 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 ULOF and UTOP calculations, 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.

Introduction

Simulating 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 complexity of the studied system 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 geometric heterogeneities of core concepts such as the low void concept ASTRID (Ref. [START_REF] Fontaine | Sodium-cooled fast reactors: the ASTRID plant project[END_REF]). 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) in the frame of a benchmark. These studies aim to define the validity domain of point kinetics models for system code applications. The one dimensional representative assembly of the ASTRID reactor used for this benchmark is described in section 2. Section 3 is devoted to the presentation of the 30 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 for ULOF (Un-35 protected Loss Of Flow) transients with different pump time constants and for both beginning of life and end of cycle configurations are presented in section 4. The UTOP (Unprotected Top Of Power accident) transient is presented in section 5.

Case presentation

The application case is a 1D assembly representative of the ASTRID sodium cooled reactor at beginning of life. Such a reactor can be 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 leakage from the superior fissile area when the sodium density is reduced. The radial distribution is neglected in this first 1D study, assuming that there is no radial power redistribution.

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. [START_REF] Sciora | Low void effect core design applied on 2400 MWth SFR reactor[END_REF]). If the sodium density decreases, the neutron absorption in the B 4 C increases so that this negative feedback leads to a decrease of the power.

Material characteristics

The material temperatures and the homogenized isotopic reference compositions of the different areas are given in tab. 2.

The calculations presented here are performed at beginning of life (by default) and at end of cycle with an atomic burnup of 12.2 and 9.92 for the top and bottom fissile areas, and 3.84 and 0.95 middle and bottom fertile areas. 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. [START_REF] Koning | The JEFF-3.1 nuclear data library[END_REF]).

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 75 heat capacity is 1265.4 J/kg/K and its expansion factor is 0.000315 K -1 . The fuel-clad exchange is 5 kW/m 2 /K, and the clad-coolant exchange is 85 kW/m 2 /K. The fuel thermal conductivity k th correlations are extracted from [START_REF] Inoue | Fuel-to-cladding gap evolution and its impact on thermal performance of high burnup fast reactor type uraniumplutonium oxide fuel pins[END_REF] (except for factor F 4 from [START_REF] Carbajo | A review of the thermophysical properties of MOX and UO 2 fuels[END_REF]). Using x the fraction of 80 oxygen in the heavy nucleus, T the fuel temperature in Kelvin, β BU the burnup in at.% and p the porosity, the thermal conductivity is calculated using:

k th = F 1 F 2 F 3 F 4 k 0 th with : k 0 th = 1 0.06059 + 0.2754 √ 2 -x + 2.011 • 10 -4 • T + 4.715 • 10 9 T 2 exp -16361 T F 1 = ω arctan(ω -1 ) F 2 = 1 + 0.019β 3 -0.019β 1 1 + e -T-1200 100 F 3 = 1 - 0.2 1 + e T-900 80
F 4 = 1 -p 1 + 2p where ω = 1.09β -3.265 + 0.0643β -0.5 T 0.5
The fuel thermal capacity C p (in J/kg/K) correlation extracted from [START_REF] Carbajo | A review of the thermophysical properties of MOX and UO 2 fuels[END_REF] is calculated using:

C p = C 0 + 2C 1 τ -3C 2 τ 2 + 4. * C 3 τ 3 -5.C 4 τ 3 -C -2 τ -2

. Thermalhydraulics model

The thermalhydraulic part of the code describes the generation and heat removal from an average fuel pin inside the core. The heat source, associated to the fission 100 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 105 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 dif-110 fusion 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 cal-115 culate 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 describe the neutron evolution under the point kinetics (PK) assumption.

Neutronics model

The deterministic calculations are based on the traditional two level lattice/core scheme. First, the self-shielded cross sections are computed at CEA by the ECCO (Ref [START_REF] Rimpault | The ERANOS code and data system for fast reactor neutronic analyses[END_REF]) cell code, using the fundamental mode assumption for each 125 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 sections collapsing to a 33 energy group mesh. For the subcritical materials such as fertile or structural parts of the 1D subassembly, the process is 130 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 135 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 140 PHENIX and SUPERPHENIX studies are the only ones available in the ECCO package.

The neutronic part of the code then solves the time dependent transport equations, written for a one dimensional slab geometry (Ref. [START_REF] Lewis | Computational methods of neutron transport[END_REF]):

1 v ∂φ(x, E, Ω, t) ∂t = L (x, E, Ω, t)φ(x, E, Ω, t) + 1 4π f χ f (E)λ f C f (x, t) + S(x, E, Ω, t), χ f 4π ∂C f (x, t) ∂t = - χ f 4π λ f C f (x, t)+ F d f (x, E, Ω, t)φ(x, E, Ω, t), f = 1, . . . , F, (1)
where F is the number of delayed neutron precursor families.

The symbol L (x, E, Ω, t) denotes the time-dependent transport operator including streaming, absorption, scattering and prompt fission, and

F d f (x, E, Ω, t)φ(x, E, Ω, t)
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 S N method of discretization for the neutron angle Ω (Ref. [START_REF] Dulla | Accuracy of a predictorcorrector quasi-static method for space-time reactor dynamics[END_REF]) 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 (Ref. [START_REF] Dulla | Accuracy of a predictorcorrector quasi-static method for space-time reactor dynamics[END_REF][START_REF] Dulla | The quasi-static method revisited[END_REF]): the time-dependent neutron flux is split into the product of two functions called amplitude and shape:

φ(x, E, Ω, t) = P (t)ψ(x, E, Ω, t) (2)
The amplitude function P (t) depends only on time and provides the bulk information about the power change, while the shape function ψ(x, E, Ω, t) 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 The simplified thermalhydraulics calculations are performed with a 1D motion using the fluid mechanics calculation code OpenFOAM (Ref. [START_REF] Jasak | OpenFOAM: A C++ library for complex physics simulations[END_REF]). Despite reference tools for thermalhydraulics applied to fuel bundle exist, we choose a fast calculation tool developed in previous work 195 (Ref. [START_REF] Laureau | Développement de modèles neutroniques pour le couplage thermohydraulique du MSFR et le calcul de paramètres cinétiques effectifs[END_REF]). 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 200 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).

205

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. The pin temperature radial distribution is reconstructed together with the av-210 erage sodium-fuel heat exchange. The pin temperature field is presented in Fig. 6 for each code, together with a reference solution (green) based on a large number of volumes (10000 bins in the fuel and 1000 bins in the steel clad). A sensitivity study on the number of discretizations 215 (see Fig. 7) leads to a discretization of 16 bins in the fuel and 4 in the steel clad to limit the error to 0.3 K which is well negligible compared to the temperature variation in the pin.

Neutronics model 220

The neutronics approach is based on a perturbative version of the Transient Fission Matrix (TFM) approach. The TFM approach described in Ref. [START_REF] Laureau | Développement de modèles neutroniques pour le couplage thermohydraulique du MSFR et le calcul de paramètres cinétiques effectifs[END_REF][START_REF] Laureau | Coupled neutronics and thermal-hydraulics transient calculations based on a fission matrix approach: application to the Molten Salt Fast Reactor[END_REF][START_REF] Laureau | Transient Fission Matrix: Kinetic calculation and kinetic parameters β eff and Λ eff calculation[END_REF] Those matrices contain the neutron propagation spatial and temporal behavior of the system, so that no more Monte Carlo calculation is required during the transient calculation thanks to an interpolation of the matrices performed on the fly to follow the system evolution.

235 Such an interpolation model is thus implemented in the TFM approach. Due to the large influence of the crossed volumes 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 245 density or fuel temperature (Doppler) modification. Considering G χxνx , 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 k on G χxνx is written 250 G den k χxνx for the density effect and G dop k χxνx for the Doppler effect. These matrices are called perturbed matrices. Finally, using T (k) and ∆ρ sodium (k) the fuel temperature and the sodium density variations in volume k, the matrix G χxνx is interpolated using a linear interpolation for 255 the density effect and a logarithmic interpolation for the Doppler effect as:

G χxνx (∆ρ sodium (k), T(k)) = G χxνx - k G den k χxνx • ∆ρ sodium (k) + k G dop k χxνx log( T (k) /T ref (k)) log( (T ref (k)+300) /T ref (k)) (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 260 the reactivity variation calculated with the temporal integration presented in Ref. [START_REF] Laureau | Coupled neutronics and thermal-hydraulics transient calculations based on a fission matrix approach: application to the Molten Salt Fast Reactor[END_REF]. The final kinetics equations solved for the prompt neutrons (N p (t)) and precursors of delayed neutrons of each family f (

f λ f P f (t)) are the fol- lowing: 265 dN p dt = G χpνp 1 l eff N p + G χ d νp f λ f P f - 1 l eff N p dP f dt = β f β 0 G χpν d 1 l eff N p + G χ d ν d f λ f P f -λ f P f (4)
where l eff is the effective prompt lifetime calculated with the time matrix, and β f/β 0 the fraction of delayed neutrons of family f. The matrix-vector multiplications correspond to the source terms, each source corresponding to a specific matrix. These matrices are updated at 270 each time step using the interpolation model with the perturbed matrices and the sodium density and fuel temperature modifications.

As mentioned, the TFM approach requires to compute these specific matrices once, prior to the transient calcu-275 lation. All the matrices are calculated with a modified version of the Serpent 2.1.21 code (Ref. [START_REF] Leppänen | Kalti-640 aisenaho, 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 ad-280 justed to reach criticality, the value obtained with Serpent is k eff = 0.99980 ± 0.00002. Concerning the point kinetics resolution, the power shape is constant and provided by the spatial kinetic calculation at nominal power. Instead of interpolating the 285 fission matrices during the transient, the reactivity weight of each perturbed matrix is performed prior to the point kinetics calculation to produce the spatial feedback coefficient distribution in the reactor used during the transient. The distributions obtained for the ASTRID representative 290 assembly considered here are presented in Fig. 8 and9. 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. 3, where the reference value corresponds to a direct Monte Carlo calculation. The point kinetics values correspond to the sum of the local contributions of Fig. 8 and9 Note that the value of PK-CEA corresponds to a sum of local contributions. There is a bias since the crossed ef-fects between these local perturbations are not taken into account. If we generate the perturbed matrix associated 310 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. 3). The difference between CEA (Monte Carlo with 315 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 gas plenum in Fig. 9, 320 also quantified in the right column of Tab. 3. 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 325 during the cross section self-shielding process.

ULOF transient calculation

Introduction

The Unprotected Loss Of Flow accident consists in a sodium flow rate reduction due to a pump failure. We con-330 sider for the 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 T1 /2 = 20s (see Fig. 10). 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 340 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 345 with equivalent neutron kinetic inputs.

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 11 shows the assembly at the beginning (top) and at the end (bottom) of the transient calculated with CEA-SK. Due to the to right of: the power, the solid (fuel and structures) temperature, the heat transfer coefficient, the sodium temperature and the power redistribution loss of flow, the sodium temperature is increasing during the transient. Due to the negative feedback effect, the power level decreases. According to this 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. The areas with the higher power production directly correspond to the positions where the fuel temperature 390 is higher (in red at the beginning of the transient), and where the sodium temperature increases. Note that the sodium temperature stays below the boiling temperature (1155 K).

We can see (Fig. 13 and 14 -right) the effect on the 395 reactivity of the temperature 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 400 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. 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.

Influence of the neutronics model

Prior to the precise comparison of the different neutronics models, one has to notice (Fig. 15, 16, 17) that the ULOF CEA-PK IAEC (green dashed line) calculation with the IAEC point kinetics 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.

Spatial kinetics versus point kinetics models

The PK and SK modelings show a very good agreement 430 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 435 is correct in such a slow transient, even for a heterogeneous geometry as the case considered here.

TFM versus S N approaches

The larger difference is obtained between the TFM (CEA-PK/SK) and the S N (IAEC-PK/SK) neutronics 440 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 and the discrepancy is around 100 K. The final power level results of this coupling be-445 tween 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).

Influence of the flow rate reduction time constant

In order to confirm these conclusions with a smaller transient time constant, the ULOF calculation has been performed with a time constant of T1 /2 = 5 s instead of 20 s. The results on the temperature, the reactivity and the power ratio are represented in Figs. 18, 19 and 20. Similar results are obtained. Due to the smaller time constant, the minimum reactivity is smaller, around -180 pcm at 22 s versus -70 pcm at 70 s in the previous ULOF calculation. The minimum flow rate being unchanged, the power stabilises at the same final power. The main difference between the results are linked to the neutronics approach (deterministic of stochastic). Finally we can see that on such a fast ULOF simulation, no important difference still appears between spatial and point kinetics.

Influence of the burnup

The ULOF transient has also been studied for the end of cycle (EOC) configuration presented in section 2.2. The We can see that the final power level is higher than at beginning of life due to the difference in the feedback effects. A global reduction of the sodium density of 1% induces a reactivity variation of -2.9 pcm versus -20.1 pcm 480 at beginning of life. Indeed the contribution of the sodium plenum in the density feedback distribution (Fig. 22) is decreased together with the importance of the sodium reflected neutrons since a larger part of the fissions occurs in the fertile areas. The difference between IAEC and CEA is Due to the reactivity insertion starting at 10 s, the first impact on the system is an increase of the power that At 100 s, the power stabilises at a maximum level, the thermal feedbacks compensating exactly the 100 pcm in-515 serted. The power redistribution is limited to 0.5 % due to the small sodium temperature increase compared to the ULOF case, the redistribution being largely linked to the neutron absorptions in the B 4 C.

Influence of the neutronics model

The calculation has been performed using both CEA/I-AEC and SK/PK tools and models. Figures 15,16 and 17 present the results of each case for the temperature, the reactivity and the power variation. A similar behaviour between each code and method is 525 obtained on each parameter. The results between the SK and the PK are very close, and a small difference exists between the CEA and IAEC results. The slight difference on the fuel maximum temperature is due to the power shape that is slightly different and impacts this maximum value. However, it does not impact the relative power stabilisation that is around 120 % in both case. A small difference also appears on the reactivity before the beginning of the transient (10 s) due to the slow SK stabilization of the system, even if this reactivity variation 535 is limited to 0.2 pcm and is negligible. We can see that the maximum value reached during the transient is only around 7 pcm at 100 s for 100 pcm inserted thanks to the feedbacks, and the exact value depends on the feedback coefficients of each approach.

Conclusion

This article presents the impact of the neutron modeling on a representative ASTRID assembly for ULOF and UTOP calculations. Different neutronics models are tested. The first neutronic model considered, the TFM 545 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 per-550 form 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 555 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 % for the ULOF calculation, the sodium density decrease at the top of the reactor inducing a power shift in the lower fissile area. The global system evolution is very similar for each code between its spatial kinetic and its point kinetics resolutions. The difference on the reactivity variation is limited to a few pcm. The main discrepancy observed concerns the neutronics TFM and S N approaches used. The difference on the reactivity is limited to around 10 pcm for the ULOF and 1 pcm for the UTOP and, and due to the discrepancy of the global feedback coefficients, the difference between the final power stabilization levels is of around 5 % for the ULOF case.

The benchmark highlights that the neutronics modeling approach (S N vs TFM) has a larger impact on the transient evolution than the spatial kinetics compared to the point kinetics even the power redistribution provides a useful information. 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.

Figure 1 :

 1 Figure 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.

Figure 2 :

 2 Figure 2: Case geometry

Figure 3 :

 3 Figure 3: Fuel thermal conductivity as a function of the fuel temperature

Figure 4 :

 4 Figure 4: Fuel heat capacity as a function of the fuel temperature

95 3 . 1 .

 31 Modeling used at IAEC 3.1.1

  120

3 . 2 .

 32 function and to the calculated effective parameters for the point kinetics equations. Then one solves the point kinetics 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, 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 the local reactivity coefficients at the beginning of the transient, which drives the point kinetics equations and then calculates at each time 185 step the reactivity using the local values for the coolant density and fuel temperature provided by the thermalhydraulic part of the code. Modeling used at CEA and single physics comparison 3.2.1. Thermalhydraulics model190

Figure 5 :

 5 Figure 5: Velocity field

  is based on the utilization of fission matrices and average time transport matrices to compute neutron kinetics. The matrices 225 are estimated using a Monte Carlo calculation once per core configuration, and prior to the transient calculation. Different matrices G, depending on the neutron spectrum

Figure 6 :

 6 Figure 6: Radial temperature distribution associated to an imposed pin-power of 50 kW/m and a sodium temperature of 700 K, with 16 bins in the fuel and 4 in the steel for the CEA calculation, and respectively 20 and 5 for IAEC

Figure 7 :

 7 Figure 7: Error on the average fuel temperature as a function of the number of discretizations

Figure 9 :

 9 Figure 8: Axial Doppler feedback coefficient

Figure 10 :

 10 Figure 10: Flow rate evolution imposed during the ULOF transient As detailed in section 3, two calculation codes from

Figure 11 :

 11 Figure 11: 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 12

 12 Figure 12 presents the evolution of the axial power in the assembly (ordinate) as a function of time (abscissa)

Figure 12 :

 12 Figure 12: Power evolution mapThe redistribution is the normalized power variation

Figure 13 :

 13 Figure 13: Fuel temperature evolution map (left) and impact on the reactivity (right)

Figure 14 :

 14 Figure 14: Sodium temperature evolution map (left) and impact on the reactivity (right)

Figures 15 ,Figure 15 :

 1515 Figure 15: Evolution of the fuel and sodium maximum temperatures -ULOF

  Figure 16: Reactivity evolution -ULOF

Figure 18 :Figure 19 :

 1819 Figure 18: Evolution of the fuel and sodium maximum temperatures -fast ULOF

Figure 20 :Figure 21 :

 2021 Figure 20: Relative power evolution -fast ULOF

485Figure 22 :Figure 23 :Figure 25 :

 222325 Figure 22: Axial density feedback coefficient at end of cycle

500 5 . 2 .

 52 Transient analysisDuring the first 10 seconds, the power, the fuel temperature and the sodium temperature are stable in the core as respectively presented on the left column of Figs. 26, 27 and 28.

  505

Figure 26 :

 26 Figure 26: Power evolution map

Figure 28 :

 28 Figure 28: Sodium temperature evolution map (left) and impact on the reactivity (right)

Figure 29 :

 29 Figure 29: Evolution of the fuel and sodium maximum temperatures -UTOP

Figure 31 :

 31 Figure 30: Reactivity evolution -UTOP

  540

Table 1 :

 1 The pin radial description is given in table1. Pin radial geometry

	Material	Radium [mm]
	Central hole (He)	1.100
	Fuel/Fertile pellet	4.225
	Pellet gap (He)	4.350
	Clad	4.850
	Spacing wire	0.500

Table 2 :

 2 Material temperature and composition -10 24 atoms per cm3

Table 3 :

 3 . The last column of Tab. 3 corresponds to a local sodium density decrease of -10 % in the upper gas plenum. Global reactivity variation in pcm

	Case	Density	Doppler	-10% density
		-1%	+300 K	in 200-210 cm
	∆ρ Ref ∆ρ CEA ∆ρ IAEC	-21.2±1.6 -171.8±1.6 -20.1 -182 -30.7 -159	-51.2 ± 2.4 -49.1 -111

Preprint submitted to Annals of Nuclear EnergyJuly 26, 2017

Acknowledgments

The authors wish to thank the IN2P3 department of the CNRS (National Center for Scientific Research) for its support during the initial 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.


Pietro Mosca

Laura Clouvel

Théophile Lacaze

Sébastien Lahaye

J.-J Lautard

A.-M Baudron

Apollo3 Mary

Helen Sparks

; K Russell Depriest

Jean Jacques Lautard

Anne Marie Baudron

Deterministic model of PWR fast

Keywords: RPV, NEUTRON FLUENCE, UNCERTAINTY, APOLLO3 ®, URANIE, AEMC

come

INTRODUCTION

The fast neutron fluence and the corresponding uncertainties are important parameters for the reactor pressure vessel life time. Usually for this kind of studies, reference tools are based on Monte Carlo codes. However, the CPU time required is still not reasonable for uncertainty and sensitivity analyses. In this perspective, this article presents the model, under development at CEA (Commissariat à l'Energie Atomique et aux Energies Alternatives), to carry out with the deterministic code APOLLO3 ® [START_REF] Schneider | APOLLO3 ® : CEA/DEN deterministic multi-purpose code for reactor physics analysis[END_REF] uncertainty calculations of the fast fluence for PWR irradiation surveillance. The new numerical scheme is based on the 3D-SN solver MINARET [START_REF] Lautard | MINARET, a deterministic neutron transport solver for nuclear core calculations[END_REF] of the APOLLO3 ® code which uses the discontinuous Galerkin finite elements approximation. The spatial mesh is unstructured and the transport calculations are parallelized with respect to the angular directions. In this scheme, the multigroup cross-sections are sub-group self-shielded and collapsed over a dedicated energy mesh optimized by the Adaptive Energy Mesh Constructor (AEMC) [START_REF] Mosca | An Adaptive Energy Mesh Constructor for Multigroup Library Generation for Transport Codes[END_REF]. The structure of the paper is as follow. First the choice of the multigroup library and the optimization of the energy mesh to collapse cross-sections are described. Second the validation of the deterministic fast fluence scheme respecting the Monte Carlo reference TRIPOLI-4 ® [START_REF]TRIPOLI-4 version 8.1, "3D general purpose continuous energy Monte Carlo Transport code[END_REF] is shown on a three-dimensional PWR geometry. Then an example of uncertainty propagation coupling MINARET to the uncertainty and sensitivity platform URANIE [START_REF] Gaudier | URANIE: The CEA/DEN Uncertainty and Sensitivity platform[END_REF] is presented and results are discussed.

VALIDATION OF THE MULTIGROUP LIBRARY

As a first step for the development of the deterministic fast fluence scheme, the choice of the multigroup library for the production of the core collapsed cross-section library is done. The considered libraries are mainly based on the JEFF3.1.1 evaluation and are declined over optimized 300, 600, 1200 group energy meshes by AEMC for fast reactor systems [START_REF] Mosca | An Adaptive Energy Mesh Constructor for Multigroup Library Generation for Transport Codes[END_REF], [START_REF] Mosca | Improvement in Transport Calculations by the Energy-Dependent Fission Spectra and Subgroup Method for Mutual Shielding[END_REF]. A simplified slab model, which corresponds to a section cut of the reactor going from the core peripheral to the inner part of the vessel, is used to evaluate the accuracy of these different multigroup libraries. Vacuum boundary conditions are applied on the right end side and an external source distributed in energy according to the fissile spectra of the core is set on the left end side. The transport calculations are carried out with the APOLLO2 [START_REF] Sanchez | APOLLO2 YEARS 2010[END_REF] discrete ordinates solver with a S32 angular discretization and a P5 approximation for the scattering cross sections. Cross-sections of all the light nuclides of the steel are subgroup self-shielded and differentiated per media of the slab. Table 1 shows the accuracy of the different deterministic solutions with respect to the reference TRIPOLI-4 ® using as discriminant parameter the total flux over 1 MeV, noted Φ 1MeV in this work. One can notice that passing from 300 to 600 group library the accuracy is almost increased of a factor two while a further energy mesh refining does not much more improve the accuracy of the solution. Fig. 1 shows that the corresponding multigroup fluxes in the vessel are within TRIPOLI-4 ® error bars in the energy range between 1 and 20 MeV. These results validate the use of the 600-group library for the production of the core collapsed cross-section library.

OPTIMIZATION OF THE ENERGY MESH FOR CROSS-SECTION COLLAPSING

Since the memory required by the 600 group library is too huge to be used in a three-dimensional fluence core calculation, a cross-section collapsing is mandatory. To address such a problem, the approach described in [START_REF] Mosca | Energy mesh optimization for multi-level calculation schemes[END_REF] has been adapted to the fluence context. The optimization is accomplished on the slab geometry of the previous section, considering a pure uranium and plutonium fission source. The energy mesh bounds are adapted by AEMC [START_REF] Mosca | An Adaptive Energy Mesh Constructor for Multigroup Library Generation for Transport Codes[END_REF] in order to minimize the errors between two successive SN transports. The reference SN calculation is carried out over the fine 600-group energy mesh and the approximate one is obtained from collapsed cross-sections by flux homogenization over a coarser energy mesh.

To accomplish the optimization, it is defined the following functional that depends on the coarse mesh 𝑀 𝑐 of 𝑁𝐺 groups:

𝐹(𝑀 𝑐 ) = ∑ 𝛼 𝑛 [𝐹 𝑛 (𝑀 𝑐 )] 2 𝑁 𝑛=1 + ∑ [𝐺 𝑘 (𝑀 𝑐 )] 2 , 𝑁𝐺 𝑘=1 (1)
where 𝑁 = 𝑁𝐼 × 𝑁𝑀 × 𝑁𝑅 × 𝑁𝐺 × 𝑁𝐶 is the total number of components, 𝑁𝐼 is the total number of isotopes, 𝑁𝑀 is the total number of media, 𝑁𝑅 represents the total number of reactions (𝑁𝑅 = 2 corresponding to absorption and scattering) and 𝑁𝐶 is the total number of problems to be considered in the optimization (NC=2 for the two types of sources). The 𝐹 𝑛 (𝑀 𝑐 ) components represent the discrepancy between the reference isotopic reaction rate and its approximate value 𝑇 𝜌,𝑥,𝑖 𝑟,𝐺 , 𝑇 𝜌,𝑥,𝑖 𝑎,𝐺 renormalized to the corresponding total isotopic reaction rates

1 ≤ 𝜌 ≤ 𝑁𝑅 , 1 ≤ G ≤ 𝑁𝐺, 1 ≤ 𝑖 ≤ 𝑁𝐼, 1 ≤ 𝑥 ≤ 𝑁𝑀 × 𝑁𝐶. (2)
The use of both the weights 𝛼 𝑛 and the renormalization is justified in [START_REF] Mosca | Energy mesh optimization for multi-level calculation schemes[END_REF]. The 𝐺 𝑘 (𝑀 𝑐 ) components measure the discrepancy between the reference and the approximate leakage terms:

𝐺 𝑛 (𝑀𝑐) = [ 𝐽 𝑟,𝐺 (𝑀 𝑐 ) 𝐽 𝑟 (𝑀 𝑐 ) - 𝐽 𝑎,𝐺 (𝑀 𝑐 ) 𝐽 𝑎 (𝑀 𝑐 ) ]. (3)
The 600-group library has 96 groups above 1 MeV, the remaining groups beneath, which are useless with regard to the problem being treated, are collapsed in a unique group. Indeed, the flux responsible for damaging the vessel is currently defined as the total neutron flux greater than 1 MeV. AEMC has therefore determined a coarse mesh consisting of a single group between 10 -11 MeV and 1MeV and (NG -1) energy groups greater than 1MeV (with NG <96 the available number of groups). Fig. 2 shows the dependency of the functional 𝐹 𝑛 (𝑀𝑐) with respect to the number of groups. One can notice that the accuracy does not improve itself over 19 groups. The same tendency has been verified comparing Φ 1MeV , which are obtained with different optimized energy meshes, and the Monte Carlo reference. These results lead to the use of the 19-group optimized energy mesh in the core fluence calculations. This energy mesh contains 1 group between 10 -11 and 1 MeV and 18 groups between 1 and 20 MeV.

VALIDATION OF THE DETERMINISTIC FLUENCE MODEL

The validation of the deterministic fluence model is carried out modeling one eights of a 900 MWe PWR core. The reflector geometry is obtained by axial extrusion of five radial basic geometries that model the baffles, the surveillance capsules, the stiffeners, the cold and the hot water part at the bottom and the top of the core (Fig. 3). The core is composed of 26 assemblies and each assembly is a standard 17x17 UO2 homogeneous pin cells. The homogenization is justified by the fact that neutron mean free path in UO2 fuel and in borated hot water is higher than the size of the cells at energies beyond 1 MeV. Only the four last peripheral assemblies of the core are taken into account in the transport calculation (Fig. 3). This simplification reduces of one third the CPU time preserving the accuracy of the solution because 99% of the fluence in the vessel and in the surveillance capsules is due to neutrons from the peripheral assemblies. The power distribution in the core is an input datum which is directly converted into a neutron source radially homogenized on the fuel cells and discretized into 56 axial planes.

The axial temperature profiles of the water in the core and in the reflector correspond to standard operating conditions, elsewhere only average operating temperature values are considered. Since Doppler effect is negligible on light nuclei and also on heavy nuclei above the resonance energy range (energy higher than 300 keV), the water density is the only considered physical parameter to be modified by the temperature. All the cross-sections are calculated at the average operating temperatures. The MINARET calculations are carried out using a S32 angular discretization, a P3 approximation for the scattering cross sections, a radial mesh size of 2 cm, an axial mesh size 10 cm. The radial and the axial polynomial order are set equal to 1 and 0 respectively and the relative error per group on the flux is imposed lower than 10 -4 . These parameters are issued of a preliminary convergence analysis of the solution.

The reference calculation is carried out with the TRIPOLI-4 ® Monte Carlo code using a geometrical description without homogenization of the core. Multigroup and pointwise crosssections are issued from the processing of the JEFF3.1.1 evaluation. From Table 2 one can observe that relative error on the total flux over 1 MeV is lower than 1%. Fig. 4 shows that the multigroup flux discrepancy in the vessel is lower than 5% between 1 and 20 MeV. These results prove that the higher complexity of the 3D core test does not impact the accuracy of the solution because it is comparable to the slab test (cf. Fig. 1). The APOLLO3 ® transport calculation took 20 minutes and required a total memory of 113 Gb with a MPI (Message Passing Interface) parallel computing on a LINUX CenOS6.4 cluster of 4 processors with 20 cores each. The corresponding TRIPOLI-4 ® reference calculations required 6 hours with the same demand of processor resources. These results support the use of MINARET to perform uncertainty studies of neutron fluence using a Monte Carlo approach.

UNCERTAINTY PROPAGATION

The deterministic fluence scheme is used in order to lead a simplified uncertainty propagation to evaluate the impact of technological parameters (manufacturing tolerance of in-vessel components and water temperature). This study is carried out coupling the MINARET solver to the CEA uncertainty and sensitivity platform URANIE [START_REF] Gaudier | URANIE: The CEA/DEN Uncertainty and Sensitivity platform[END_REF] which pilots the design of experiment (DOE).

Table 3 shows the perturbed parameters and the corresponding variations in the study. In our case, we consider that the variation of the uncertain variables follows a uniform law centered on the nominal value; that is, the realization of the geometric values in the finite set, defined by the uncertainties, is equally probable. For temperature sampling, the perturbation is the same at each axial step of the core, since total power of the reactor is not an uncertain variable of the study. The Latin hypercube sampling method is used by URANIE for the sampling. The DOE is composed of 2000 samples for each uncertainty variables in the study. The choice of the uniform distribution law is justified by the fact that only the extreme values for each uncertain variable are provided and any other laws would give supplementary and arbitrary information. With URANIE a sensitivity analysis is also carried out, in order to evaluate if the variations of Φ 1MeV can be approached by a linear, additive or nonlinear form of the uncertainty variables, and also to determine the most influential ones in the point of interest. Fig. 5 illustrates the impact of the nine variables on Φ 1MeV in the surveillance capsules and the vessel. One can observe that the relation between Φ 1MeV and the uncertainty variables is well fitted by a linear relation because the coefficient of determination R 2 is higher than 0.99. The sensitivities are thus simply proportional to the coefficient of the linear regression. Fig. 6 classes the sensitivities of the uncertainty variables in the surveillance capsules. In the 17° capsule, the angular position, the water temperature, the internal radius of the shield and of the envelope have more important impact on Φ 1MeV than the other parameters because of much higher range of variation. A similar analysis shows that the vessel is very sensitive to its internal radius, and much less to the internal radii of the envelope and the shield. It is almost unaffected by other uncertain variables because of its distance from the core and the shielding effect produced by the interposed water. Fig. 7 shows the impact of the uncertain variables on the multigroup flux in the surveillance capsules. The angular position is anti-correlated with the multigroup flux because the capsule gets closer to the vessel point of highest flux ("hot spot") when the angle decreases. The most impacted energetic range is between 1 and 3 MeV because by approaching the peripheral assemblies the weight of the fission spectrum increases (energy range between 1 and 3 MeV) and the probability is higher that neutrons with lower energy reach the capsule. The shield internal radius and the multigroup flux are positively correlated in the range between 2 and 20 MeV, because the slowing-down power of the steel decreases due to the reduction of its corresponding thickness. On the other hand, the anti-correlation below 3 MeV is due to a high neutron transfer below 1 MeV. This effect is produced by the increase of the water thickness and by the high cross-section value of the water in this energy range (Fig. 8).

Figure 1

 1 Figure 1 Multigroup flux over 1 MeV in the vessel medium of the slab obtained running APOLLO2 with 300 and 600-group libraries and compared to TRIPOLI-4 ® reference. A2 and T4 stand for APOLLO2 and TRIPOLI-4 ® respectively.

Figure 2

 2 Figure 2 Funtional value of the energy mesh optimization versus the number of groups

Figure 3 :

 3 Figure 3: Two-dimensional geometry of the PWR core

Figure 4 :

 4 Figure 4 : Multigroup flux over 1 MeV in the vessel of 900 MWe PWR obtained with APOLLO3 ® and compared to the reference TRIPOLI-4 ® . A3 and T4 stand for APOLLO3 ® and TRIPOLI-4 ® respectively.

Figure 5 :Figure 6 :

 56 Figure 5 : DOE for the total flux over 1 MeV in surveillance capsules, the vessel cladding and the vessel

Figure 7 :

 7 Figure 7: Correlation between variables and multigroup flux in the surveillance capsule.

Table 1 : Total flux over 1 MeV in a slab obtained running APOLLO2 with different multigroup libraries and compared to TRIPOLI-4 ® reference. A2 and T4 stand for APOLLO2 and TRIPOLI-4 ® respectively. 1-Φ 1MeV A2 /Φ 1MeV T4 300 groups 600 groups 1200 groups 3σ (T4)

 1

	Capsule	-1.48%	-0.81%	-0.69%	0.38%
	Vessel	-1.25%	-0.69%	-0.44%	1.20%

Table 2 : Total flux over 1 MeV for the 900 MWe PWR obtained with APOLLO3 ® and compared to the reference TRIPOLI-4 ® . A3 and T4 stand for APOLLO3 ® and TRIPOLI-4 ® respectively.

 2

	Medium	1-Φ 1MeV A3 /Φ 1MeV T4	
	Vessel	0.61 %	0.24%
	Capsule 20°	0.86 %	0.09%
	Capsule 17°	0.30 %	0.09%

Table 3 : Uncertainty variables

 3

	Parameter	Perturbation
	Azimuth	α ∓ 0.33°
	Internal Vessel Radius	𝑅 𝑖𝑛𝑡𝑣𝑒𝑠 ±5.8𝑚𝑚
	Internal Envelope Radius	𝑅 𝑖𝑛𝑡𝑒𝑛𝑣 ±4𝑚𝑚
	External Envelope Radius	𝑅 𝑒𝑥𝑡𝐸𝑛𝑣 ±1𝑚𝑚
	Internal Shield Radius	𝑅 𝑖𝑛𝑡𝑠ℎ ±4.8𝑚𝑚
	Capsule thickness	𝐸 𝑐𝑎𝑝𝑡ℎ𝑖𝑐𝑘 ±0.38𝑚𝑚
	Water gap	𝐸 𝑤𝑎𝑡𝑔𝑎𝑝 ±0.5𝑚𝑚
	Capsule position	𝐸 𝑐𝑎𝑝𝑝𝑜𝑠 ±0.2𝑚𝑚
	Water temperature	𝑇 𝑤𝑎𝑡 ±1°

Capsule 17°Capsule 20°Vessel Cladding Vessel Temperature Azimut Vessel Radius Int Envelope Radius Ext Envelope Radius Int Shield Radius Capsule thickness Water gap

	Capsule position
	R²			0,99711			0,994609		0,999351		0,999507
	Mean		1,82E+11			1,56E+11		2,50E+10		9,73E+09
	Sigma		3,43E+09			1,84E+09		1,18E+09		3,44E+08
	Sigma/Mean (%)			1,89			1,18			4,71			3,54
	Nominal Value		1,82E+11			1,56E+11		2,50E+10		9,72E+09
	Ecart (%)			-0,13			-0,11			-0,15		-0,11
							Capsule 17°
	h	Azmut	Rintves	v	Rinten	Twat	Cappos	v	Rexten	Watgap	Capthick	s	Rintve

DEN/DANS/DM2S/SERMA/LLPR, Bâtiment 470, 91191 Gif sur Yvette Cedex, France, pietro.mosca@cea.fr

DEN/DANS/DM2S/SERMA/LLPR, Bâtiment 470, 91191 Gif sur Yvette Cedex, France

CONCLUSIONS

This paper shows the new deterministic fast fluence scheme with the APOLLO3 ® code developed at CEA. The use of sub-group self-shielded cross-section collapsed over an optimized collapsing energy mesh by AEMC allows to carry out very precise deterministic calculations of neutron flux in the vessel. Results of this deterministic model are encouraging, the integral of the flux over 1 MeV in the location of interest is calculated in less than 20 minutes with an error lower than 1% with respect to the TRIPOLI-4 ® Monte Carlo reference. An example of uncertainty propagation with the CEA uncertainty and sensitivity platform URANIE is also provided. This example offers an alternative approach to the classic first order sensitivity analysis to perform uncertainty studies like in [START_REF] Kodeli | Multidimensional Deterministic Nuclear Data Sensitivity and Uncertainty Code System : Method and Application[END_REF]. It will be pursued in order to estimate the impact of nuclear data (cross-sections, fission spectra) and of other technological parameters (power distribution, power intensity …). The authors acknowledge EDF (Electricité De France) for the financial and technical support.