# SUSD3D Computer Code as Part of the XSUN-2017 Windows Interface Environment for Deterministic Radiation Transport and Cross-Section Sensitivity-Uncertainty Analysis

Authors: Ivan A. Kodeli; Slavko Slavič
Journal: Science and Technology of Nuclear Installations
 (2017)

---

## Abstract

A Windows interface XSUN-2017 facilitating the deterministic radiation transport and cross-section sensitivity-uncertainty (S/U) calculation is presented. The package was developed to assist the users in the preparation of input cards, rapid modification, and execution of the complete chain of codes including TRANSX, PARTISN, and SUSD3D, all available from the OECD/NEA Data Bank and RSICC. The objective of this work was to make the input and output handling for these codes as user-friendly as possible, passing information among codes internally. XSUN-2017 allows a user-friendly viewing of results obtained from the PARTISN and SUSD3D programs. The first version of the Windows interface XSUN-2013 was developed in 2013 and submitted to OECD/NEA Data Bank Computer Code Collection and RSICC in early 2014. An updated version, XSUN-2017, will be released in 2017. The package includes also the new version of the SUSD3D code. The XSUN-2017 and SUSD3D code systems and recent improvements and updates are described. Examples of the use and validation are presented, including the S/U intercomparison exercise using the SNEAK-7 benchmarks involving the XSUN-2017 code system comparison with the codes such as TSUNAMI, SERPENT, and MCNP6, and the S/U analysis of thekeff and βeff parameters for the MYRRHA accelerator driven system (ADS).


---

## Body

## 1. Introduction

In order to facilitate the use of the deterministic particle transport codes, in decline as compared to the more and more widely used Monte Carlo codes, the development of a modern Windows interface environment was started in 2012 initially sponsored by the NEA Data Bank. Several of the high quality deterministic neutron transport codes, such as ANISN [1], DORT-TORT [2], DANTSYS [3], and PARTISN [4], date back to the 60s and 70s and have difficulties to meet the expected “comfort of use” of nowadays users and are therefore today less and less used. However, the deterministic codes are still attractive and useful in many possible applications, for fast radiation transport providing a complete flux distributions, and for the cross-section sensitivity and uncertainty analyses, deep penetration problems, and validation of Monte Carlo codes and methods. The motivation for this work was to facilitate the use of deterministic transport codes within a modern environment and thus bring these tools closer to nowadays users.The XSUN-2013 (cross-section sensitivity and uncertainty) [5, 6] code system was released in spring 2014 to OECD Nuclear Energy Agency Data Bank (Package-ID NEA-1882) and RSICC (Code Number CCC-825). The package can be obtained from the two data centers, by contacting the RSICC or NEA Data Bank Computer Program Services or one of the nominated establishments in your country. An updated version of the code system, XSUN-2017, is under final testing and will be released at the beginning of 2017.


## 2. XSUN Computer Code System Description

The XSUN-2017 code package is a user–computer interface environment integrating several deterministic codes and organized in a way to make the input and output handling for these codes as user-friendly as possible, passing information among the codes internally to assure input data preparation in a consistent way. The package is based on the Xbase++ (R) Compiler 1.90.331 and Alaska 32-Bit Linker [7], for the pre- and post-processing of the input and output data. It supports for now a limited, but complete and self-consistent set of deterministic codes, all available from the OECD/NEA Data Bank and RSICC, with the following utility functions (see Figure 1):(1)
Multigroup nuclear cross-section preparation (TRANSX-2-15 [8] code): TRANSX is a computer code that reads nuclear data from a library in MATXS format and produces transport tables compatible with many discrete ordinates (SN) and diffusion codes. MATXS format libraries are prepared using the NJOY-99 [9] code (or more recent versions). Tables can be produced for neutron, photon, charged-particle, or coupled transport. Options include adjoint tables, mixtures, homogeneous or heterogeneous self-shielding, group collapse, homogenization, thermal upscatter, prompt or steady-state fission, transport corrections, elastic removal corrections, and flexible response function edits. The output self-shielded multigroup cross-sections for the mixtures according to the user’s mix instructions can be written in DTF-style card images, FIDO, ISOTXS, or the binary group-ordered GOXS format.(2)
Transport of neutral particles for criticality and shielding calculations (codePARTISN  [4]): PARTISN (PARallel, TIme-Dependent SN) is a discrete ordinates (SN) transport code for shielding and criticality calculations. The code is a modular computer program package designed to solve the time-independent or dependent multigroup discrete ordinates form of the Boltzmann transport equation in several different geometries. The modular construction of the package separates the input processing, the transport equation solving, and the post-processing (or edit) functions into distinct code modules: the Input Module, the Solver Module, and the Edit Module, respectively. In addition to the diamond-differencing method, the Solver Module also has Adaptive Weighted Diamond-Differencing, Linear Discontinuous, and Exponential Discontinuous spatial differencing methods. The spatial mesh may consist of either a standard orthogonal mesh or a block adaptive orthogonal mesh. The Solver Module may be run in parallel for two and three-dimensional problems. One can run 1D problems in parallel using Energy Domain Decomposition. Both the static (fixed source or eigenvalue) and time-dependent forms of the transport equation are solved in forward or adjoint mode. In addition, PARTISN has a probabilistic mode for Probability of Initiation (static) and Probability of Survival (dynamic) calculations. Vacuum, reflective, periodic, white, or inhomogeneous boundary conditions are solved. General anisotropic scattering and inhomogeneous sources are permitted. PARTISN solves the transport equation on orthogonal (single level or block-structured AMR) grids in 1-dimensional (slab, two-angle slab, cylindrical, or spherical), 2-dimensional (X-Y, R-Z, or R-T), and 3-dimensional (X-Y-Z or R-Z-T) geometries.(3)
Nuclear data sensitivity and uncertainty calculations (codeSUSD3D [10, 11]): the development of SUSD3D started in the early 1990s in the scope of the French pressure vessel surveillance programme and the EC fusion project. XSUN-2017 includes the latest improved and extended version of the SUSD3D multidimensional nuclear cross-section sensitivity and uncertainty code, based on the first-order generalized perturbation theory. The code calculates the sensitivity coefficients and standard deviation in the calculated detector responses or design parameters of interest due to the input cross-sections and their uncertainties. Complex one-, two-, and three-dimensional transport problems can be studied. Several types of uncertainties can be considered, that is, those due to (1) neutron/gamma multigroup cross-sections, (2) energy-dependent response functions, and (3) secondary angular distribution (SAD) or secondary energy distribution (SED) uncertainties.The particle transport calculations are done externally using the existing codes (such as PARTISN in this case), and the information on the direct and adjoint fluxes is passed to SUSD3D via the neutron/gamma flux moment files calculated by PARTISN. The sensitivity profiles are folded with the cross-section covariance matrices to determine the variance in an integral response of interest. The input deck is prepared by XSUN-2017 system taking the geometry description from the PARTISN input and the isotopic cross-section composition as defined in the TRANSX input as well as the self-shielding factors calculated from the latter code. 1-, 2-, and 3-dimensional test cases are included in the package (Figure2).More details on the SUSD3D code are given in a separate chapter below.(4)
Several nuclear cross-section and covariance matrix libraries are available with the XSUN-2017 package. The following multigroup cross-section libraries are to be utilized in TRANSX/PARTISN code suite:
(i)
Nuclear cross-section data library in the ECCO 33-neutron energy group structure, in MATXS format: the library includes data for 109 isotopes at different temperatures (300 K, 550 K, 600 K, 650 K, 700 K, and 800 K) and for different self-shielding factors. The data were produced using the NJOY-99 code and are based mostly on the ENDF/B-VII.1 nuclear data evaluation (with few minor dates still from the -VII.0 release).(ii)
FENDL-3.1 [12] 211 neutron group cross-sections can be extracted from IAEA site https://www-nds.iaea.org/fendl3/ for fusion application. The data were processed at the IAEA.The above neutron cross-sections are also available in the GROUPR format to be used in the code SUSD3D for the proper normalization of the sensitivity profiles.(5)
The following multigroup covariance matrix libraries are available with the XSUN-2017 package:
(i)
ENDF/B-VII.1 [13] covariance matrices processed using the NJOY-99 code in the ECCO 33 energy groups(ii)
JENDL-4.0 [14] covariances processed using the NJOY-99 code in the ECCO 33 energy groups(iii)
SCALE-6.0m covariances [15] processed using the ANGELO code [16] from the original 44-groups to the ECCO 33 energy groups(iv)
Detector response functions from the IRDFF evaluation processed using NJOY-99 (33-group structure).The standard thermal + 1/E + fission + fusion shape of the weighting function was used in the NJOY cross-section and covariance matrix processing.Figure 1
XSUN-2017 computer code system.Figure 2
Screen captures of the PARTISN 1-dimensional (a) and 2-dimensional (b) test case running sequences.
(a)
(b)Note that the XSUN-2017 package distributed through OECD/NEA Data Bank and RSICC includes the SUSD3D code and the multigroup cross-sections and the corresponding covariance matrices. The codes TRANSX and PARTISN have to be requested separately and integrated in the code system by the users. Modifications to the TRANSX code are included in the XSUN-2017 package in the form of an UPD update file.


### 2.1. Modifications with respect to the XSUN-2013 Release

XSUN is an active project and is continuously being developed and updated at JSI. Since the beginning of the project in 2012 altogether 154 internal versions of the XSUN package were prepared. Some of the recently implemented features since the XSUN-2013 version include the following:TRANSX code modifications:
(i)
An error was corrected to make the code compatible with modern (and less permissive) FORTRAN-77 compilers(ii)
Memory size was increased for problems with large number of isotopes (parameter “mixs” > 1000)(iii)
Updates allowing the self-shielding print on a separate file to be used subsequently in XSUN-2017 for the SUSD3D input preparation.CROSS-SECTION LIBRARIES:
(i)
The default ENDF/B-VII.1 cross-sections library was updated with new materials, several temperatures, and self-shielding tables(ii)
Other libraries such as FENDL-3 (211 energy groups) were tested.PARTISN input deck interactive preparation utility:
(i)
Corrections for input preparation of 3-dimensional problems (space mash and zones)(ii)
Multilines in case of lines exceeding 80 characters (e.g., specification ofx-y-z mash coordinates, zone assignment)(iii)
External source specifications are allowed.SUSD3D code modifications, corrections, and updates, such as the following:
(i)
Treatment of the self-shielding factors (σ0) of the cross-sections used in SUSD3D for the normalization of the sensitivities: instead of a single fixed energy independent σ0 factor defined by the user for the previous version, the new code now allows the interpolation among the σ0 tables in the file. By default, the energy-dependent σ0 values, as calculated in the TRANSX code, are used allowing the full consistency between the cross-section used in the transport calculation (PARTISN) and for the SUSD3D sensitivity calculations(ii)
A bug in the code for the 3D PARTISN calculations was corrected(iii)
The uncertainty in the secondary angular distributions (SAD) can be calculated using the covariance matrices provided in the MT = 251 file format for the P1 Legendre terms.SUSD3D interactive input preparation: many modifications, corrections, and updates such as the following:
(i)
Preparation of the Overlays 2 and 3 specifying the list of materials and reactions for which the sensitivities and uncertainties are calculated was completely reorganized and automatized using internally the TRANSX and PARTISN input decks. The material and reaction list is constructed according to the selected cross-section and covariance files selected and requires little user’s intervention. The user can select among several temperatures available in the cross-section library(ii)
Energy-dependent self-shielding factors are taken directly from the TRANSX output(iii)
Consistency checks between SUSD3D/TRANSX and SUSD3D/PARTISN inputs(iv)
Major corrections and updates in the space mash conversion from the PARTISN input, in particular for 3D geometries(v)
Multilines in case of lines exceeding 80 characters (e.g., specification ofx-y-z mash coordinates)(vi)
For the visualization of the sensitivity profiles, a link to the user-friendly “Sensitivities Plots” [17] utility was introduced.Many other modifications and improvements were introduced, such as project test case save/open features, file history cleaning button, and so on.


