Benchmarking of the Serpent 2 Monte Carlo Code for Fusion Neutronics Applications

Analyses of radiation fields resulting from a deuterium-tritium (DT) plasma in fusion devices is a critical input to the design and validation of many aspects of the reactor design, including, shielding, material lifetime and remote maintenance requirements/scheduling. Neutronics studies, which perform in-depth analysis are typically performed using radiation transport codes such as MCNP, TRIPOLI, Serpent, FLUKA and OpenMC. The Serpent 2 Monte-Carlo code, developed by VTT in Finland, is the focus of this work which seeks to benchmark the code for fusion applications. The application of Serpent 2 in fusion specific analysis requires validation of the codes performance in an energy range, and a geometrical description, which significantly differs to conventional nuclear fission analysis, for which the code was originally developed. A Serpent model of the Frascati Neutron Generator (FNG) Helium Cooled Pebble Bed (HCPB) mock up experiment has been prepared and the calculated results compared against experimental data, as well as the reference Monte Carlo code MCNP. The analysis is extended to a model of DEMO with HCPB blanket concept. For this model, the flux, nuclear heating, tritium production and DPA are calculated, all of which are integral nuclear responses in fusion reactor analysis. In general, a very good agreement is demonstrated for both of the benchmarks, with any discrepancies pinpointed to different physics models implemented.


INTRODUCTION
Monte-Carlo N-particle (MCNP) [1] is a comprehensive code for reactor safety and shielding analysis, and is, at present, the most commonly adopted code by the fusion neutronics community. This code has been extensively validated and actively maintained over its˜40 year history. There is a specific set of requirements that a radiation transport code must demonstrate to be applicable to fusion relevant analysis: neutron and photon coupled radiation transport using point-wise cross section libraries, able to provide a geometric representation of the modelled system in all its complexity, capable of employing acceleration techniques, accommodate complex plasma neutron source definitions, and have parallelisation capability for deployment on high performance computer architectures stemming from the computer intensive nature of Monte-Carlo simulations [2]. MCNP complies with this list of fundamental requirements hence its prevalence in fusion neutronics and as the benchmark code for the work presented herein.
The development of Serpent began in 2004 as part of a doctoral thesis at the VTT technical research centre, Finland. The code now supports many of the fundamentals highlighted above, prompting detailed investigations and extensive benchmarking of the code in the nuclear fusion domain. The current version of the code, Serpent 2 [3], has been applied to ITER and DEMO analysis in previous studies [4] [5]. To generate Serpent models, the current methodology is to use a script which automates the conversion from an MCNP input [6]. The ensuing process of Monte Carlo code validation is then performed using published experimental data. SINBAD * [7] provides a large database of experiments which cover accelerator shielding, fission reactor shielding and nuclear fusion. In this paper, focus is given to one such experiment [8], the Frascati Neutron Generator (FNG), Helium Cooled Pebble Bed (HCPB) mock up. FNG is a 14 MeV neutron source used to simulate the deuterium tritium reaction, the favoured nuclear fusion reaction. The HCPB blanket concept is one of the selected designs under Test Blanket Module (TBM) framework in ITER, which will test the feasibility of this particular concept against others. The SINBAD distribution includes an MCNP model of the experimental set up with tallies for direct comparison to the experimental data.
Subject to the demonstration of its viability, the HCPB blanket concept will potentially be adopted in EU DEMO, the European successor to ITER. As such, radiation transport models of DEMO with this blanket geometry have been created for simulation purposes. In the mapping of the radiation environment during both operation and in the shutdown phase, nuclear responses must be demonstrated to adhere to specified limits dictated at a higher level by the regulating authority. The neutronics assessment is therefore an integral part of the complete integration design process. At the time of writing, validation of Serpent 2 for this application is confined to within the vacuum vessel beyond which the statistical uncertainty is too high with analog transport.

