Investigation of discretization uncertainty in Monte Carlo neutron transport simulations of the Molten Salt Fast Reactor (MSFR)

In the present work, an assessment of the Neutronic Benchmark of the Molten Salt Fast Reactor (MSFR) was performed using mesh based Monte Carlo Neutron Transport (MCNT) calculations with numerical uncertainty quantification due to discretization in neutronic parameters. Calculations with Constructive Solid Geometry (CSG) models where made as a baseline for the developed mesh based models. The numerical uncertainty given by the mesh utilization is evaluated using an extended version of the Grid Convergence Index (GCI). The fuel salt reprocessing is evaluated regarding a constant reprocessing rate. The fuel salt inventory variation with time for the developed models (CSG and meshed) is presented. The differences caused by the discretization procedure are noticeable, which shows that mesh based MCNT require careful mesh sensitivity evaluation and further validation.


INTRODUCTION
Molten Salt Reactors (MSRs) are a type of nuclear reactors with specific features being the most important one that it works with liquid fuel.This characteristic brings several challenges to their development.That said, its feasibility was proved by the research conducted in Oak Ridge National Laboratory (ORNL) during the 50s and 60s, resulting in the construction and operation of the Molten Salt Reactor Experiment (MSRE) [1].
Recently, in the Generation IV International Forum (GIF-IV) the Molten Salt Fast Reactor (MSFR) was selected as a reference of reactors with circulating fuel.The EURATOM EVOL (Evaluation and Viability for Liquid Fuel Fast Reactor System) Project was responsible for the development of this project [2].
The MSFR is intended to produce 3 GWth using a fast neutron spectrum.The fuel is a melted salt composed by lithium fluoride salt with addition of fissile and fertile isotopes.The plant is being projected to be composed by three major circuits: the fuel circuit, which is the core cavity with the fuel and fertile salt and the adjacent pumps and heat exchangers (16 loops); the intermediate circuit, which is responsible for extracting the heat generated in the fuel circuit; the power conversion circuit, that is in charge of the electricity generation [3].The conceptual design of the MSFR is presented in Figure 1.Source: [3] Most of the codes currently used for neutronic were developed for solid fuel reactors, therefore, its usage in liquid fuel reactors requires systematic testing and assessment.There are codes that solves neutron distribution using Finite Volume Method (FVM) or Monte Carlo Neutron Transport (MCNT).MCNT codes are stochastic, which implies that the result, will fluctuate around a value, being the magnitude of this fluctuation dependent of the sample size.Generally, MCNT codes are developed to work with Constructive Solid Geometry (CSG) [4].
Recent developments in MCNT codes made some of the available codes work with mesh based geometries [5,6].This feature must be systematically addressed due to the impact in a modeled reactor.The impact is expected to happen mainly due to the approximations done by the calculation methods on curved regions in the geometry [7].Therefore, a grid sensitivity evaluation ensures the quality of the achieved result [8].
The use of the meshed geometries by MCNT codes are usually to perform coupled calculations with thermal-hydraulics.However, there are geometrical complexities in several parts of nuclear reactors that make their modeling with built-in CSG capabilities in MCNT codes difficult, if practical.That happens because of the logical aspect of building a complex CSG model, because these models are constructed from Boolean operations.Perform Boolean operations from basic shapes (lines, planes, cylinders, etc) to build complex geometries tends to consume a lot of time.
Therefore, using mesh based geometries in MCNT codes are an asset for modeling complex geometries.
In this context, the concept of mesh quality becomes relevant.Mesh Quality is related to the mesh configuration that allows a particular simulation of a numerical partial differential equation (PDE) to be efficiently performed, with fidelity to its physics [9].Therefore, a poor mesh stands for meshes with low elements quality (small or large angles, poor shape) and low quantity of elements to capture a specific phenomenon, for example, viscous layer effects.However, MCNT simulations do not solves any PDE, but the concept of mesh quality is still relevant.Regarding MCNT mesh based calculations, a mesh is poor when is not able to capture the curvatures of the desired geometry, which implies in volume differences between the discretized (meshed) and continuous domains.This volume differences causes mass differences as well, inserting error in the calculations.
The MCNT mesh based simulations do not crash due to poor mesh.This leads to the utilization of coarse meshes for the sake of decreasing computational demand [10].Coarse meshes in curved shapes insert mass error in the calculations regardless of deterministic or stochastic approaches [7,8,11].These errors will unintentionally give discrepancies in the values of the reactor state (effective multiplication factor -keff) of the system [12].Hence, the mesh based MCNT studies must be performed under a mesh sensitivity evaluation, that is, the level of refinement required where the mass error can be neglected.
The lack of sensitivity evaluation studies on MCNT calculations makes the researchers to perform their MCNT calculations neglecting mass errors due to mesh utilization, since the available works are related with deterministic methods [8,11,13].
The main purpose of the present work is to assess the impact of mesh usage by MCNT simulations with their respective uncertainties discretization.In order to perform the evaluation, the results from MSFR Neutronic Benchmark [3] is compared with the calculated data in the present work.The numerical uncertainty given by the mesh utilization is evaluated using an extended version of the Grid Convergence Index (GCI) [4].At first, calculations with CSG model is conducted to serve as baseline for the developed mesh based models.The fuel salt reprocessing is evaluated with a constant reprocessing rate and its influence due to of the mesh usage.The mass differences inserted by mesh utilization are presented with their impact on the fuel salt inventory.