## 2.1. Modifications with respect to the XSUN-2013 Release

XSUN is an active project and is continuously being developed and updated at JSI. Since the beginning of the project in 2012 altogether 154 internal versions of the XSUN package were prepared. Some of the recently implemented features since the XSUN-2013 version include the following:TRANSX code modifications:
(i)
An error was corrected to make the code compatible with modern (and less permissive) FORTRAN-77 compilers(ii)
Memory size was increased for problems with large number of isotopes (parameter “mixs” > 1000)(iii)
Updates allowing the self-shielding print on a separate file to be used subsequently in XSUN-2017 for the SUSD3D input preparation.CROSS-SECTION LIBRARIES:
(i)
The default ENDF/B-VII.1 cross-sections library was updated with new materials, several temperatures, and self-shielding tables(ii)
Other libraries such as FENDL-3 (211 energy groups) were tested.PARTISN input deck interactive preparation utility:
(i)
Corrections for input preparation of 3-dimensional problems (space mash and zones)(ii)
Multilines in case of lines exceeding 80 characters (e.g., specification ofx-y-z mash coordinates, zone assignment)(iii)
External source specifications are allowed.SUSD3D code modifications, corrections, and updates, such as the following:
(i)
Treatment of the self-shielding factors (σ0) of the cross-sections used in SUSD3D for the normalization of the sensitivities: instead of a single fixed energy independent σ0 factor defined by the user for the previous version, the new code now allows the interpolation among the σ0 tables in the file. By default, the energy-dependent σ0 values, as calculated in the TRANSX code, are used allowing the full consistency between the cross-section used in the transport calculation (PARTISN) and for the SUSD3D sensitivity calculations(ii)
A bug in the code for the 3D PARTISN calculations was corrected(iii)
The uncertainty in the secondary angular distributions (SAD) can be calculated using the covariance matrices provided in the MT = 251 file format for the P1 Legendre terms.SUSD3D interactive input preparation: many modifications, corrections, and updates such as the following:
(i)
Preparation of the Overlays 2 and 3 specifying the list of materials and reactions for which the sensitivities and uncertainties are calculated was completely reorganized and automatized using internally the TRANSX and PARTISN input decks. The material and reaction list is constructed according to the selected cross-section and covariance files selected and requires little user’s intervention. The user can select among several temperatures available in the cross-section library(ii)
Energy-dependent self-shielding factors are taken directly from the TRANSX output(iii)
Consistency checks between SUSD3D/TRANSX and SUSD3D/PARTISN inputs(iv)
Major corrections and updates in the space mash conversion from the PARTISN input, in particular for 3D geometries(v)
Multilines in case of lines exceeding 80 characters (e.g., specification ofx-y-z mash coordinates)(vi)
For the visualization of the sensitivity profiles, a link to the user-friendly “Sensitivities Plots” [17] utility was introduced.Many other modifications and improvements were introduced, such as project test case save/open features, file history cleaning button, and so on.


## 3. SUSD3D Code

The 3-dimensional version of the SUSD3D code was developed in the early 1990s in the scope of the European fusion project. The SUSD3D code was distributed through the OECD/NEA Data Bank and RSICC since 2000 as a stand-alone code package NEA-1628 [18]. The last version was released in 2008 (NEA-1628/03). The XSUN-2017 package includes the latest improved and extended version of the SUSD3D multidimensional nuclear cross-section sensitivity and uncertainty code. Based on the first-order generalized perturbation theory [19], the code calculates the sensitivity coefficients and standard deviation in the calculated detector responses or design parameters of interest (keff, βeff, reaction rates) due to the input cross-sections and their uncertainties. Complex one-, two-, and three-dimensional transport problems can be studied. Several types of uncertainties can be considered, that is, those due to the following:(i)
Neutron/gamma multigroup cross-sections(ii)
Energy-dependent response functions(iii)
Secondary angular distribution (SAD) or secondary energy distribution (SED) uncertainties.Either relative σ/RdR/dσ or absolute 1/RdR/dσ sensitivities can be calculated.The particle transport calculations are done externally using the existing codes such as DOORS, DANTSYS, and PARTISN codes, and the direct and adjoint fluxes are passed to SUSD3D via the neutron/gamma flux files calculated by these codes. This guarantees great flexibility in the choice of the transport solver and allows the use of the most up-to-date transport codes. At present, SUSD3D can use the neutron/gamma flux moment files produced by the DORT, TORT [2], ONEDANT, TWODANT, THREEDANT [3], and PARTISN [4] discrete ordinates codes or the angular flux files from the ANISN [1] and DOT-III [20] codes. Extensions to other codes such as DRAGON, ATTILA, and DENOVO are planned.The sensitivity profiles are folded with the cross-section covariance matrices to determine the variance in an integral response of interest. Uncertainties due to the fission spectra uncertainties can be calculated using either the classical or the constrained sensitivity method [21–23], useful particularly in case the fission spectra covariance matrices do not comply exactly with the ENDF-6 Format Manual rules. The code runs under DOS Windows and LINUX using FORTRAN-95 compiler. Several improvements were included in the code since the last stand-alone version of SUSD3D in 2008, such as the following:(i)
Self-shielding treatment: the covariance data are provided for the differential cross-sections. To calculate the corresponding sensitivities, the partial cross-sections at different background cross-sections (σ0) are prepared using the NJOY code. The latest SUSD3D uses the energy-dependent σ0 values (i.e., the same as those used in the transport calculations) instead of an average (effective) σ0 used in the previous versions, which required some expertise and was time consuming for the user. To some surprise, this modification was found to have relatively small impact on the calculated sensitivities and uncertainties comparing to the old method for properly selected average σ0 values; however it greatly simplifies an automatic preparation of the SUSD3D input data.(ii)
Secondary angular distributions (SAD): the latest versions of NJOY-99 and -2012 can process the average cosine covariance data (MT = 251) and recently these data became available (in the ENDF “MF34” format) in some evaluations (JENDL-4.0, TENDL, ENDF/B-VII.1). Note that MT = 251/MF34 data are equivalent to the Legendre P1 term of the elastic cross-section. SUSD3D was therefore updated to read these data. Alternatively, SUSD3D can also read the higher PN terms of Legendre coefficients, such as those prepared in 2000 in the scope of the EFF-2.4 project for the 56Fe, 58,60Ni isotopes. However, these data can only be processed using the SUNJOY code (part of SUSD3D [18]) and no new data in this format were released since 2000.(iii)
SUSD3D was extended to read the direct and adjoint flux moments produced by the PARTISN code.


### 3.1. Special Features Developed in SUSD3D

Several innovative mathematical methods were developed in the scope of the SUSD3D project, and some of them were later introduced also in other S/U computational sequences (e.g., TSUNAMI [24], SERPENT [25], and MCNP [26]) such as the following:(i)
Use of flux angular moments instead of much more voluminous angular fluxes [10, 27](ii)
Demonstration of the equivalence of theconstrained sensitivity method and covariance matrix fix-up in the calculation of the sensitivity to the secondary energy distributions, such as fission spectra sensitivities [21–23](iii)
Sensitivity and uncertainty calculations of theeffective delayed neutron fraction (βeff) [28–30].Future development includes the implementation of the implicit sensitivity method, not yet available in the SUSD3D code and coupling to other transport codes.


#### 3.1.1. Flux Moments versus Angular Fluxes: How Sensitive the Sensitivity Calculations Are on PN Order Approximations?

