A stochastic method to propagate uncertainties along large cores deterministic calculations

Deterministic uncertainty propagation methods are certainly powerful and time-sparing but their 
access to uncertainties related to the power map remains difficult due to a lack of numerical convergence. On the 
contrary,stochastic methods do notface such an issue and they enable a more rigorous access to uncertainty related 
to the PFNS. Our method combines an innovative transport calculation chain and a stochastic way of propagating 
uncertainties on nuclear data: first, our calculation scheme consists in the calculation of assembly self-shielded cross 
sections and a pin-by-pin flux calculation on the whole core. Validation was done and the required CPU time is 
suitable to allow numerous calculations. Then, we sample nuclear cross sections with consistent probability 
distribution functions with a correlated optimized Latin Hypercube Sampling. Finally, we deduce the power map 
uncertainties from the study of the output response functions. We performed our study on the system described in 
the framework of the OECD/NEA Expert Group in Uncertainty Analysis in Modelling. Results show the $^{238}$U 
inelastic scattering crosssection,the $^{235}$U PFNS,the elastic scattering crosssection of $^1$H and the $^{56}$Fe crosssections as major contributors to the total uncertainty on the power map: the power tilt between central and peripheral 
assemblies using COMAC-V2 covariance library amountsto 5.4% (1$\sigma$)(respectively 7.4% (1$\sigma$) using COMAC-V0).


Context: GEN-III cores
The improvements in reactor technology of the so-called GEN-III reactors are intended to result in a longer operational lifetime (at least 60 years) compared with currently used GEN-II reactors (designed for 40 years of operation). In particular, they take advantage from a simpler and more rugged design, making them easier to operate and less vulnerable to operational upsets. From a neutronic point of view, a higher burn-up is aimed to reduce fuel consumption as well as the amount of corresponding waste. More specifically, we will study advanced GEN-III cores which are bigger than current PWRs and characterized by a heavy reflector.

Objectives
Whilst powerful and timeÀ sparing methods have been largely used to propagate uncertainties in core calculations [1], a great endeavour is made to brush up on brute force methods. Actually, thanks to growing calculation means, stochastic methods become more attractive to calculate the uncertainty introduced by simulation codes [2]. These methods consist in taking into account the uncertainties either from the very beginning of the calculation chain [3,4] by sampling nuclear model parameters or by sampling nuclear cross sections with consistent probability distribution [5,6] for burn-up calculations [7], for the TMI core power map calculations [8] with different uncertainty propagation systems [9,10]. Having regard to the nuclear model parameters uncertainty ranges one can produce random nuclear data evaluation (for example with TALYS [11]). Finally, the uncertainties are deduced from the study of the output distribution functions of interest.
Here, we propose a similar method which combines an innovative calculation chain and a stochastic way of propagating uncertainties on nuclear data. Given that large cores are more sensitive to a modification of the neutronic balance (for example, a local modification of the moderator properties), the uncertainties associated with calculation parameters are worth studying. We propose then here to propagate uncertainties due to nuclear data on LWR key parameters such multiplication factor and core power map.

Theoretical background
The Boltzmann equation, which translates the neutron balance in a nuclear reactor, can be written in a compact form as where A 0 is the disappearance operator, F 0 the neutron production operator, ' 0 is the unperturbed neutron flux and l 0 = 1/k 0 with k 0 the first eigenvalue associated with the fundamental mode flux ' 0 . Theoretical arguments lead to express a perturbed flux sensitivity coefficient a 1 as following [12]: where l 1 = 1/k 1 is the inverse of the first harmonic eigenvalue of the flux, ' the neutron flux and ' † the adjoint neutron flux. While the scalar product is dependent on the external perturbation only, the first term is directly related to the size of the core. 1=ðl 1 À l 0 Þ is called the eigenvalue separation factor (EVS) and grows like the size of the core. Thus with the same initial perturbation, the larger the core, the higher the resulting flux perturbation amplitude. The deterministic nuclear data uncertainty propagation on a manifold sample of french PWRs (from 900 MWe to 1700 MWe) showed that the central assembly power uncertainty increases from 1.5% to 4% (1s) [13] mainly due to the uncertainty on 238 U inelastic scattering cross section.

The core calculation scheme 2.1 Physical model
Our model is based on a realistic pattern proposed in the framework of the Uncertainty Analysis in Modelling Expert Group of the OECD/NEA [14]. An international numerical benchmark has been proposed to study the uncertainties related to large cores: a fresh core with 241 assemblies, each of them containing 265 pins at hot zero power (HZP). Under these operating conditions, no thermohydraulics feedback flattens the power map. Thus the HZP radial power uncertainty is expected to reach its maximum. A whole 2D core calculation is undertaken, with a refined pin description and a flux calculation scheme in two steps. First, each individual assembly geometry in an infinite lattice is described. After that, self-shielded 281 SHEM [15] energy groups cross sections are produced. Then, the neutron flux is calculated thanks to the method of characteristics onto the whole geometry. Even though the computational power has been steadily growing with time, yet the CPU time needed in order to have flux convergence is still too high. That is why several assumptions are made in order to reduce CPU time cost. Given that the steady state Boltzmann equation is discretized in space and energy, we decided to vary the calculation parameters from APOLLO2.8 reference calculation scheme used at CEA [16]: -We present in Figure 1 the spatial mesh of our study. The interface between the reflector and the peripheral assembly has been finely meshed to keep a good description of the impact of the reflector on the power map. Due to this refinement, the distance between each characteristic was settled to Dr = 0.04 cm.
-We used a 20 groups energy mesh [16] (cf. Tab. 1) instead of 26 groups [16]: we removed the six groups devoted to the description of the 240 Pu 1 eV resonance, that is not needed in LWR-UOx cores. -In order to obtain an accurate neutron migration in the core, a P3 Legendre expansion was used to describe the anisotropic scattering.
Our method allowed to reduce the CPU cost for a pinby-pin transport core calculation from 1 day to 1 h, on a single Intel 3 GHz processor with 11Go of RAM use.