THE MSFR NEUTRONIC BENCHMARK
The MSFR Neutronic Benchmark was performed by six groups from the Euratom FP7 project EVOL (Evaluation and Viability of Liquid Fuel Fast Reactor Systems) and the ROSATOM project MARS (Minor Actinides Recycling in Molten Salt) [3].Multiple codes and cross section libraries were used in the calculations.The benchmark had the main goal of comparing several neutronics results using different codes in static and evolution simulations.
The geometry used in the Benchmark and its respective dimensions are presented in Figure 2.
The proposed model by the MSFR Neutronic Benchmark is a cylindrical reactor that does not separate the 16 pump and heat exchanger loops [3].A simple geometry was chosen to reduce computational time and allow strait forward comparisons of the numerical results from different groups.
The core has a cylindrical shape with three volumes comprising the active core, the upper and lower plenum.The fuel salt considered in the benchmark is a binary salt.Two compositions of salts were tested, 233 U-started salt and transuranic (TRU)-started salt.In this study we focus on the 233 Ustarted salt composed of Li-232 ThF4-233 UF3.Radial reflectors are composed of Li-232 ThF4 fertile material, that operate as a blanket with an external layer of 20 cm thick B4C that protects the heat exchangers.

Materials initial compositions and properties
The fuel salt composition contains 233 U and 232 Th as fissile and fertile materials, respectively.
The initial fuel salt composition can be seen in Table 1.The proportion of the heavy elements is presented in Table 2.The MSFR Neutronic Benchmark gives several thermophysical properties of the fuel and blanket salts [3].However, only the density (ρ) is needed for reproducing the benchmark results.
This property is required for the calculation of the thermal feedback coefficient.The relation for ρ variation with temperature is given by Equation 1.For nominal temperature (900 K), the fuel salt density is 4.1249 g/cm 3 .
In this equation, ρ is given in g/cm 3 and the temperature (T) in K.The validity range for this equation is between 893 and 1123 K [3].
The blanket salt composition is presented in Table 3.The expected thermal behavior of the blanket salt is similar to the fuel salt, therefore, the thermophysical properties were considered the same for both salts [3].

22.5
The composition for the reflector is shown in Table 4 and the used density was 10 g/cm 3 [3].The composition of the neutron protection is given in Table 5.The used protection density was 2.52 g/cm 3 [3].