In the mathematical expression of the nuclear data sensitivity function as derived from the first-order perturbation theory [19] both the angular fluxes and their expansion in the Legendre flux moments are required. Both terms are calculated by the traditional discrete ordinates transport codes and stored in the files. This leaves the developers of the nuclear data sensitivity and uncertainty codes several possible choices of the method to account for the neutron flux anisotropy. Three possibilities available in the SUSD3D code were studied here and compared using the FLATTOP-Pu benchmark exercise. This benchmark was particularly suitable due to its high flux anisotropy.Using the first-order perturbation theory, the relative cross-section sensitivity functionSk,x(E), which expresses the fractional change of an integral response, such as the multiplication factor per fractional variation of a cross-section for reaction x in a particular nuclide k, and in energy group g, is expressed in the following form:(1)Sk,xE=σk,xEkeff∂keff∂σk,xE=1R∫dr⃑∫Ω⃑dΩ⃑ρkr⃑-σk,xEΦr⃑,Ω⃑,EΦ∗r⃑,Ω⃑,E︷loss term+∫Ω⃑′dΩ⃑′∫dE′σsk,xΩ⃑⟶Ω⃑′,E⟶E′Φr⃑,Ω⃑,EΦ∗r⃑,Ω⃑′,E′︷gain term,where R is the normalization factor andρkr⃑ = atomic number density of nuclide k,σk,xE = microscopic cross-section for reaction x, in nuclide k, at energy E,σsk,x(Ω⃑→Ω⃑′,E→E′) = microscopic scattering cross-section from angle Ω⃑ to Ω⃑′, and from energy E to E′, for reaction x, in nuclide k,Φr⃑,Ω⃑,E, Φ∗r⃑,Ω⃑,E = direct and adjoint angular fluxes at position r⃑, with direction Ω⃑ and energy E.The angular dependence of the scattering cross-section is usually described in terms of the cosine of the scattering angle and expanded in Legendre polynomials. In practice the Legendre expansion is truncated after finite number of terms, with L representing here the largest term taken into account. Using in addition the multigroup approximation and finite-difference method the above equation (1) becomes(2)Sk,x,g=1R∑iViρk,i-σk,x,g∑mΦg,i,mΦg,i,m∗ΔΩm+∑g′∑l=0L2l+1σsk,x,l,g→g′∑n=-llMg,il,nM∗g′,il,n.In one- and two- dimensional calculations due to symmetry conditions the summation over n goes only from 0 to l, and in three dimensions all the terms from -l to l must be considered. Note that the (4π) term is included in the integration over angles, that is, ∫4πdΩ=1.In the above equation,Vi is the volume of the space mesh interval i, and Φg,i,m and Φg,i,m∗ are the direct and adjoint angular fluxes, corresponding to energy group g, averaged over space interval i and solid angle ΔΩm. Similarly, ρk,i,σk,x,g,σsk,x,l,g→g′ are space and group averaged quantities defined above. Mg,il,n and M∗g′,il,n are direct and adjoint angular flux moments with indices defined above.We see immediately that both angular fluxes and angular flux moments appear in (2). Although in principle the deterministic transport codes usually calculate and provide on output both quantities, the use of both the angular fluxes and angular flux moments is unpractical and unnecessary, since this requires large computer storage. Indeed, the two quantities are linked by the following expressions: (3)Φg,i,m=∑l=0L∑n=-llYl,nΩm·Mg,il,n,Φg,i,m∗=∑l=0L∑n=-llYl,n+Ωm·M∗g,il,n,(4)Mg,il,n=∑mYl,nΩm·Φg,i,m·ΔΩm′,M∗g,il,n=∑mYl,n+Ωm·Φg,i,m∗·ΔΩm.Several alternative methods are consequently available for the computation of the sensitivity profiles using the above equations (1) to (4), namely,(1)
The loss and gain terms of (3) are expressed in terms of the angular flux moments calculated by the discrete ordinates (SN) transport code; that is, the loss term is calculated by (4):(5)Sk,x,g=1R∑iViρk,i∑l=0L2l+1-σk,x,g∑n=-llMg,il,nM∗g′,il,n+∑g′σsk,x,l,g→g′∑n=-llMg,il,nM∗g′,il,n.Note that thePN expansion order is the same for the loss and gain terms; that is, the one used in the transport calculations. Equation (5) is an approximation of (2), with the two being equivalent if high Legendre order expansion is used for the loss term.(2)
The loss term in (2) is expressed in terms of the angular flux calculated using the SN code and the flux moments needed for the gain term are computed using (4) from these angular fluxes.(3)
Angular flux and flux moments from theSN code output are used to calculate the loss and gain terms of the sensitivity by (3).Most perturbation codes make use of the angular fluxes output from the discrete ordinates transport codes and the angular moments are expressed in terms of these angular fluxes (Method (2)). However, the angular flux files may be in some cases quite bulky. On the other hand, the flux moments files are much more reasonably sized, as shown in Table 1 which compares the number of terms (for each energy group and space interval) required in case of using the angular fluxes and flux moments. Using angular moments (the above method (1)) becomes particularly attractive in three dimensions where, for example, a typically used S16-P3 approximation requires only 16 flux moment terms, but as many as 288 angular flux terms (note that this number is to be multiplied by the number of energy groups and space points to obtain the total number of flux terms needed in a sensitivity calculation).Table 1
Number of flux moment terms used in PARTISN for different Legendre expansion orderL and the number of discrete directions for different SN fully symmetrical direction sets.
Number of flux moments
Number of discrete directions
P0
P1
P3
P3
P3
P5
S2
S4
S8
S16
1D
1
2
4
4
4
6
2
4
8
16
1D-(r, z)
1
2
6
6
6
12
2
6
20
72
2D
1
3
10
10
10
21
4
12
40
144
3D
1
4
16
16
16
26
8
24
80
288An additional advantage of using flux moments, besides a considerable reduction of computer space requirements and CPU time gain, resides in the fact that manySN codes (e.g., ANISN [1]) print the angular fluxes at the space cell edges, and sophisticated negative flux fix-up procedures are subsequently adopted in the transport codes to construct the flux at the center of the cell. Flux moments are conveniently calculated at the center of space mesh.In case of using the flux moments it is on the other hand expected that higherPN Legendre expansion order may be needed as compared to using the angular fluxes, since also the angular fluxes for the loss term are expressed in terms of the angular moments.The SUSD3D code offers the possibility of using the above described options. The mostly used approach is based on (5) (method (1)) calculating the relative sensitivity profiles Sk,x,g from the reasonably sized angular flux moments, produced by the SN codes DORT-TORT, DANTSYS, and PARTISN. Additionally, for one- and two-dimensional analysis, the angular flux files produced by the ANISN [1] and DOT-3.5 [20] codes can be used, thus providing a convenient way to compare and test the results using the two methods.The performance of the above described methods was tested on a highly anisotropic problem while using a relatively low Legendre order. The FLATTOP-Pu benchmark was found to be very sensitive to the scattering anisotropy and as such very suitable for studying the advantages and inconveniences of different approaches. The benchmark, sometimes also referred to as Popsy, consists of ~20-cm natural U reflected 94 wt.%Pu239 sphere, with a radius of 4.533-cm. It is described in the ICSBEP criticality benchmark handbook [31] under the reference PU-MET-FAST-006.FLATTOP-Pu was modeled in a simple 1D spherical geometry using the PARTISN and ANISN transport and SUSD3D S/U codes. S16 fully symmetrical quadratures, ENDF/B-VII.0 33-energy group cross-sections, and the covariance matrices from SCALE-6.1 [15] were used. Legendre expansion orders were varied (P5, P3, and P1). Bell-Hansen-Sadmeister (BHS) transport correction was used for the P1 and P3 calculations.The comparison of the calculated multiplication factors is presented in Table2 showing good consistency between the PARTISN and ANISN, as well as between direct and adjoint results. Small differences may be due to the use of a single fission spectra in ANISN contrary to the zone-dependent one used in PARTISN. BHS correction is shown to improve considerably the results (note ~1700 pcm difference with respect to the uncorrected P1 calculations). P3 results are reasonably close to those obtained using P5.Table 2
Multiplication factors calculated using the PARTISN and ANISN codes using P1 and P5 Legendre expansion. BHS indicates that the Bell-Hansen-Sadmeister transport correction was used. Note an excellent agreement between the direct and adjoint calculations, within few pcm for PARTISN, and within ~15 pcm for ANISN calculations.
Code/options
k
e
f
f-direct
k
e
f
f-adjoint
PARTISN/P5
1.00030
1.00028
PARTISN/P3
1.00013
1.00011
PARTISN/P1
0.99902
0.99900
ANISN/P5
1.00337
1.00320
ANISN/P1
1.00208
1.00193
ANISN/P1 (no BHS correction)
0.98463
0.98448As mentioned, SUSD3D uses angular moment files produced by PARTISN to resolve (3) (Method (1)), or alternatively the angular flux files from ANISN to resolve (2) by Method (2). In case of using the angular fluxes from ANISN, the flux moments are first calculated from the angular fluxes. However, PN order in principle does not need to be the same in the ANISN transport and in the SUSD3D sensitivity calculations. Note that using a high PN order in the loss and the same PN order as in the transport calculation for the gain term is approximately equivalent to Method (3), that is, the exact formulation of the sensitivity as defined in (2).The following computational procedures were intercompared and the corresponding results are presented in Table3 and Figures 3 and 4:(i)
PARTISN/SUSD3D with P5 approximation(ii)
PARTISN using P3 with BHS cross-section correction/SUSD3D using P3(iii)
PARTISN using P1 with BHS cross-section correction/SUSD3D using P1(iv)
ANISN/SUSD3D both using P5(v)
ANISN using P5 and SUSD3D using P1(vi)
ANISN using P1 with BHS cross-section correction/SUSD3D using P1(vii)
ANISN using P1 with BHS cross-section correction/SUSD3D using P5(viii)
ANISN using P1 with BHS correction/SUSD3D using P5 for loss & P1 for gain.The choice of the PN expansion has little impact on the sensitivities of the keff factor to the absorption, fission, and nu-bar data. As expected, the keff factor was found to be particularly sensitive to 238U elastic and inelastic scattering cross-sections and these reactions are also very sensitive to the PN anisotropic scattering approximations; therefore only these results are presented in Table 3 and Figures 3 and 4.Table 3
Comparison ofkeff sensitivity and uncertainty results using different approaches.
Code/options
Sensitivity (%/%)
Uncertainty (%)
238U el
238U inel
238U(el + inel)
Total
PARTISN/P5
0.1399
0.06451
0.406%
1.20%
PARTISN P3
0.1431
0.06776
0.426%
1.21%
PARTISN P1
0.1721
0.09731
0.681%
1.32%
ANISN-P5/SUSD3D-P5
0.1395
0.06325
0.388%
1.19%
ANISN-P5/SUSD3D-P1
0.1710
0.09543
0.654%
1.30%
ANISN P1/SUSD3D-P1
0.1712
0.09550
0.656%
1.31%
ANISN P1/SUSD3D-P5
0.1407
0.064445
0.397%
1.19%
ANISN P1/SUSD3D-P1_gain
0.1180
0.06466
0.633%
1.29%Figure 3
Sensitivity ofkeff to U238 elastic cross-sections calculated by SUSD3D/ANISN codes.P5-P5 = P5 transport & sensitivity calculation,P1-P1 = P1 transport using BHS cross-sections and P1 sensitivity,P1-P5 = P1 transport using BHS, P5 sensitivity calculations,P1-P1 gain = P1 transport using BHS cross-sections, and sensitivity calculation using (2) with P1 gain term.Figure 4
Sensitivity ofkeff to 238U inelastic cross-sections calculated by SUSD3D/ANISN codes.It is important to note however that the BHS correction was not (yet) applied to the cross-sections used in the SUSD3D sensitivity calculations making the two calculations at this stage slightly inconsistent.Several conclusions can be drawn from the results presented in Table3 and on Figures 3 and 4:(i)
Good consistency is observed between the SUSD3D results based on PARTISN angular moments and ANISN angular fluxes, both for P5 and P1 results.(ii)
The analysis of the FLATTOP-Pu critical benchmark revealed large differences in calculatedkeff and the corresponding sensitivity profiles when using low Legendre order expansion such as P1 as compared to the reference P5 (and P3) calculation. The use of BHS correction results in a good agreement with keff factors calculated using the P1 BHS and P5 expansions. It is expected that the relatively large differences between the sensitivities based on P5 and P1 BHS transport calculations could be resolved to large extent by applying the BHS correction to the cross-sections used in the SUSD3D sensitivity code.(iii)
Surprisingly, the use of highPN order for the loss and P1 order for the gain term (SUSD3D-P1_gain), which is equivalent to the “correct” (2) (Method (3)), does not perform necessarily much better than Method (1) (SUSD3D-P1) routinely used in SUSD3D. It can be seen to underpredict the sensitivities to the 238U elastic (Figure 3) and slightly overpredict those to the 238U inelastic (Figure 4) cross-sections. Taking into account the computer space penalties linked to this method and the limited improvements in the results, it seems that the use of (2) (Method (3)) is therefore not justified. Lower PN order calculations would be possible using the angular fluxes and BHS corrected cross-sections, but again at high computer space requirements.


#### 3.1.2. Constrained Sensitivities: Demonstration of the Equivalence of the Constrained Sensitivity Method and the Covariance Matrix Fix-Up