Methodology
All calculation results were performed using MCNP6v1 [1] and Serpent v2.1.31 beta [3]. Serpent was originally developed as reactor physics code for performing burnup calculations, spatial homogenisation and modelling of research reactors. Its application now extends to coupled neutronphoton multiphysics, and more recently radiation transport, particularly in the sphere of nuclear fusion. As the scope of the code evolved, new features were implemented aligned with the expanding list of requirements. One such example is support for CAD based geometry types. This is potentially a more efficient method for CAD-centric work flows with complex models such as those in fusion reactor analysis. The conversion of such models to constructive solid geometry (CSG) can be error prone and time consuming. It is nonetheless the representation used here on the basis that the MCNP CSG model is the reference starting point.
The reference nuclear data library used for neutron transport for the FNG HCPB benchmark is FENDL-2.1 [9]. For all DEMO analysis, the neutron transport library is JEFF-3.2 [10]. Dosimetry cross section libraries were used for the activation foils, namely IRDFFv1.05 [11]. The adopted photon library in all cases is MCPLIB04/84 [12].
The parametric plasma source description for DEMO was re-written as a C routine. This is called in Serpent as a user defined source. Similarly, the MCNP FNG source routine is written in Fortran and therefore can not be directly invoked by Serpent. Instead, a list of starting source particles with position, direction, energy and weight was generated in MCNP and subsequently called in Serpent. All simulations were performed to 10 9 histories on an Intel Xeon E5-2665 cluster, using coupled neutron-photon particle transport.

FNG HCPB mock up
In the mock-up experiment, the tritium production in breeding cassettes and reaction rate in a set of activation foils was determined. A comparison of the computed results in MCNP and Serpent against the experimental data, aims to validate both the implemented particle transport model as well as the nuclear data. An illustration of the Serpent model detailing the experimental configuration of the mock up can be seen in Figure 1. The Serpent model of the HCPB mock up was produced through automatic conversion of the MCNP input file. Beryllium is included in the model as the neutron multiplying medium, shown in light blue/cyan in Figure 1. Nb, Au, Al and Ni activation foils were located at increasing distance from the source through the mock-up. The reaction rates for the following activation reactions were determined experimentally: 93 Nb(n,2n); 197 Au(n,γ); 27 Al(n,α); and 58 Ni(n,p). The tritium production rate (TPR) was also determined in breeder cassettes, again positioned at different depths in the mock-up, each containing 12 individual pellets of Li 2 CO 3 powder. At each of the foil locations, the reaction rate was calculated in MCNP and Serpent. The most notable differences between the calculated and experimental results, as plotted in Figure 2, occur in the centre of the block. The experimental errors include counting statistics, detector uncertainty and FNG source calibration. In general there is very good agreement between the two calculated responses, with the largest deviation equal tõ 1.5%. To try and understand the differences to the experimental results and assess the sensitivity to different dosimetry cross sections, the calculations were repeated using the previous IRDFF revision, IRDFF 2002, and the LLDOS library [13]. For all isotopes, the two IRDFF libraries agreed within the bounds of uncertainty. For LLDOS, only Aluminium demonstrated agreement with the IRDFF libraries, with on average a factor of 2 discrepancy for all other elements. The total activity of tritium in each of the Li 2 CO 3 pellets is the sum of the activities of the respective nuclides, 6 Li and 7 Li. The production channel for tritium following neutron irradiation of these two nuclides is via two key pathways: The calculation results are appropriately normalised to factor the total neutron production rate, mass of each pellet and the decay of tritium. For the purposes of this comparison, the lower stack of pellets was chosen, nominally ENEA 2 through ENEA 8, corresponding to increasing distance from the source. Figure 3 compares the computed results in MCNP and Serpent against the experimental data. Excellent agreement in the calculated responses is seen for all of the pellet stacks. As for the activation foils, larger deviations are seen towards the geometric centre of the mock up, which requires more detailed investigation. The codes are consistent in under predicting the experimental data, which is a possible artefact of the transport library used -more recent revisions of the FENDL library are available [14].