Geometry and mesh
In the present work, the MSFR geometry is represented using two different approaches: Constructive Solid Geometry (CSG) and Boundary Representation (B-REP).The CSG geometry model is constructed using the Serpent Nuclear Code (version 2.1.32)built-in CSG capabilities [14].The B-REP geometry is generated using a CAD software.A B-REP geometry must be converted into a volumetric mesh composed of a finite number of elements to be used in the simulation.That way the B-REP approach is Mesh based and will be referred to in this fashion henceforward.
The B-REP model and meshes are generated using GMSH code (version 3.0.6)[15].GMSH is a free finite element mesh generator constructed to build solid geometry representation from B-REP.
The used CSG geometry in the present work is presented in Figure 3.This developed CSG model is used for establishing the comparison baseline with the benchmark (validation), the steady state composition and the burnup simulations.After the CSG model validation, the mesh base models follow the same simulation procedure.The mesh based results are compared with the results from the CSG model and the benchmark.The parameters of the used meshes are shown in Table 6.In this table it can be seen the decrease of the representative mesh size (h) with the increment of elements number.Also, the mesh refinement ratio (r) is shown.The parameters h and r are calculated using the Equations 2 and 3, respectively.shown in Table 6.The generated meshes are presented in Figure 4.In this figure, meshes 3, 2 and 1 are shown and a volume round off is clearly seen between mesh 3 and the other ones.
The volumes of the different components used in the developed models can be seen in Table 7.
In this table it can be observed that as the mesh is refined, the fuel salt volume of the meshed models gets closer to the CSG model.However, the reflector shows a different trend, since more material is present with mesh refinement.From this table is possible to see the mass differences induced by the discretization procedure.For instance, the difference of the fuel salt volume between CSG and Mesh 1 is 0.01 m 3 , which gives around more 40 kg of fuel salt to the CSG model.

The extended GCI method
In order to quantify the numerical uncertainty from the mesh usage, the Grid Convergence Index (GCI) method is applied in an extended mode.This method was originally designed for Computational Fluid Dynamics (CFD) calculations due to the FVM requirement to work with meshed geometries.However, as the complexity of the geometry rises, built-in capabilities of MCNT codes for geometry generation become unpractical, leading to other geometric approaches, such as meshed geometries.GCI method was extrapolated for mesh based Monte Carlo simulations previously, leading to promising results [4].Therefore, this work presents an application of the GCI utilization for MCNT mesh based calculations using MSFR.A detailed explanation of the GCI method is beyond the scope of the present work, nevertheless, it can be found in several papers [16][17][18][19][20].
The GCI method relies on a generalized Richardson Extrapolation method, which assumes that discrete solutions (γ) holds a power series representations in successive mesh refinements based on a representative mesh size (h) [17].The used value of the safety factor (Fs) was 1.25, which is the recommended in a study with at least three meshes [16].The Equation 4 is used for the calculation of the GCI values.In this equation, p is the apparent order of mesh convergence and ε is the absolute difference of a chosen variable values between subsequent mesh results, r is the refinement ratio between successive meshes.The subscript of ε21 and r21 p means that these quantities are relative to the meshed models 2 and 1.
The chosen parameters to be evaluated using extended GCI were the keff and the fuel salt composition evolution.At the end of extended GCI method, the 95% confidence interval of the chosen parameters is obtained for the finer mesh (in the group of three meshes).
The GCI method provides an extrapolation alternative to calculate the parameter evaluated at an infinite discretization level (h = 0).This can be performed using Equation 5.In this equation, φ 21 ext is the extrapolated value of the chosen parameter and φn is the value of the chosen parameter for the correspondent refined mesh (nth meshed model).
The MCNT simulations provide a 95% confidence interval regardless of the used geometry (mesh or CSG).This is accounted as an additional uncertainty component and combined with the GCI value in the numerical uncertainty using the root sum of squares method.

Fuel reprocessing
The operation of the MSFR produces several undesired fission products (FP) inside the core (such as xenon and krypton), if uncontrolled.The presence of such FP adversely affect power production.In nominal operation condition, these FP form and decay at the same rate.However, according to the electricity demand, the power production can vary.In a decreasing power shift, the FP starts to build at a higher rate than they decay, accumulating inside the reactor.After a while, if power production must be raised, the response rate is impaired by the time required for the FP to reach their equilibrium level.Hence, buildup of these undesired FP should be avoided according to Doligez et al. [21].Opportunely, the fuel salt reprocessing is responsible for extraction of unwanted FP and insertion of more fissile/fertile materials [22].
In this context, molten salt reactors modelling requires specific features in the existing depletion codes.These codes should be capable of capturing the online fuel reprocessing and refueling.Most contemporary nuclear reactor physics software do not have these features, since they were originally designed to simulate solid-fuel reactors [23,24].However, Serpent Nuclear Code has this capability of online reprocessing of the used materials, limited to constant reprocessing through the depletion time [25].
Fortunately, Serpent Nuclear Code has a capability that allows liquid-fuel reprocessing and refueling.This feature has the limitation of a constant insertion/removal rate per isotope.The time constant (λ) for each element is given in the input.
The λ values for each element removed and inserted in the core is presented in Table 8.The used reprocessing parameters (λ values, day steps, reprocessing materials composition) is the same for the developed CSG and mesh base models.This implies that the long term impact of discretization level can be assessed.