The SUSD3D analysis revealed several differences between the results based on different prompt fission neutron spectra (PFNS) covariances, in some cases resulting in rather large (unphysical) uncertainties. The problem was studied in the scope of the Working Party on Evaluation Cooperation (WPEC) of Nuclear Energy Agency (NEA), Subgroup 26 entitled “Nuclear Data Needs for Advanced Reactor Systems” [32]. Two sensitivity methods were compared, the classical and the “constrained” method. It was finally demonstrated in [21–23] that in the case of the covariance matrix of the constrained (normalized) quantities, which have to obey the zero-sum rule, the constrained sensitivity method is equivalent to correcting (“normalizing”) the matrix. Indeed, if the zero-sum constraint has not been applied in the derivation of the secondary distribution covariance matrix, the covariance matrix is corrected as follows:(6)M~=Sχt·M·Sχ,where the relative transformation matrix defined as Sχ=I-χ⃑·1⃑t with I the unitary matrix, and 1⃑t a transposed unitary vector.However, rather than “correcting” the covariance matrix, the correction can be directly applied to the sensitivity coefficients in the following way:(7)ΔRR2=St·M~·S=St·S⃑χt·M·S⃑χ·S=Sχ·St·M·Sχ·S=Sconst·M·Scons,where S is the unconstrained, Scons the constrained sensitivity method, which automatically assures the normalization of the matrix. An example is shown in Figure 5.Figure 5
Examples of the sensitivities with respect to the PFNS of239Pu calculated by the constrained (/Norm.) and unconstrained sensitivity methods (SNEAK-7 and KRITZ benchmarks).Both constrained and unconstrained sensitivity methods are available in the SUSD3D code and are used, for example, for the sensitivity with respect to the prompt and delayed fission neutron spectra. The constrained and unconstrained sensitivity methods will clearly result in different uncertainties in cases where the fission spectra covariance matrices do not comply exactly with the ENDF-6 Format Manual rules. More information can be found in [21–23].


#### 3.1.3. Beta-Effective Sensitivity and Uncertainty Analysis

The effective delayed neutron fraction is a key reactor safety parameter involved in the control rods worth calculations and transient (reactivity feedback effect) studies. In SUSD3D the sensitivity coefficients ofβeff with respect to the basic nuclear data were calculated by deriving the Bretscher’s k-ratio formula as initially proposed in 2011 (see [28–30]). As demonstrated in these papers (e.g., [29]) the definition of βeff as given in Keepin [33] is equivalent to the 1st order sensitivity of keff to the delayed neutron yield. βeff and all its components by actinides can be therefore obtained using the 1st order perturbation codes such as SUSD3D as a sum of the keff sensitivities with respect to the delayed neutron yields (ν-d) for all fissile isotopes i:(8)βeff=∑i∫Sν-d,iEdE,where: Sν-d,i=νd,i/keff∂keff/∂νd,i.In the Monte Carlo codes the Bretscher’s approximation, sometimes called also the promptk-ratio method, was often used to calculate the beta-effective from(9)βeff=1-kpk,where kp is the keff taking into account only prompt neutrons and k is the total (prompt and delayed neutron) keff.Since the sensitivity ofβeff would require the calculation of the second derivative (βeff being already equal the 1st derivative) another approach was developed in the SUSD3D code [10, 11] based on the derivation of the Bretscher’s prompt k-ratio method:(10)Sβ=σβeff∂βeff∂σ=kpkβeffSk-Skp.The two terms Sk and Skp correspond to the sensitivities of the k and kp which can be obtained using the standard linear perturbation theory. Note also that using the M/C technique βeff sensitivity can relatively easily evaluated as the 2nd derivative of keff.The method implemented in the SUSD3D code was applied to several fast neutron benchmark experiments from the ICSBEP and IRPhE databases [31, 34] and to the MYRRHA reactor.


## 3.1. Special Features Developed in SUSD3D

Several innovative mathematical methods were developed in the scope of the SUSD3D project, and some of them were later introduced also in other S/U computational sequences (e.g., TSUNAMI [24], SERPENT [25], and MCNP [26]) such as the following:(i)
Use of flux angular moments instead of much more voluminous angular fluxes [10, 27](ii)
Demonstration of the equivalence of theconstrained sensitivity method and covariance matrix fix-up in the calculation of the sensitivity to the secondary energy distributions, such as fission spectra sensitivities [21–23](iii)
Sensitivity and uncertainty calculations of theeffective delayed neutron fraction (βeff) [28–30].Future development includes the implementation of the implicit sensitivity method, not yet available in the SUSD3D code and coupling to other transport codes.


### 3.1.1. Flux Moments versus Angular Fluxes: How Sensitive the Sensitivity Calculations Are on PN Order Approximations?