The cross section sampling method
The values of our input nuclear data are assumed to be described by a Gaussian distribution with the standard deviation given by the covariance library. Given that the number of calculations is limited, the population of our statistical sample must be small. This so-called design of experiments must be wisely chosen in order to fulfill the three following properties : optimal covering of the input parameter space, robustness of projections over 2D subspaces and sequentiality. Now, we will show that the Latin hypercube sampling (LHS) is the optimum fitting to our need. The LHS comes down to equally chop off all the dimensions, and thus make sub-intervals of equal bin width. Each sampling point coordinate is the only one in each sub-interval. To get the best compromise between a not large confidence interval and a limited population, we chose a population of 50 with a L2-star optimized LHS [17]. Once we sampled each cross section in each energy group as a Gaussian N ð0; 1Þ (the corresponding vector is noted X), covariance libraries have to be taken into account. This was made by using a single value decomposition. Here we note the covariance matrix C, the dimension of our problem d, the vector containing all the means m. We are looking for a multivariate Gaussian vector Y ∈ ℝ d , whose probability density function (pdf) is N ðm;CÞ. Let us note Y = QX+m. The problem is actually equivalent to choose Q so that EððY À mÞðY À mÞ T Þ ¼ C. We took eventually Q as Q = PD 1/2 where D is the eigenvalues diagonal matrix of C and P the corresponding transfer matrix. In the framework composed by the eigenvectors, the linear application corresponding to D 1/2 is a dilatation of the distribution and P corresponds to a rotation. Then the distribution is shifted by m.

Results
The cross-section values calculated after the self-shielding step are modified according to Y. Then, for each sampling, our calculation routine is run. Finally we study the final distribution function of the multiplication factor and the one of the assembly power map in order to deduce the related uncertainties. The Gaussian output function is infered by fitting the probability plot, as shown in Figure 2. For this work, we chose to use two sets of covariances matrices taken from COMAC: the first one, called here COMAC-V0 [18], was released in 2012 and the second one, called COMAC-V2 contains major results obtained until 2016 [19,20]. Thus we compare the impact of the two covariance libraries on the power map and highlight how the library change has contributed to reduce the contribution of several major isotopes to the total uncertainty. Table 2 spots the seven major contributions to the total uncertainty on the k eff , the center of the power map, and the peripheral assemblies. 235 U x and 238 U (n,n') contain the highest uncertainties of the power map.
By combining quadratically the major contributions, we obtain for the k eff a total uncertainty of 669 pcm. On the plus side, we decided to take into account more correlations between cross sections thanks to a dedicated design of experiments for the total uncertainty. So, the corresponding uncertainty raises up to 688 pcm.
Concerning the center of the power map and the outer ring of assemblies these total uncertainties stand respectively for 4.2% and 3.2%.
Similarly, Table 3 presents the last results obtained with the new covariance library, COMAC-V2: -The contribution due to the 238 U fission cross section has dramatically been reduced: above 1 MeV, the standard deviation in COMAC-V2 is around 2%-3%, the same order of magnitude as the standards. That leads to a reduction of its contribution by a factor of 4 on the k eff . -Concerning the contribution of the inelastic cross section of 238 U: in COMAC-V0, the uncertainty value on the plateau is around 10%.
In COMAC-V2, the uncertainty value above the threshold has been strongly reduced to 5%-6%. However, we assume this value to be optimistic. For further calculations, we propose an uncertainty level of 15% on the plateau. This would include an overestimation of the cross section value in JEFF-3.1.1 by 10% on the plateau. Even so, we point out that these values stay much lower than the ones given example by ENDF/BVII.1 (30%).
-The contribution of the 235 U PFNS uncertainty has been reduced from COMAC-V0 to COMAC-V2: the level of the input uncertainty has been reduced by more than a factor of 2 on the whole energy range, reducing drastically the uncertainty on the power map. -Interesting is the contribution of the iron scattering in the power map uncertainty. Given that most of the iron is contained in the heavy reflector, an increase of the iron scattering cross section will enhance the ability of the reflector to send back neutrons in the core. Then, at the periphery of the core more neutrons will be moderated and more fissions will occur. Finally, this will add a radial swing to the power map. That is why this cross section  held a lot of attention with the PERLE program [21,22] which helped to produce the covariances included in COMAC V0.
Finally, Figure 3 presents the whole uncertainty map on our core. It clearly shows the radial swing between the centre and the outer ring of assemblies.  We proved that uncertainties on light water reactors parameters due to nuclear data can be propagated through a brute force method thanks to modern computation power. This method gives access to all of the needed uncertainties without developing any dedicated perturbation theory or using special hypothesis. We applied this method to our PWR large core NEA benchmark and showed that the overall k eff uncertainty reaches 634 pcm, 3.1% for the centre of the power map and 2.3% for the outer ring of assemblies, thus a potential power swing of ± 5.4% (1s). The main contributors are n, the capture and fission cross sections of 235 U, the capture cross section of 238 U for the k eff . Inelastic scattering cross section of 238 U, PFNS of 235 U, elastic scattering cross section of 56 Fe and 1 H are the main contributors to the assembly power map uncertainty.
This method could be applied to propagate other uncertainties, especially design and technological uncertainties, whose analytical expressions are difficult to derive with the usual perturbation theory.