RESULTS AND DISCUSSION
The first calculation in the MSFR Neutronic Benchmark was a criticality one, given an initial composition for fuel salt and the other materials of the reactor (see Tables 1, 2, 3, 4 and 5).The extracted results from these initial calculations were the effective multiplication factor, delayed neutron fraction, generation time and thermal feedback coefficients.
Additional calculations took place using the initial composition as a start and the aim was to find the fuel salt composition which resulted in steady state behavior (keff = 1).These additional calculations were performed in order to establish the base composition for depletion calculations, and in this case, steady state was necessary to allow the long operation time used in the benchmark (200 years).The results from the depletion calculations are the heavy nuclei of the fuel salt inventory.
The MSFR Neutronic Benchmark was performed by different groups from the Euratom FP7 project EVOL.Multiple codes and cross section libraries were used in the calculations.The calculations presented in this paper are compared with the Benchmark results using the same cross section library (ENDF/B-VII.1).From now on, we refer to the benchmark results regarding the specific results using ENDF/B-VII.1.

Steady state calculations
The first criticality calculations are performed to find the fuel composition which resulted in a steady state condition (keff = 1).These initial calculations are performed using the fuel salt composition shown in Table 2.The keff results from these calculations and the MSFR Benchmark are shown in Table 10.In this table, the uncertainty values are related with the MCNT component, neglecting the discretization effect.The CSG model result presented good agreement with the benchmark one, being validated.The results from the mesh based models showed consistency between each other, that is, as the mesh refinement increases, the result approximates to the CSG one.The results from the meshed models showed the expected discrepancies in comparison with the CSG model and the benchmark, demonstrating the effect of the discretization in the materials volume.
The results shown in Table 10 are complemented in Figure 5.In this figure, the developed CSG and mesh based models regarding keff and the discretization level are shown, with the combined uncertainty values (discretization and MCNT) and the extrapolated value for infinite discretization.
The uncertainty bars are calculated from the squares root sum of MCNT uncertainty plus GCI method.The results of mesh based models shows an asymptotic behavior, being suitable for the calculation of the extrapolated result of GCI method (h = 0).
As described for Table 10, as the mesh is refined, the keff values shows a clear trend to reach the CSG value.The most refined meshed model (Mesh 1) shows no statistically significant difference with the CSG model, hence being validated.The same behavior is observed with the extrapolated GCI value (h = 0).Nevertheless, the uncertainty bars show that the contribution of the discretization component is less impacting than the MCNT one, since GCI was calculated only for the two finest points (lower h) and the magnitude of the uncertainty bars are equivalent for all points.
It must be noted that GCI method suited well the calculations of discretization uncertainty for keff, which was expected because of the literature results [4].After the initial criticality calculations, an iterative process was carried and the steady state composition was found.The achieved proportion of heavy elements for the developed models and the benchmark one can be observed in Table 11.The results from CSG model and the meshes 1 and 2 showed good agreement with the benchmark, being validated.In this table it can be seen the decrease of 233 U with mesh refinement, which is an indication that a more refined geometry gets closer to the actual geometry (CSG model) in terms of volume differences.In order to keep the proportion of the fuel salt, as 233 U decreases, Th quantity is incremented, and the other salt components are kept constant.
However, meshed model 1 result showed less 233 U than the CSG model.This could be due to the volume differences on the other materials of the reactor, such as the reflector and the blanket (see Table 7).The meshed model 1 showed a higher volume of the reflector and the blanket, while the density was kept the same for both CSG and the meshed models.This causes a virtual enhanced in the density of these materials, leading to more neutrons to be in the fuel salt cavity, instead of interacting with the other materials.The β0, βeff and generation time given by the MSFR Benchmark and its comparison with the developed models for steady state composition are presented in Table 12.The achieved results showed reasonable agreement with the benchmark ones, however the uncertainty values are not given in the MSFR Benchmark for the used cross section library.The calculated reactivity coefficients of the developed models are presented in Table 13.The results from the developed models showed good agreement in comparison to the benchmark.