In the mathematical expression of the nuclear data sensitivity function as derived from the first-order perturbation theory [19] both the angular fluxes and their expansion in the Legendre flux moments are required. Both terms are calculated by the traditional discrete ordinates transport codes and stored in the files. This leaves the developers of the nuclear data sensitivity and uncertainty codes several possible choices of the method to account for the neutron flux anisotropy. Three possibilities available in the SUSD3D code were studied here and compared using the FLATTOP-Pu benchmark exercise. This benchmark was particularly suitable due to its high flux anisotropy.Using the first-order perturbation theory, the relative cross-section sensitivity functionSk,x(E), which expresses the fractional change of an integral response, such as the multiplication factor per fractional variation of a cross-section for reaction x in a particular nuclide k, and in energy group g, is expressed in the following form:(1)Sk,xE=σk,xEkeff∂keff∂σk,xE=1R∫dr⃑∫Ω⃑dΩ⃑ρkr⃑-σk,xEΦr⃑,Ω⃑,EΦ∗r⃑,Ω⃑,E︷loss term+∫Ω⃑′dΩ⃑′∫dE′σsk,xΩ⃑⟶Ω⃑′,E⟶E′Φr⃑,Ω⃑,EΦ∗r⃑,Ω⃑′,E′︷gain term,where R is the normalization factor andρkr⃑ = atomic number density of nuclide k,σk,xE = microscopic cross-section for reaction x, in nuclide k, at energy E,σsk,x(Ω⃑→Ω⃑′,E→E′) = microscopic scattering cross-section from angle Ω⃑ to Ω⃑′, and from energy E to E′, for reaction x, in nuclide k,Φr⃑,Ω⃑,E, Φ∗r⃑,Ω⃑,E = direct and adjoint angular fluxes at position r⃑, with direction Ω⃑ and energy E.The angular dependence of the scattering cross-section is usually described in terms of the cosine of the scattering angle and expanded in Legendre polynomials. In practice the Legendre expansion is truncated after finite number of terms, with L representing here the largest term taken into account. Using in addition the multigroup approximation and finite-difference method the above equation (1) becomes(2)Sk,x,g=1R∑iViρk,i-σk,x,g∑mΦg,i,mΦg,i,m∗ΔΩm+∑g′∑l=0L2l+1σsk,x,l,g→g′∑n=-llMg,il,nM∗g′,il,n.In one- and two- dimensional calculations due to symmetry conditions the summation over n goes only from 0 to l, and in three dimensions all the terms from -l to l must be considered. Note that the (4π) term is included in the integration over angles, that is, ∫4πdΩ=1.In the above equation,Vi is the volume of the space mesh interval i, and Φg,i,m and Φg,i,m∗ are the direct and adjoint angular fluxes, corresponding to energy group g, averaged over space interval i and solid angle ΔΩm. Similarly, ρk,i,σk,x,g,σsk,x,l,g→g′ are space and group averaged quantities defined above. Mg,il,n and M∗g′,il,n are direct and adjoint angular flux moments with indices defined above.We see immediately that both angular fluxes and angular flux moments appear in (2). Although in principle the deterministic transport codes usually calculate and provide on output both quantities, the use of both the angular fluxes and angular flux moments is unpractical and unnecessary, since this requires large computer storage. Indeed, the two quantities are linked by the following expressions: (3)Φg,i,m=∑l=0L∑n=-llYl,nΩm·Mg,il,n,Φg,i,m∗=∑l=0L∑n=-llYl,n+Ωm·M∗g,il,n,(4)Mg,il,n=∑mYl,nΩm·Φg,i,m·ΔΩm′,M∗g,il,n=∑mYl,n+Ωm·Φg,i,m∗·ΔΩm.Several alternative methods are consequently available for the computation of the sensitivity profiles using the above equations (1) to (4), namely,(1)
The loss and gain terms of (3) are expressed in terms of the angular flux moments calculated by the discrete ordinates (SN) transport code; that is, the loss term is calculated by (4):(5)Sk,x,g=1R∑iViρk,i∑l=0L2l+1-σk,x,g∑n=-llMg,il,nM∗g′,il,n+∑g′σsk,x,l,g→g′∑n=-llMg,il,nM∗g′,il,n.Note that thePN expansion order is the same for the loss and gain terms; that is, the one used in the transport calculations. Equation (5) is an approximation of (2), with the two being equivalent if high Legendre order expansion is used for the loss term.(2)
The loss term in (2) is expressed in terms of the angular flux calculated using the SN code and the flux moments needed for the gain term are computed using (4) from these angular fluxes.(3)
Angular flux and flux moments from theSN code output are used to calculate the loss and gain terms of the sensitivity by (3).Most perturbation codes make use of the angular fluxes output from the discrete ordinates transport codes and the angular moments are expressed in terms of these angular fluxes (Method (2)). However, the angular flux files may be in some cases quite bulky. On the other hand, the flux moments files are much more reasonably sized, as shown in Table 1 which compares the number of terms (for each energy group and space interval) required in case of using the angular fluxes and flux moments. Using angular moments (the above method (1)) becomes particularly attractive in three dimensions where, for example, a typically used S16-P3 approximation requires only 16 flux moment terms, but as many as 288 angular flux terms (note that this number is to be multiplied by the number of energy groups and space points to obtain the total number of flux terms needed in a sensitivity calculation).Table 1
Number of flux moment terms used in PARTISN for different Legendre expansion orderL and the number of discrete directions for different SN fully symmetrical direction sets.
Number of flux moments
Number of discrete directions
P0
P1
P3
P3
P3
P5
S2
S4
S8
S16
1D
1
2
4
4
4
6
2
4
8
16
1D-(r, z)
1
2
6
6
6
12
2
6
20
72
2D
1
3
10
10
10
21
4
12
40
144
3D
1
4
16
16
16
26
8
24
80
288An additional advantage of using flux moments, besides a considerable reduction of computer space requirements and CPU time gain, resides in the fact that manySN codes (e.g., ANISN [1]) print the angular fluxes at the space cell edges, and sophisticated negative flux fix-up procedures are subsequently adopted in the transport codes to construct the flux at the center of the cell. Flux moments are conveniently calculated at the center of space mesh.In case of using the flux moments it is on the other hand expected that higherPN Legendre expansion order may be needed as compared to using the angular fluxes, since also the angular fluxes for the loss term are expressed in terms of the angular moments.The SUSD3D code offers the possibility of using the above described options. The mostly used approach is based on (5) (method (1)) calculating the relative sensitivity profiles Sk,x,g from the reasonably sized angular flux moments, produced by the SN codes DORT-TORT, DANTSYS, and PARTISN. Additionally, for one- and two-dimensional analysis, the angular flux files produced by the ANISN [1] and DOT-3.5 [20] codes can be used, thus providing a convenient way to compare and test the results using the two methods.The performance of the above described methods was tested on a highly anisotropic problem while using a relatively low Legendre order. The FLATTOP-Pu benchmark was found to be very sensitive to the scattering anisotropy and as such very suitable for studying the advantages and inconveniences of different approaches. The benchmark, sometimes also referred to as Popsy, consists of ~20-cm natural U reflected 94 wt.%Pu239 sphere, with a radius of 4.533-cm. It is described in the ICSBEP criticality benchmark handbook [31] under the reference PU-MET-FAST-006.FLATTOP-Pu was modeled in a simple 1D spherical geometry using the PARTISN and ANISN transport and SUSD3D S/U codes. S16 fully symmetrical quadratures, ENDF/B-VII.0 33-energy group cross-sections, and the covariance matrices from SCALE-6.1 [15] were used. Legendre expansion orders were varied (P5, P3, and P1). Bell-Hansen-Sadmeister (BHS) transport correction was used for the P1 and P3 calculations.The comparison of the calculated multiplication factors is presented in Table2 showing good consistency between the PARTISN and ANISN, as well as between direct and adjoint results. Small differences may be due to the use of a single fission spectra in ANISN contrary to the zone-dependent one used in PARTISN. BHS correction is shown to improve considerably the results (note ~1700 pcm difference with respect to the uncorrected P1 calculations). P3 results are reasonably close to those obtained using P5.Table 2
Multiplication factors calculated using the PARTISN and ANISN codes using P1 and P5 Legendre expansion. BHS indicates that the Bell-Hansen-Sadmeister transport correction was used. Note an excellent agreement between the direct and adjoint calculations, within few pcm for PARTISN, and within ~15 pcm for ANISN calculations.
Code/options
k
e
f
f-direct
k
e
f
f-adjoint
PARTISN/P5
1.00030
1.00028
PARTISN/P3
1.00013
1.00011
PARTISN/P1
0.99902
0.99900
ANISN/P5
1.00337
1.00320
ANISN/P1
1.00208
1.00193
ANISN/P1 (no BHS correction)
0.98463
0.98448As mentioned, SUSD3D uses angular moment files produced by PARTISN to resolve (3) (Method (1)), or alternatively the angular flux files from ANISN to resolve (2) by Method (2). In case of using the angular fluxes from ANISN, the flux moments are first calculated from the angular fluxes. However, PN order in principle does not need to be the same in the ANISN transport and in the SUSD3D sensitivity calculations. Note that using a high PN order in the loss and the same PN order as in the transport calculation for the gain term is approximately equivalent to Method (3), that is, the exact formulation of the sensitivity as defined in (2).The following computational procedures were intercompared and the corresponding results are presented in Table3 and Figures 3 and 4:(i)
PARTISN/SUSD3D with P5 approximation(ii)
PARTISN using P3 with BHS cross-section correction/SUSD3D using P3(iii)
PARTISN using P1 with BHS cross-section correction/SUSD3D using P1(iv)
ANISN/SUSD3D both using P5(v)
ANISN using P5 and SUSD3D using P1(vi)
ANISN using P1 with BHS cross-section correction/SUSD3D using P1(vii)
ANISN using P1 with BHS cross-section correction/SUSD3D using P5(viii)
ANISN using P1 with BHS correction/SUSD3D using P5 for loss & P1 for gain.The choice of the PN expansion has little impact on the sensitivities of the keff factor to the absorption, fission, and nu-bar data. As expected, the keff factor was found to be particularly sensitive to 238U elastic and inelastic scattering cross-sections and these reactions are also very sensitive to the PN anisotropic scattering approximations; therefore only these results are presented in Table 3 and Figures 3 and 4.Table 3
Comparison ofkeff sensitivity and uncertainty results using different approaches.
Code/options
Sensitivity (%/%)
Uncertainty (%)
238U el
238U inel
238U(el + inel)
Total
PARTISN/P5
0.1399
0.06451
0.406%
1.20%
PARTISN P3
0.1431
0.06776
0.426%
1.21%
PARTISN P1
0.1721
0.09731
0.681%
1.32%
ANISN-P5/SUSD3D-P5
0.1395
0.06325
0.388%
1.19%
ANISN-P5/SUSD3D-P1
0.1710
0.09543
0.654%
1.30%
ANISN P1/SUSD3D-P1
0.1712
0.09550
0.656%
1.31%
ANISN P1/SUSD3D-P5
0.1407
0.064445
0.397%
1.19%
ANISN P1/SUSD3D-P1_gain
0.1180
0.06466
0.633%
1.29%Figure 3
Sensitivity ofkeff to U238 elastic cross-sections calculated by SUSD3D/ANISN codes.P5-P5 = P5 transport & sensitivity calculation,P1-P1 = P1 transport using BHS cross-sections and P1 sensitivity,P1-P5 = P1 transport using BHS, P5 sensitivity calculations,P1-P1 gain = P1 transport using BHS cross-sections, and sensitivity calculation using (2) with P1 gain term.Figure 4
Sensitivity ofkeff to 238U inelastic cross-sections calculated by SUSD3D/ANISN codes.It is important to note however that the BHS correction was not (yet) applied to the cross-sections used in the SUSD3D sensitivity calculations making the two calculations at this stage slightly inconsistent.Several conclusions can be drawn from the results presented in Table3 and on Figures 3 and 4:(i)
Good consistency is observed between the SUSD3D results based on PARTISN angular moments and ANISN angular fluxes, both for P5 and P1 results.(ii)
The analysis of the FLATTOP-Pu critical benchmark revealed large differences in calculatedkeff and the corresponding sensitivity profiles when using low Legendre order expansion such as P1 as compared to the reference P5 (and P3) calculation. The use of BHS correction results in a good agreement with keff factors calculated using the P1 BHS and P5 expansions. It is expected that the relatively large differences between the sensitivities based on P5 and P1 BHS transport calculations could be resolved to large extent by applying the BHS correction to the cross-sections used in the SUSD3D sensitivity code.(iii)
Surprisingly, the use of highPN order for the loss and P1 order for the gain term (SUSD3D-P1_gain), which is equivalent to the “correct” (2) (Method (3)), does not perform necessarily much better than Method (1) (SUSD3D-P1) routinely used in SUSD3D. It can be seen to underpredict the sensitivities to the 238U elastic (Figure 3) and slightly overpredict those to the 238U inelastic (Figure 4) cross-sections. Taking into account the computer space penalties linked to this method and the limited improvements in the results, it seems that the use of (2) (Method (3)) is therefore not justified. Lower PN order calculations would be possible using the angular fluxes and BHS corrected cross-sections, but again at high computer space requirements.


### 3.1.2. Constrained Sensitivities: Demonstration of the Equivalence of the Constrained Sensitivity Method and the Covariance Matrix Fix-Up

The SUSD3D analysis revealed several differences between the results based on different prompt fission neutron spectra (PFNS) covariances, in some cases resulting in rather large (unphysical) uncertainties. The problem was studied in the scope of the Working Party on Evaluation Cooperation (WPEC) of Nuclear Energy Agency (NEA), Subgroup 26 entitled “Nuclear Data Needs for Advanced Reactor Systems” [32]. Two sensitivity methods were compared, the classical and the “constrained” method. It was finally demonstrated in [21–23] that in the case of the covariance matrix of the constrained (normalized) quantities, which have to obey the zero-sum rule, the constrained sensitivity method is equivalent to correcting (“normalizing”) the matrix. Indeed, if the zero-sum constraint has not been applied in the derivation of the secondary distribution covariance matrix, the covariance matrix is corrected as follows:(6)M~=Sχt·M·Sχ,where the relative transformation matrix defined as Sχ=I-χ⃑·1⃑t with I the unitary matrix, and 1⃑t a transposed unitary vector.However, rather than “correcting” the covariance matrix, the correction can be directly applied to the sensitivity coefficients in the following way:(7)ΔRR2=St·M~·S=St·S⃑χt·M·S⃑χ·S=Sχ·St·M·Sχ·S=Sconst·M·Scons,where S is the unconstrained, Scons the constrained sensitivity method, which automatically assures the normalization of the matrix. An example is shown in Figure 5.Figure 5
Examples of the sensitivities with respect to the PFNS of239Pu calculated by the constrained (/Norm.) and unconstrained sensitivity methods (SNEAK-7 and KRITZ benchmarks).Both constrained and unconstrained sensitivity methods are available in the SUSD3D code and are used, for example, for the sensitivity with respect to the prompt and delayed fission neutron spectra. The constrained and unconstrained sensitivity methods will clearly result in different uncertainties in cases where the fission spectra covariance matrices do not comply exactly with the ENDF-6 Format Manual rules. More information can be found in [21–23].


### 3.1.3. Beta-Effective Sensitivity and Uncertainty Analysis

The effective delayed neutron fraction is a key reactor safety parameter involved in the control rods worth calculations and transient (reactivity feedback effect) studies. In SUSD3D the sensitivity coefficients ofβeff with respect to the basic nuclear data were calculated by deriving the Bretscher’s k-ratio formula as initially proposed in 2011 (see [28–30]). As demonstrated in these papers (e.g., [29]) the definition of βeff as given in Keepin [33] is equivalent to the 1st order sensitivity of keff to the delayed neutron yield. βeff and all its components by actinides can be therefore obtained using the 1st order perturbation codes such as SUSD3D as a sum of the keff sensitivities with respect to the delayed neutron yields (ν-d) for all fissile isotopes i:(8)βeff=∑i∫Sν-d,iEdE,where: Sν-d,i=νd,i/keff∂keff/∂νd,i.In the Monte Carlo codes the Bretscher’s approximation, sometimes called also the promptk-ratio method, was often used to calculate the beta-effective from(9)βeff=1-kpk,where kp is the keff taking into account only prompt neutrons and k is the total (prompt and delayed neutron) keff.Since the sensitivity ofβeff would require the calculation of the second derivative (βeff being already equal the 1st derivative) another approach was developed in the SUSD3D code [10, 11] based on the derivation of the Bretscher’s prompt k-ratio method:(10)Sβ=σβeff∂βeff∂σ=kpkβeffSk-Skp.The two terms Sk and Skp correspond to the sensitivities of the k and kp which can be obtained using the standard linear perturbation theory. Note also that using the M/C technique βeff sensitivity can relatively easily evaluated as the 2nd derivative of keff.The method implemented in the SUSD3D code was applied to several fast neutron benchmark experiments from the ICSBEP and IRPhE databases [31, 34] and to the MYRRHA reactor.


## 3.1.1. Flux Moments versus Angular Fluxes: How Sensitive the Sensitivity Calculations Are on PN Order Approximations?