DEMO HCPB
To reduce the computational load, fusion reactor analysis takes advantage of the toroidal symmetry of tokamaks by modelling a single repeating sector. Serpent requires that the sector is housed in a universe which is unfolded into a 360 • model. This is analogous to the reflecting boundary implementation in MCNP (see Figure 4). The DEMO HCPB MCNP neutronics model is a 10 • sector and includes a heterogeneous representation of the blanket modules, with the breeder zone, cooling plates, and neutron multiplier modelled explicitly. A universe lattice implementation in MCNP was replicated in Serpent to model this repeating structure. The conventional method to quantitatively validate the conversion of the model, and ensure that transport critical masses are preserved, is to perform a volume comparison. One can then be confident that all deviations in nuclear responses are inherent to the physics models (and nuclear data called) in the code itself. Both codes include a stochastic method for estimating the geometric volumes. The largest deviation was 3.5% in the total mass of SS3316. On looking more closely at a cell-wise comparison of the volumes, it was determined that the deviation is an artefact of the statistical error on the volume of smaller features and not an error in conversion. Serpent calculations were performed using a hybrid MPI-OMP parallelisation to optimise memory consumption and simulation time. A map of the relative error is given in Figure 5, where it is visible that results are converged for the blankets and divertor region only. Low statistical uncertainty for nuclear responses in the ex-vessel region necessitates the use of variance reduction methodsan active area of development of Serpent of which a robust implementation is integral to fusion applications.
The neutron and photon spectrum in 175 energy groups was determined for a cell in the plasma facing layer of the divertor, plotted in Figures 6 and 7 respectively. For the the ratio of the two calculated values, a metric is provided for highlighting data points where the deviation exceeds two times the standard deviation on the result.  For both neutrons and photons, the general shape of the raw spectrum is in good agreement with no conspicuous outliers. There is however a notable number of photon spectrum data points where the deviation is not within the 2σ bounds, particularly for energies below 1 MeV. This is a result of the different photon physics models used in the two codes. An example is for the thick target bremsstrahlung approximation, where Serpent correctly accounts for a higher photon yield.
Further, both codes have different predictions for photon annihilation.
Further in-vessel nuclear responses ( Figure 8) have been tallied and compared to the MCNP calculated results. The nuclear heating is significant for many tokamak components, the DPA for material lifetime, and the tritium production assesses the rate at which tritium fuel can be recycled in the reactor. These responses are typically calculated for the blanket modules where globally, peak values are observed. Large deviations are seen in the neutron nuclear heating in the tungsten layers of the divertor. After further investigation it was noticed that there is a significant number of negative cross section values in the ACE files in the energy range, 1-15 MeV. This stems from non-physical data in the original data evaluation prior to any processing. The issue of negative entries in ACE files is reported [15] to extend to DPA cross sections, however only iron was considered here. The two codes employ different routines for handling negative values which translates to the different reported values, however, significantly, only Serpent notifies the user of the erroneous data.

Conclusions and Further Work
The Serpent 2 Monte-Carlo code has been benchmarked against MCNP and experimental data in the SINBAD database. Serpent inputs for the FNG HCPB mock up experiment and the EU DEMO HCPB model were generated using a script which converts the MCNP model. The comparison of nuclear responses in DEMO serves an important demonstration of the applicability and validity of Serpent to fusion reactor analysis. Typical responses were selected in the breeder blanket including neutron spectrum, nuclear heating, displacements per atom (DPA) and the tritium production. In general, very good agreement is observed between MCNP and Serpent for all nuclear responses considered. The one anomalous value was observed in the nuclear heating in tungsten. The origin of this is non-physical data at the nuclear data source, the processed result of which is handled differently by MCNP and Serpent. The erroneous data was only flagged by Serpent, which serves as a minor motive for exploring alternative transport codes.
The statistical uncertainty increased significantly through the vacuum vessel necessitating variance reduction methods to tally ex-vessel responses. Weight window based variance reduction techniques are currently an active area of development in Serpent and a robust implementation is a fundamental inclusion for most radiation transport applications. Current testing is focused on the built in response matrix solver method [16] which can potentially provide a highly efficient means of generating weight windows.