Burnup calculations
In general, the achieved results for steady calculations are in good agreement with the ones available in the MSFR Benchmark, being validated.In order to fully assess the reliability of the developed models, the burnup calculations are performed and compared to the benchmark.
Regarding the keff, the MSFR did not present its evolution with time, since this parameter was kept constant due to the variable reprocessing scheme [3,29].However, the evolution with time of the keff is presented in Figure 6.In this figure the uncertainty bars are given only by the combination of the components from the extrapolated GCI method and from the MCNT standard deviation.The results from CSG and mesh based models showed no statistically significant difference.In terms of behavior, it can be notice an overshoot in the beginning of the operation time.This behavior is expected and reported by Aufiero [29] when variable reprocessing is not applied.Nevertheless, after approximately 50 years the keff reaches a constant value.The results showed that even with constant reprocessing parameters, the reactor tends to reach equilibrium condition.
The reprocessing scheme in the present work is constant for the fuel salt and the blanket salt, and even so, an equilibrium condition is achieved.In the MSFR Neutronic Benchmark, the reprocessing scheme changed over time, although keff was not shown.The benchmark reprocessing scheme was intended to control the reactivity, keeping keff practically constant [3,30].Hence, the results differences accumulate between the present work and the benchmark, in terms of the burned materials.The time evolution (burnup) of the fuel salt inventory is presented in Figure 7 for the benchmark, CSG and mesh based models (Mesh 1).In this figure it can be seen that the developed models lead to mass evolution with good similarity to the MSFR Benchmark results.
The mass of total Uranium remains almost constant for all simulation time, which denotes that the chosen reprocessing scheme is effective, although it is constant.The developed models showed good agreement for the evolution of mass of total Uranium.
The mass of Pa isotopes is also kept almost constant the entire time.A small difference of Pa isotopes mass evolution between the developed models and the benchmark can be observed through almost all accounted time, which is due to the difference in the reprocessing scheme.
The masses of Pu, Np, Am and Cm isotopes are built slowly in comparison to Uranium and Pa isotopes, which is expected given its decay path.There are discrepancies related to the number of time steps for these isotopes in comparison to the benchmark.This is due to the fact that the present work used more time steps than the benchmark.However, the final masses show good agreement between the benchmark and the developed models (CSG and mesh based).These discrepancies can be observed between 15 and 50 years of operation for Pu and Np isotopes.For Am and Cm isotopes the differences can be seen around 150 of operation time.
The differences observed between the developed models and the benchmark are probably due to the fuel reprocessing differences as attested earlier.The extended GCI method produced such small values that it cannot be seen due to the magnitude of the elements masses.This happens due to a strong convergence of the method and small differences between the results of the meshes, meaning that the results are in consistent between each other.Besides, the method was stable for all the elements and all operation time of the Figure 7.
The fuel salt inventory for all Uranium isotopes is shown in Figure 8.In this figure, the developed models and the benchmark are in good agreement with small discrepancies.The initial amount of U isotopes was smaller for the developed models; however, this was changing through the simulation time.At the end, the developed models (CSG and meshed) accumulated more Uranium mass than the benchmark.Even with the differences between the developed models and the benchmark, the general behavior of U isotopes was captured.The extended GCI method was applied and produced such small values that the vertical scale is unable to highlight these uncertainty values.Again, this occur because the results are in asymptotic convergence and presents negligible results differences between the successively refined meshes.The fuel salt inventory for 233 U is shown in Figure 9. Regarding the developed models, the initial mass of 233 U presented in the benchmark calculation was higher and as for the total Uranium mass (see Figure 8), at the end, the developed models slightly over predicted the 233 U mass.The 233 U mass curve for the developed models shows that around 100 years, the small amount of 233 U accumulated starts to be consumed by the reactor.In general, one can say that 233 U is kept almost constant during the evaluated operation time, which means that the reprocessing scheme kept the reactor close to a steady state condition.The differences between the reprocessing schemes of the developed models and the MSFR Benchmark causes the discrepancies between the results, since the present work used constant reprocessing values and the benchmark changed these values to keep keff constant.The extended GCI method showed the same behavior for 233 U isotope and total Uranium (see Figure 8), due to minor differences between the results of the meshes and strong convergence of the method.The mass evolution of 231 Pa is presented in Figure 10.In this figure, the results from the developed models and the benchmark ones are in good agreement, except for the last two time step points.In these points, there is a discrepancy between the presented models and the benchmark.It seems that 231 Pa starts to be consumed at the end of the simulation time due to the reprocessing differences.The extend GCI method was convergent and consistent for the mass evolution of 231 Pa isotope presented in Figure 10.However, as it happened for 233 U isotope and total Uranium, the The Figure 11 shows the mass evolution of Pu isotopes in the fuel salt inventory.The results from the developed models are in good agreement with the ones from MSFR Benchmark for all Pu isotopes.This was expected, since in Figure 7, the sum of Pu isotopes was practically identical between CSG model and the benchmark.The extended GCI method produced small convergent and consistent values for the mass evolution of all Pu isotopes presented in Figure 11.This happens due to negligible differences between the results of the meshes and asymptotic convergence.In general, the main difference between the results of the developed models and the benchmark was in the range of the total Uranium mass, the 233 U mass and the last two compared points for the 231 Pa mass.This is due to the reprocessing scheme, even with the reactor reaching a stable state (see Figure 6).The GCI method is convergent and consistent for the mass evolution of all monitored elements presented in the present work.