In the mathematical expression of the nuclear data sensitivity function as derived from the first-order perturbation theory [19] both the angular fluxes and their expansion in the Legendre flux moments are required. Both terms are calculated by the traditional discrete ordinates transport codes and stored in the files. This leaves the developers of the nuclear data sensitivity and uncertainty codes several possible choices of the method to account for the neutron flux anisotropy. Three possibilities available in the SUSD3D code were studied here and compared using the FLATTOP-Pu benchmark exercise. This benchmark was particularly suitable due to its high flux anisotropy.Using the first-order perturbation theory, the relative cross-section sensitivity functionSk,x(E), which expresses the fractional change of an integral response, such as the multiplication factor per fractional variation of a cross-section for reaction x in a particular nuclide k, and in energy group g, is expressed in the following form:(1)Sk,xE=σk,xEkeff∂keff∂σk,xE=1R∫dr⃑∫Ω⃑dΩ⃑ρkr⃑-σk,xEΦr⃑,Ω⃑,EΦ∗r⃑,Ω⃑,E︷loss term+∫Ω⃑′dΩ⃑′∫dE′σsk,xΩ⃑⟶Ω⃑′,E⟶E′Φr⃑,Ω⃑,EΦ∗r⃑,Ω⃑′,E′︷gain term,where R is the normalization factor andρkr⃑ = atomic number density of nuclide k,σk,xE = microscopic cross-section for reaction x, in nuclide k, at energy E,σsk,x(Ω⃑→Ω⃑′,E→E′) = microscopic scattering cross-section from angle Ω⃑ to Ω⃑′, and from energy E to E′, for reaction x, in nuclide k,Φr⃑,Ω⃑,E, Φ∗r⃑,Ω⃑,E = direct and adjoint angular fluxes at position r⃑, with direction Ω⃑ and energy E.The angular dependence of the scattering cross-section is usually described in terms of the cosine of the scattering angle and expanded in Legendre polynomials. In practice the Legendre expansion is truncated after finite number of terms, with L representing here the largest term taken into account. Using in addition the multigroup approximation and finite-difference method the above equation (1) becomes(2)Sk,x,g=1R∑iViρk,i-σk,x,g∑mΦg,i,mΦg,i,m∗ΔΩm+∑g′∑l=0L2l+1σsk,x,l,g→g′∑n=-llMg,il,nM∗g′,il,n.In one- and two- dimensional calculations due to symmetry conditions the summation over n goes only from 0 to l, and in three dimensions all the terms from -l to l must be considered. Note that the (4π) term is included in the integration over angles, that is, ∫4πdΩ=1.In the above equation,Vi is the volume of the space mesh interval i, and Φg,i,m and Φg,i,m∗ are the direct and adjoint angular fluxes, corresponding to energy group g, averaged over space interval i and solid angle ΔΩm. Similarly, ρk,i,σk,x,g,σsk,x,l,g→g′ are space and group averaged quantities defined above. Mg,il,n and M∗g′,il,n are direct and adjoint angular flux moments with indices defined above.We see immediately that both angular fluxes and angular flux moments appear in (2). Although in principle the deterministic transport codes usually calculate and provide on output both quantities, the use of both the angular fluxes and angular flux moments is unpractical and unnecessary, since this requires large computer storage. Indeed, the two quantities are linked by the following expressions: (3)Φg,i,m=∑l=0L∑n=-llYl,nΩm·Mg,il,n,Φg,i,m∗=∑l=0L∑n=-llYl,n+Ωm·M∗g,il,n,(4)Mg,il,n=∑mYl,nΩm·Φg,i,m·ΔΩm′,M∗g,il,n=∑mYl,n+Ωm·Φg,i,m∗·ΔΩm.Several alternative methods are consequently available for the computation of the sensitivity profiles using the above equations (1) to (4), namely,(1)
The loss and gain terms of (3) are expressed in terms of the angular flux moments calculated by the discrete ordinates (SN) transport code; that is, the loss term is calculated by (4):(5)Sk,x,g=1R∑iViρk,i∑l=0L2l+1-σk,x,g∑n=-llMg,il,nM∗g′,il,n+∑g′σsk,x,l,g→g′∑n=-llMg,il,nM∗g′,il,n.Note that thePN expansion order is the same for the loss and gain terms; that is, the one used in the transport calculations. Equation (5) is an approximation of (2), with the two being equivalent if high Legendre order expansion is used for the loss term.(2)
The loss term in (2) is expressed in terms of the angular flux calculated using the SN code and the flux moments needed for the gain term are computed using (4) from these angular fluxes.(3)
Angular flux and flux moments from theSN code output are used to calculate the loss and gain terms of the sensitivity by (3).Most perturbation codes make use of the angular fluxes output from the discrete ordinates transport codes and the angular moments are expressed in terms of these angular fluxes (Method (2)). However, the angular flux files may be in some cases quite bulky. On the other hand, the flux moments files are much more reasonably sized, as shown in Table 1 which compares the number of terms (for each energy group and space interval) required in case of using the angular fluxes and flux moments. Using angular moments (the above method (1)) becomes particularly attractive in three dimensions where, for example, a typically used S16-P3 approximation requires only 16 flux moment terms, but as many as 288 angular flux terms (note that this number is to be multiplied by the number of energy groups and space points to obtain the total number of flux terms needed in a sensitivity calculation).Table 1
Number of flux moment terms used in PARTISN for different Legendre expansion orderL and the number of discrete directions for different SN fully symmetrical direction sets.
Number of flux moments
Number of discrete directions
P0
P1
P3
P3
P3
P5
S2
S4
S8
S16
1D
1
2
4
4
4
6
2
4
8
16
1D-(r, z)
1
2
6
6
6
12
2
6
20
72
2D
1
3
10
10
10
21
4
12
40
144
3D
1
4
16
16
16
26
8
24
80
288An additional advantage of using flux moments, besides a considerable reduction of computer space requirements and CPU time gain, resides in the fact that manySN codes (e.g., ANISN [1]) print the angular fluxes at the space cell edges, and sophisticated negative flux fix-up procedures are subsequently adopted in the transport codes to construct the flux at the center of the cell. Flux moments are conveniently calculated at the center of space mesh.In case of using the flux moments it is on the other hand expected that higherPN Legendre expansion order may be needed as compared to using the angular fluxes, since also the angular fluxes for the loss term are expressed in terms of the angular moments.The SUSD3D code offers the possibility of using the above described options. The mostly used approach is based on (5) (method (1)) calculating the relative sensitivity profiles Sk,x,g from the reasonably sized angular flux moments, produced by the SN codes DORT-TORT, DANTSYS, and PARTISN. Additionally, for one- and two-dimensional analysis, the angular flux files produced by the ANISN [1] and DOT-3.5 [20] codes can be used, thus providing a convenient way to compare and test the results using the two methods.The performance of the above described methods was tested on a highly anisotropic problem while using a relatively low Legendre order. The FLATTOP-Pu benchmark was found to be very sensitive to the scattering anisotropy and as such very suitable for studying the advantages and inconveniences of different approaches. The benchmark, sometimes also referred to as Popsy, consists of ~20-cm natural U reflected 94 wt.%Pu239 sphere, with a radius of 4.533-cm. It is described in the ICSBEP criticality benchmark handbook [31] under the reference PU-MET-FAST-006.FLATTOP-Pu was modeled in a simple 1D spherical geometry using the PARTISN and ANISN transport and SUSD3D S/U codes. S16 fully symmetrical quadratures, ENDF/B-VII.0 33-energy group cross-sections, and the covariance matrices from SCALE-6.1 [15] were used. Legendre expansion orders were varied (P5, P3, and P1). Bell-Hansen-Sadmeister (BHS) transport correction was used for the P1 and P3 calculations.The comparison of the calculated multiplication factors is presented in Table2 showing good consistency between the PARTISN and ANISN, as well as between direct and adjoint results. Small differences may be due to the use of a single fission spectra in ANISN contrary to the zone-dependent one used in PARTISN. BHS correction is shown to improve considerably the results (note ~1700 pcm difference with respect to the uncorrected P1 calculations). P3 results are reasonably close to those obtained using P5.Table 2
Multiplication factors calculated using the PARTISN and ANISN codes using P1 and P5 Legendre expansion. BHS indicates that the Bell-Hansen-Sadmeister transport correction was used. Note an excellent agreement between the direct and adjoint calculations, within few pcm for PARTISN, and within ~15 pcm for ANISN calculations.
Code/options
k
e
f
f-direct
k
e
f
f-adjoint
PARTISN/P5
1.00030
1.00028
PARTISN/P3
1.00013
1.00011
PARTISN/P1
0.99902
0.99900
ANISN/P5
1.00337
1.00320
ANISN/P1
1.00208
1.00193
ANISN/P1 (no BHS correction)
0.98463
0.98448As mentioned, SUSD3D uses angular moment files produced by PARTISN to resolve (3) (Method (1)), or alternatively the angular flux files from ANISN to resolve (2) by Method (2). In case of using the angular fluxes from ANISN, the flux moments are first calculated from the angular fluxes. However, PN order in principle does not need to be the same in the ANISN transport and in the SUSD3D sensitivity calculations. Note that using a high PN order in the loss and the same PN order as in the transport calculation for the gain term is approximately equivalent to Method (3), that is, the exact formulation of the sensitivity as defined in (2).The following computational procedures were intercompared and the corresponding results are presented in Table3 and Figures 3 and 4:(i)
PARTISN/SUSD3D with P5 approximation(ii)
PARTISN using P3 with BHS cross-section correction/SUSD3D using P3(iii)
PARTISN using P1 with BHS cross-section correction/SUSD3D using P1(iv)
ANISN/SUSD3D both using P5(v)
ANISN using P5 and SUSD3D using P1(vi)
ANISN using P1 with BHS cross-section correction/SUSD3D using P1(vii)
ANISN using P1 with BHS cross-section correction/SUSD3D using P5(viii)
ANISN using P1 with BHS correction/SUSD3D using P5 for loss & P1 for gain.The choice of the PN expansion has little impact on the sensitivities of the keff factor to the absorption, fission, and nu-bar data. As expected, the keff factor was found to be particularly sensitive to 238U elastic and inelastic scattering cross-sections and these reactions are also very sensitive to the PN anisotropic scattering approximations; therefore only these results are presented in Table 3 and Figures 3 and 4.Table 3
Comparison ofkeff sensitivity and uncertainty results using different approaches.
Code/options
Sensitivity (%/%)
Uncertainty (%)
238U el
238U inel
238U(el + inel)
Total
PARTISN/P5
0.1399
0.06451
0.406%
1.20%
PARTISN P3
0.1431
0.06776
0.426%
1.21%
PARTISN P1
0.1721
0.09731
0.681%
1.32%
ANISN-P5/SUSD3D-P5
0.1395
0.06325
0.388%
1.19%
ANISN-P5/SUSD3D-P1
0.1710
0.09543
0.654%
1.30%
ANISN P1/SUSD3D-P1
0.1712
0.09550
0.656%
1.31%
ANISN P1/SUSD3D-P5
0.1407
0.064445
0.397%
1.19%
ANISN P1/SUSD3D-P1_gain
0.1180
0.06466
0.633%
1.29%Figure 3
Sensitivity ofkeff to U238 elastic cross-sections calculated by SUSD3D/ANISN codes.P5-P5 = P5 transport & sensitivity calculation,P1-P1 = P1 transport using BHS cross-sections and P1 sensitivity,P1-P5 = P1 transport using BHS, P5 sensitivity calculations,P1-P1 gain = P1 transport using BHS cross-sections, and sensitivity calculation using (2) with P1 gain term.Figure 4
Sensitivity ofkeff to 238U inelastic cross-sections calculated by SUSD3D/ANISN codes.It is important to note however that the BHS correction was not (yet) applied to the cross-sections used in the SUSD3D sensitivity calculations making the two calculations at this stage slightly inconsistent.Several conclusions can be drawn from the results presented in Table3 and on Figures 3 and 4:(i)
Good consistency is observed between the SUSD3D results based on PARTISN angular moments and ANISN angular fluxes, both for P5 and P1 results.(ii)
The analysis of the FLATTOP-Pu critical benchmark revealed large differences in calculatedkeff and the corresponding sensitivity profiles when using low Legendre order expansion such as P1 as compared to the reference P5 (and P3) calculation. The use of BHS correction results in a good agreement with keff factors calculated using the P1 BHS and P5 expansions. It is expected that the relatively large differences between the sensitivities based on P5 and P1 BHS transport calculations could be resolved to large extent by applying the BHS correction to the cross-sections used in the SUSD3D sensitivity code.(iii)
Surprisingly, the use of highPN order for the loss and P1 order for the gain term (SUSD3D-P1_gain), which is equivalent to the “correct” (2) (Method (3)), does not perform necessarily much better than Method (1) (SUSD3D-P1) routinely used in SUSD3D. It can be seen to underpredict the sensitivities to the 238U elastic (Figure 3) and slightly overpredict those to the 238U inelastic (Figure 4) cross-sections. Taking into account the computer space penalties linked to this method and the limited improvements in the results, it seems that the use of (2) (Method (3)) is therefore not justified. Lower PN order calculations would be possible using the angular fluxes and BHS corrected cross-sections, but again at high computer space requirements.


## 3.1.2. Constrained Sensitivities: Demonstration of the Equivalence of the Constrained Sensitivity Method and the Covariance Matrix Fix-Up

The SUSD3D analysis revealed several differences between the results based on different prompt fission neutron spectra (PFNS) covariances, in some cases resulting in rather large (unphysical) uncertainties. The problem was studied in the scope of the Working Party on Evaluation Cooperation (WPEC) of Nuclear Energy Agency (NEA), Subgroup 26 entitled “Nuclear Data Needs for Advanced Reactor Systems” [32]. Two sensitivity methods were compared, the classical and the “constrained” method. It was finally demonstrated in [21–23] that in the case of the covariance matrix of the constrained (normalized) quantities, which have to obey the zero-sum rule, the constrained sensitivity method is equivalent to correcting (“normalizing”) the matrix. Indeed, if the zero-sum constraint has not been applied in the derivation of the secondary distribution covariance matrix, the covariance matrix is corrected as follows:(6)M~=Sχt·M·Sχ,where the relative transformation matrix defined as Sχ=I-χ⃑·1⃑t with I the unitary matrix, and 1⃑t a transposed unitary vector.However, rather than “correcting” the covariance matrix, the correction can be directly applied to the sensitivity coefficients in the following way:(7)ΔRR2=St·M~·S=St·S⃑χt·M·S⃑χ·S=Sχ·St·M·Sχ·S=Sconst·M·Scons,where S is the unconstrained, Scons the constrained sensitivity method, which automatically assures the normalization of the matrix. An example is shown in Figure 5.Figure 5
Examples of the sensitivities with respect to the PFNS of239Pu calculated by the constrained (/Norm.) and unconstrained sensitivity methods (SNEAK-7 and KRITZ benchmarks).Both constrained and unconstrained sensitivity methods are available in the SUSD3D code and are used, for example, for the sensitivity with respect to the prompt and delayed fission neutron spectra. The constrained and unconstrained sensitivity methods will clearly result in different uncertainties in cases where the fission spectra covariance matrices do not comply exactly with the ENDF-6 Format Manual rules. More information can be found in [21–23].


## 3.1.3. Beta-Effective Sensitivity and Uncertainty Analysis

The effective delayed neutron fraction is a key reactor safety parameter involved in the control rods worth calculations and transient (reactivity feedback effect) studies. In SUSD3D the sensitivity coefficients ofβeff with respect to the basic nuclear data were calculated by deriving the Bretscher’s k-ratio formula as initially proposed in 2011 (see [28–30]). As demonstrated in these papers (e.g., [29]) the definition of βeff as given in Keepin [33] is equivalent to the 1st order sensitivity of keff to the delayed neutron yield. βeff and all its components by actinides can be therefore obtained using the 1st order perturbation codes such as SUSD3D as a sum of the keff sensitivities with respect to the delayed neutron yields (ν-d) for all fissile isotopes i:(8)βeff=∑i∫Sν-d,iEdE,where: Sν-d,i=νd,i/keff∂keff/∂νd,i.In the Monte Carlo codes the Bretscher’s approximation, sometimes called also the promptk-ratio method, was often used to calculate the beta-effective from(9)βeff=1-kpk,where kp is the keff taking into account only prompt neutrons and k is the total (prompt and delayed neutron) keff.Since the sensitivity ofβeff would require the calculation of the second derivative (βeff being already equal the 1st derivative) another approach was developed in the SUSD3D code [10, 11] based on the derivation of the Bretscher’s prompt k-ratio method:(10)Sβ=σβeff∂βeff∂σ=kpkβeffSk-Skp.The two terms Sk and Skp correspond to the sensitivities of the k and kp which can be obtained using the standard linear perturbation theory. Note also that using the M/C technique βeff sensitivity can relatively easily evaluated as the 2nd derivative of keff.The method implemented in the SUSD3D code was applied to several fast neutron benchmark experiments from the ICSBEP and IRPhE databases [31, 34] and to the MYRRHA reactor.


## 4. Use and Validation of XSUN-2017

The SUSD3D code was extensively used and tested in the past ~25 years within several projects such as European fusion program, PWR pressure vessel surveillance programme, and medical and industrial projects. In the scope of the European fusion programme several benchmark experiments performed at the FNG facility in Frascati were analyzed [35, 36]. SUSD3D was used to validate the Monte Carlo cross-section sensitivity-uncertainty code MCSEN, including the secondary angular distribution sensitivity options. Recently, the performance of the XSUN-2017 code system was compared with other S/U computer codes in the scope of the following:(i)
S/U intercomparison exercise using the SNEAK-7A and -7B benchmark experiments [37, 38] from the International Reactor Physics Benchmark Experiments (IRPhE) database. The performance of the XSUN-2017 code system was compared with the codes such as TSUNAMI-3D, XSUSA, SERPENT2 GPT, MCNP6, and SNATCH.(ii)
Transport, sensitivity, and uncertainty analysis of thekeff and βeff parameters for the MYRRHA accelerator driven system (ADS) [39].


### 4.1. Sensitivity Code Intercomparison Using SNEAK-7 Benchmarks

A sensitivity benchmark exercise was organized within the scope of the Uncertainty Analysis in Modeling (UAM) project [40, 41] of the OECD/Nuclear Energy Agency (NEA) to compare the available and stimulate the development of new, deterministic, and stochastic methods and codes for the sensitivity and uncertainty computations of the effective multiplication factor (keff) and the effective delayed neutron fraction (βeff). Several solutions were received from participants using both deterministic (SUSD3D, SNATCH) and Monte Carlo codes (TSUNAMI-3D, XSUSA, SERPENT2 GPT, and MCNP6). The sensitivity and uncertainty codes were applied to the SNEAK-7A and -7B [42] fast neutron benchmark experiments from the IRPhE database and the results of the keff and βeff S/U intercomparisons are described in [37] and [38], respectively. The SNEAK facility is a fixed vertical assembly with fuel elements suspended from a grid plate. The fuel element (lattice) pitch is 5.44 cm. Platelets of various thicknesses are stacked horizontally within square fuel element tubes. The cross-section of the platelets is 5.077 × 5.077 cm2. In SNEAK-7A, the core unit cell consists of one PuO2-UO2 platelet (26.6% PuO2 and 73.4% UO2) and one graphite platelet. Radial and axial blankets are loaded with depleted UO2 plates. In SNEAK-7B the graphite platelet in the cell of 7A is replaced by UnatO2 platelet resulting in an average Pu-enrichment of about 13%.The criticality of SNEAK-7A and 7B was determined with all control rods in their most reactive position, that is, with the fueled portion of the rod in the core. The measurements were performed on critical eigenvalue (keff), and material buckling, reaction rate ratios, material worth, fission rate, and capture rate distributions as well as effective delayed neutron fraction βeff were measured. The description of the benchmark is available in the IRPhE handbook as the SNEAK-LMFR-EXP-001 evaluation. Two- and three-dimensional models for the TWODANT/THREEDANT (similar to PARTISN) and MCNP codes are included in the IRPhE evaluation. The 3D models of the benchmark models are shown on Figure 6.Figure 6
Cross cuts of detailed and simplified 3D benchmark models of the SNEAK 7A (a) and SNEAK-7B (b) fast reactor assemblies.
(a)
(b)ENDF/B-VII cross-sections (versions .0 or .1) were used in all the calculations, but in very different energy representations, that is, as multigroup with 33 to 238 groups (in deterministic codes) and pointwise cross-sections (M/C codes). Some differences were also in the computational models developed by different participants. Nevertheless, an excellent agreement between different codes was observed, both for integral values and sensitivity profiles. Thekeff sensitivities are compared in Figure 7. For most of the important reactions, the integral sensitivities agree within 1-2%. Slightly larger differences of up to around 10% were observed for the elastic and inelastic scattering, traditionally difficult since evaluated as a difference between two large (gain and loss) components causing numerical difficulties. Note also an excellent match of the sensitivities to the secondary angular distributions (SAD) calculated using SUSD3D and SERPENT codes for P1 to P3 terms of the 238U elastic scattering.Figure 7
Example of sensitivities ofkeff to 235,238U, 239Pu cross-sections, prompt fission spectra (chi), and P1 to P3 elastic scattering angular distributions (SAD) calculated by the SUSD3D, SERPENT, TSUNAMI3D, MCNP6, and SNATCH codes, including numerically demanding elastic and inelastic data sensitivities.For theβeff S/U calculations the solutions obtained using the codes SUSD3D, SERPENT, and XSUSA are compared in Figure 8. Differences are slightly larger than in the case of keff due to numerically demanding calculations of (10) and M/C statistical errors, but still very satisfactory for the validation of the βeff S/U methods. More details on the comparison of different codes and methods are available in [19, 20].Figure 8
Sensitivity ofβeff to 238U inelastic and 238U, Pu239,240 fission cross-sections, calculated using the SUSD3D, XSUSA, and SERPENT codes for the SNEAK-7A benchmark.