CONCLUSIONS
In the present paper, an assessment of the impact of MCNT mesh based simulations and its numerical uncertainty are performed comparing the achieved results to the data available in the MSFR Neutronic Benchmark.The developed CSG model is used as a baseline to the mesh based models.The fuel salt reprocessing scheme is evaluated using a stable reprocessing rate and its influence in the developed models is discussed.The numerical uncertainty given by the mesh utilization is evaluated using an extended version of the Grid Convergence Index (GCI), which is combined with the MC uncertainty.
The developed models are validated against the benchmark results.The achieved results showed consistent agreement with the benchmark, regardless of the differences in the reprocessing scheme.
Regarding the calculations using mesh based models, the main parameters affected by the mesh refinement are the fuel salt proportion of heavy elements and the keff.This sets a flag on the current utilization of MCNT simulations with mesh based models, since generally no mesh sensitivity is performed, potentially leading to erratic values.
Concerning the mesh based models, they captured the behavior of the reactor, as their validation against the benchmark and the CSG models showed.However, the mesh impact in the initial composition is clear because of the differences induced by the discretization process.
The reprocessing scheme seems to be the key parameter to reproduce integrally the benchmark results.However, if the constant reprocessing parameters are well calibrated it is possible to partially reproduce the benchmark, in terms of mass evolution of the monitored elements, as presented here.The used version of Serpent Nuclear Code does not allow variable reprocessing schemes through time.

Figure 2 :
Figure 2: The dimensions of the geometry used in the MSFR Neutronic Benchmark (mm).

Figure 3 :
Figure 3: Side and half top view of the used CSG geometry.

Figure 4 :
Figure 4: Side and half top view of the used CSG geometry.

Figure 5 :
Figure 5: keff evolution with discretization level for the initial composition.

Figure 6 :
Figure 6: Evolution of the keff with time.

Figure 7 :
Figure 7: Time evolution of the trans-Th elements.

Figure 11 :
Figure 11: Evolution of the plutonium isotopes in the fuel salt.

Table 2 :
[3]tial composition of the heavy elements in the fuel salt[3].

Table 4 :
[3]based alloy composition used in the structural material of the core[3].

Table 6 :
Parameters of the developed mesh based models.

Table 7 :
Volumes of the components of the developed models.

Table 8 :
[26]ocessing/refueling table[3].In the present work, the used daysteps setup is presented in Table9.For the static composition criticality calculations, the number of neutron histories is based on Fiorina et al.[26], being 2.5E7,