### 4.2. MYRRHA ADS Reactor S/U Analysis

MYRRHA (Multipurpose hYbrid Research Reactor for High-tech Applications) is being designed at SCK•CEN, Mol, Belgium, conceived to operate both in subcritical (Accelerator Driven System-ADS) and in critical modes. In the scope of the CHANDA project of the European Commission a sensitivity and uncertainty analysis of the neutron multiplication factorkeff and the effective delayed neutron fraction βeff was performed to identify the most important nuclear data for neutron-induced reactions for criticality calculations of the latest MYRRHA designs. The sensitivity profiles were calculated using different, both stochastic and deterministic, codes (i.e., SCALE/TSUNAMI, MCNP, and XSUN/SUSD3D). MCNP and TSUNAMI calculations were performed using a detailed 3D geometry model. A simplified 2D r-z geometrical model was prepared for SUSD3D (and TSUNAMI) based on the detailed MCNP model of MYRRHA. The problem was still relatively large, with the dimensions of 5 m (r-direction) and 8.5 m (z-direction), 24 material mixtures with 1257 materials at different temperatures ranging from 300 K to 800 K. 33-group cross-sections were used in PARTISN and SUSD3D calculations, pointwise data in MCNP, and 238 energy group library in SCALE/TSUNAMI codes.The XSUN-2017 sequence with the geometry is shown on Figure9. The neutron flux distributions in the MYRRHA reactor generated using the MeshTal Viewer [17] utility called from XSUN-2017 are shown in Figures 10–12.Figure 9
XSUN-2017 computational running segment with the view of the MYRRHA geometry modeling.Figure 10
Distribution of the thermal neutron flux in MYRRHA calculated using the PARTISN code.Figure 11
Distribution of the 1-2 keV neutron flux in MYRRHA calculated using the PARTISN code.Figure 12
Distribution of ~3 MeV neutron flux in MYRRHA calculated using the PARTISN code.In spite of the differences in S/U methodologies, different cross-section energy representations, and cross-section processing and treatment, the neutron-induced nuclear data sensitivity analysis resulted in differences less than 4% between codes, with few exceptions (239Pu fission, 238U elastic scattering, and 56Fe capture reactions) (Figures 13 and 14). Nuclear data covariance matrices of different libraries (SCALE-6, ENDF/B-VII.1, and JENDL-4.0u) were used to derive the uncertainty in keff based on the calculated sensitivities. This study concluded that, depending on the covariance data used, the keff and βeff uncertainties due to the uncertainty in nuclear data were about 0.5–1% and ~2.2%, respectively.Figure 13
Sensitivity ofkeff calculated using the SUSD3D and TSUNAMI codes for the MYRRHA reactor.Figure 14
Sensitivity ofβeff to the total neutron yield ofU235, 238U, 239Pu, 240Pu and Pu241 calculated using the SUSD3D code for the MYRRHA reactor.


## 4.1. Sensitivity Code Intercomparison Using SNEAK-7 Benchmarks

A sensitivity benchmark exercise was organized within the scope of the Uncertainty Analysis in Modeling (UAM) project [40, 41] of the OECD/Nuclear Energy Agency (NEA) to compare the available and stimulate the development of new, deterministic, and stochastic methods and codes for the sensitivity and uncertainty computations of the effective multiplication factor (keff) and the effective delayed neutron fraction (βeff). Several solutions were received from participants using both deterministic (SUSD3D, SNATCH) and Monte Carlo codes (TSUNAMI-3D, XSUSA, SERPENT2 GPT, and MCNP6). The sensitivity and uncertainty codes were applied to the SNEAK-7A and -7B [42] fast neutron benchmark experiments from the IRPhE database and the results of the keff and βeff S/U intercomparisons are described in [37] and [38], respectively. The SNEAK facility is a fixed vertical assembly with fuel elements suspended from a grid plate. The fuel element (lattice) pitch is 5.44 cm. Platelets of various thicknesses are stacked horizontally within square fuel element tubes. The cross-section of the platelets is 5.077 × 5.077 cm2. In SNEAK-7A, the core unit cell consists of one PuO2-UO2 platelet (26.6% PuO2 and 73.4% UO2) and one graphite platelet. Radial and axial blankets are loaded with depleted UO2 plates. In SNEAK-7B the graphite platelet in the cell of 7A is replaced by UnatO2 platelet resulting in an average Pu-enrichment of about 13%.The criticality of SNEAK-7A and 7B was determined with all control rods in their most reactive position, that is, with the fueled portion of the rod in the core. The measurements were performed on critical eigenvalue (keff), and material buckling, reaction rate ratios, material worth, fission rate, and capture rate distributions as well as effective delayed neutron fraction βeff were measured. The description of the benchmark is available in the IRPhE handbook as the SNEAK-LMFR-EXP-001 evaluation. Two- and three-dimensional models for the TWODANT/THREEDANT (similar to PARTISN) and MCNP codes are included in the IRPhE evaluation. The 3D models of the benchmark models are shown on Figure 6.Figure 6
Cross cuts of detailed and simplified 3D benchmark models of the SNEAK 7A (a) and SNEAK-7B (b) fast reactor assemblies.
(a)
(b)ENDF/B-VII cross-sections (versions .0 or .1) were used in all the calculations, but in very different energy representations, that is, as multigroup with 33 to 238 groups (in deterministic codes) and pointwise cross-sections (M/C codes). Some differences were also in the computational models developed by different participants. Nevertheless, an excellent agreement between different codes was observed, both for integral values and sensitivity profiles. Thekeff sensitivities are compared in Figure 7. For most of the important reactions, the integral sensitivities agree within 1-2%. Slightly larger differences of up to around 10% were observed for the elastic and inelastic scattering, traditionally difficult since evaluated as a difference between two large (gain and loss) components causing numerical difficulties. Note also an excellent match of the sensitivities to the secondary angular distributions (SAD) calculated using SUSD3D and SERPENT codes for P1 to P3 terms of the 238U elastic scattering.Figure 7
Example of sensitivities ofkeff to 235,238U, 239Pu cross-sections, prompt fission spectra (chi), and P1 to P3 elastic scattering angular distributions (SAD) calculated by the SUSD3D, SERPENT, TSUNAMI3D, MCNP6, and SNATCH codes, including numerically demanding elastic and inelastic data sensitivities.For theβeff S/U calculations the solutions obtained using the codes SUSD3D, SERPENT, and XSUSA are compared in Figure 8. Differences are slightly larger than in the case of keff due to numerically demanding calculations of (10) and M/C statistical errors, but still very satisfactory for the validation of the βeff S/U methods. More details on the comparison of different codes and methods are available in [19, 20].Figure 8
Sensitivity ofβeff to 238U inelastic and 238U, Pu239,240 fission cross-sections, calculated using the SUSD3D, XSUSA, and SERPENT codes for the SNEAK-7A benchmark.


## 4.2. MYRRHA ADS Reactor S/U Analysis

MYRRHA (Multipurpose hYbrid Research Reactor for High-tech Applications) is being designed at SCK•CEN, Mol, Belgium, conceived to operate both in subcritical (Accelerator Driven System-ADS) and in critical modes. In the scope of the CHANDA project of the European Commission a sensitivity and uncertainty analysis of the neutron multiplication factorkeff and the effective delayed neutron fraction βeff was performed to identify the most important nuclear data for neutron-induced reactions for criticality calculations of the latest MYRRHA designs. The sensitivity profiles were calculated using different, both stochastic and deterministic, codes (i.e., SCALE/TSUNAMI, MCNP, and XSUN/SUSD3D). MCNP and TSUNAMI calculations were performed using a detailed 3D geometry model. A simplified 2D r-z geometrical model was prepared for SUSD3D (and TSUNAMI) based on the detailed MCNP model of MYRRHA. The problem was still relatively large, with the dimensions of 5 m (r-direction) and 8.5 m (z-direction), 24 material mixtures with 1257 materials at different temperatures ranging from 300 K to 800 K. 33-group cross-sections were used in PARTISN and SUSD3D calculations, pointwise data in MCNP, and 238 energy group library in SCALE/TSUNAMI codes.The XSUN-2017 sequence with the geometry is shown on Figure9. The neutron flux distributions in the MYRRHA reactor generated using the MeshTal Viewer [17] utility called from XSUN-2017 are shown in Figures 10–12.Figure 9
XSUN-2017 computational running segment with the view of the MYRRHA geometry modeling.Figure 10
Distribution of the thermal neutron flux in MYRRHA calculated using the PARTISN code.Figure 11
Distribution of the 1-2 keV neutron flux in MYRRHA calculated using the PARTISN code.Figure 12
Distribution of ~3 MeV neutron flux in MYRRHA calculated using the PARTISN code.In spite of the differences in S/U methodologies, different cross-section energy representations, and cross-section processing and treatment, the neutron-induced nuclear data sensitivity analysis resulted in differences less than 4% between codes, with few exceptions (239Pu fission, 238U elastic scattering, and 56Fe capture reactions) (Figures 13 and 14). Nuclear data covariance matrices of different libraries (SCALE-6, ENDF/B-VII.1, and JENDL-4.0u) were used to derive the uncertainty in keff based on the calculated sensitivities. This study concluded that, depending on the covariance data used, the keff and βeff uncertainties due to the uncertainty in nuclear data were about 0.5–1% and ~2.2%, respectively.Figure 13
Sensitivity ofkeff calculated using the SUSD3D and TSUNAMI codes for the MYRRHA reactor.Figure 14
Sensitivity ofβeff to the total neutron yield ofU235, 238U, 239Pu, 240Pu and Pu241 calculated using the SUSD3D code for the MYRRHA reactor.


## 5. Conclusions

A modern computer code interface XSUN-2017 was developed for the preparation and execution of the deterministic neutron-gamma computer codes in a user-friendly way suitable for today standards and users. At present, the system integrates the codes for the nuclear cross-section preparation (TRANSX-2.15), the code PARTISN for 1-, 2-, and 3-dimensional neutron and gamma transport calculations, and the SUSD3D code for the nuclear data sensitivity and uncertainty analysis. User-friendly plotting codes are also available for the 3-dimensional visualization of neutron fluxes, spectra, and sensitivity profiles. XSUN-2013 system is available through the OECD/NEA Data Bank and RSICC since early 2014. An updated and improved version, XSUN-2017, is under final testing and will be released during 2017. The code system was thoroughly validated against other sensitivity and uncertainty codes.Future development is planned including the implementation of the implicit sensitivity method (e.g., the one developed in the scope of the PhD thesis [43]) and coupling of the SUSD3D methodology to other transport codes such as ATTILE and DENOVO. Organization of workshops and training courses, such as [44], is also planned.



---

*Source: 1264736-2017-10-17.xml*