APPLICATION OF TWO STAGE METHOD OF CHARACTERISTICS / SP3 METHODOLOGY TO TRISO-FUELLED LEU SPACE REACTOR IN WIMS 11

In order to minimise the mass of a 1MWe LEU space fission power system design, a rapid neutronics analysis tool is sought. A two-stage deterministic analysis routine has been constructed using a core-plane method of characteristics calculation followed by a full-core SP3 calculation, within the ANSWERS code WIMS11. This is compared to a faster route that skips the core-plane calculation and also the Monte Carlo code Serpent. Results suggest sufficiently good agreement for the WIMS-based methods to be useful in a full system massminimising optimisation routine.


INTRODUCTION
There has been a renewed interest in the last decade in using nuclear fission power in space exploration. It has the potential to serve as a means of propulsion or as a source of electrical power for scientific apparatus, communications and in some cases, life support. The system considered here is intended primarily as a means of enabling nuclear electric in-space propulsion. Electric propulsion is now a relatively mature technology, being used as a means of orbital station-keeping and adjustment in satellites. It has also been demonstrated as a viable means of propulsion for space exploration, most notably in the Dawn mission [1] which explored the asteroid Vesta and dwarf planet Ceres. In this mission, solar photovoltaics served as the power source for the electric propulsion system which drew a power of 2.6kW. However, for larger payload missions, especially manned missions, the required power must necessarily increase. In addition, for missions venturing further than the asteroid belt, the reduction in solar flux imposes a challenge for highpower missions; for example, the flux in orbit around Jupiter is around 4% of that in low earth orbit. A potential alternative to solar power is a radioisotope power system, however the relatively low power density and high cost of Pu-238 reduces the attractiveness of this option. Certainly, for the relatively high levels of power considered here (1MWe), nuclear fission emerges as the preferred choice of power source.
Concurrent with the recent interest in space fission applications, there has been an increased interest in 'micro-reactors' for terrestrial use, with intended applications such as powering remote communities, mining installations and military bases. Therefore, methods of modelling very small reactor cores are perhaps more relevant now than at any time in the past.
The present study aims to quantify the impact of choosing to use a LEU (<20% enriched) fuel rather than the common choice of HEU. Such a choice has impacts on many important aspects of space power system design [2], but the main metric of interest here is the power system mass, which should be as low as possible. In order to determine the mass penalty of adopting LEU it is necessary to find the lowest mass power system for both the LEU and HEU cases so that a fair comparison between them can be made. This necessarily involves the analysis of a large number of potential reactor designs in a global mass-minimising optimisation scheme. Therefore, a low computational cost neutronics analysis tool is highly desirable.

REACTOR CORE DESIGN
The reactor design under consideration uses TRISO particles in a graphite matrix, with additional moderation provided by a metal hydride or beryllium. The core is prismatic, composed of hexagonal fuel blocks with seven interior coolant channels. The moderator is placed interstitially between the fuel blocks but separated from them by a hexagonal-annular coolant channel. This arrangement can be seen in Figure  1. The coolant is a 40g/mole mixture of helium and xenon, a common choice in direct-cycle Brayton space power systems.
The core is hexagonal overall with a beryllium reflector. For this type of reactor design, the design space is defined by the following variables: x Fuel block diameter x TRISO particle packing fraction x Coolant channel diameter x Width of interstitial moderator x Width of coolant channel between fuel block and moderator x Number of rings of fuel blocks x Core aspect ratio (fuel block height vs core radius) x Axial reflector thickness x Radial reflector thickness x Coolant mass flow rate x Coolant inlet temperature x Core thermal power In searching for a minimum mass power system, the core imposes some constraints that limit the design space. These include not breaching temperature limits in the core materials, the most restrictive of which is the temperature of the moderator in the case of metal hydrides. However, the constraint which is relevant here is the requirement for core criticality at the end of the mission. For the purposes of this investigation, this lifetime has been taken from the DEMOCRITOS [3] project specification, which calls for seven years equivalent at full power, where the power is 1MWe. As well as assessing criticality, the neutronics tool will need to be capable of producing a core power distribution in a format which enables coupling with a separate thermal hydraulics analysis tool.

COMPUTATIONAL ROUTES
The requirement for low computational cost makes a deterministic, rather than Monte Carlo, code the most appropriate choice. The ANSWERS © code WIMS11 [4] contains a range of mostly deterministic neutronics analysis tools, with a highly flexible environment that enables a number of possible calculation routes. Two such routes are presented here, tested in a development version of WIMS11 provided by Wood.
The two calculation routes are shown in Figure 2 in simplified form. What now follows is a more detailed description of the calculation steps.

Inner Core
The first step is the self-shielding calculation for the TRISO particles in the fuel block. WIMS11 has a module, PROCOL, for fuels randomly dispersed in spherical or annular geometries. Therefore, the fuel block, coolant and moderator are represented in simplified annular geometry. From this, collision probabilities are calculated. These are then used in a subgroup resonance treatment using the modules PRES and RES.
The second stage is to homogenise the fuel block and coolant into a single effective material. The collision probability calculation is repeated with the resonance-corrected cross sections, which are used to create a flux solution using the module PIP. The flux solution can then be used to smear the TRISO constituents, graphite matrix and coolant together. At this point, there are only effectively two materials in the inner core, a 'fuel' and the moderator.
Next the energy group structure is reduced from the XMAS 172 structure to fewer groups, in this case six. Whilst it would be possible to do this using the flux calculation used from the material homogenisation stage (based on collision probabilities in annular geometry), it was found that a relatively low-cost alternative gave somewhat better results. This involves constructing a very simple triangular geometry with rotational symmetry, holding the two materials 'fuel' and moderator. A flux solution is produced using the method of characteristics (MoC) with the module CACTUS. This flux solution is then used as the basis for energy group condensation.
Ultimately, a single homogenised inner core (everything inside the reflector) material is sought for the final full-core calculation. However, homogenising the fuel and moderator at the single fuel element level will not reflect the possibly very different environments at the core centre versus the periphery, next to a somewhat moderating reflector. Therefore, a 2-D core plane geometry is created in CACTUS. A few-group flux solution is produced which can then be used as the basis for the final material homogenization of the fuel and moderator. This is done for each core ring separately so that there are as many versions of the inner core material as there are rings in the core.
The final stage is the full-core eigenvalue calculation. Reflector material data in the few group energy structure is prepared separately (see below) and placed in an RZ geometry with the inner core material. This full-core calculation is done using the SP3 method within the MERLIN module.
Given the large core design space, there is value in having a faster route available to carry out low fidelity calculations that narrow the design space. In the main route described above, the core-plane CACTUS calculation is the most computationally expensive element. Therefore, an alternative route has been constructed that bypasses this step. This is also shown in Figure 2. In this route, all inner core materials are homogenised spatially at the earliest possible stage and converted to the few group energy structure. This single material is then used to fill the entire inner core region in the SP3 full core calculation.

Reflector Data Preparation
The reflector is only a single material so does not require any spatial homogenisation. However, a suitable flux solution is required to carry out energy group condensation. It was found that a relatively simple geometry gives final results that differ negligibly from more complex geometries. The approximate geometry used is a 1-D MoC calculation in CACTUS for the core which captures the relative sizes of the inner core and the axial reflector. Data for the reflector is produced for five concentric reflector regions to capture some of the differences in neutron energy spectra that would be expected to appear across the reflector's width.

Burnup
When carrying out burnup calculations, the resonance self-shielding calculations must of course be carried out at each burnup step. One possible approach to this is to combine all burned materials into a single 'average' material, which is used to carry out a self-shielding calculation. The microscopic cross sections arising from this are then combined with burned isotopic compositions to form the basis of the next flux calculation. This approach is typically used in WIMS for assembly-level calculations, which is its traditional application. However, in a full core calculation (especially for a very small core) the assumption that an average effect of depletion on self-shielding is suitable is less tenable. Therefore, a separate selfshielding calculation should take place for each burnable region, along with the subsequent steps which prepare few-group, spatially homogenised information for the full core calculation. This scheme is shown in Figure 3. Following the full-core calculation, the calculated fluxes in each burnable region must then then be translated back onto the most detailed level of geometry. This is the reverse process of the prior energy and spatial homogenisation and is carried out using the same flux solutions that were used to carry out the original homogenisation. This requires storing a separate version, one for each burnable region, of the finest level geometry, onto which the appropriate flux from the full-core calculation is applied. Once this is done, the burnup module ('BURNUP') is applied to each of these versions. The burned materials are then fed back to the self-shielding calculation for the next step.
It should be noted that in the WIMS models, a separate burnable region for each core ring is assigned but there is no axial division. For the 'quick' route it is relatively trivial to include additional axial burnable regions. For the main WIMS route, however, adaptation of the 2-D core plane MoC calculation to different axial burn regions requires further investigation. One possibility is having a separate core plane calculation for each burnable layer. This will increase computation time noticeably as this stage generally consumes the majority of the overall calculation time. The layers nearest the ends of the core might also need to be adapted to include the effect of the axial reflectors. This additional development will be pursued if the main route proves useful in the overall optimisation scheme. The continuous-energy Monte Carlo code Serpent [5] is used to produce reference results against which the WIMS-based methods can be compared. Note that these comparator Serpent runs have assembly-wise burnable regions but only one axial burnable region for each of these, to provide a fair comparison.

Computational cost comparisons
The purpose of pursuing a deterministic neutronics analysis method was to achieve a rapid calculation time compatible with a global optimisation method seeking to find a minimum mass space power system. The main comparison of interest here is between the two WIMS routes but Serpent run times are also given, with the caveat that different machines are used for the two codes and Serpent was configured with an emphasis on accuracy rather than speed. WIMS runs are carried out using a single thread on an Intel Xeon E5 processor at 2.6GHz whilst Serpent runs are carried out using 32 threads on Cambridge University's CSD3 High Performance Computing system. Sample run times for the case of a 0.4cm interstitial moderator are shown in Table 1.

keff comparisons
The results in terms of keff at a variety of moderation levels, using a zirconium hydride moderator, at the end of life are shown in Figure 4. These results are all for a core with eight rings of hexagonal fuel/moderator assemblies, a fuel block radius of 2.0cm, 8cm radial and axial reflector, and a core radius/height aspect ratio of 0.5. The fuel kernels are 20% enriched UO2. This configuration results in a maximum burnup of approximately 100GW/t. The beginning of life keff estimates are broadly similar between the two WIMS routes up to a moderator thickness of around 0.4cm. With larger amounts of moderation, the main WIMS route shows its utility, maintaining an error vs Serpent of below 1,000pcm whilst the quick route diverges to a maximum error of 3,000pcm.

CONCLUSIONS
Two possible computational routes within the WIMS11 suite have been explored for use in modelling a small LEU-fuelled thermal spectrum reactor intended for space applications. These show sufficiently good agreement with Serpent to be useful in a full power system mass-minimising optimisation routine. For heavily moderated designs, the more computationally expensive route, which features an intermediate coreplane MoC calculation performs better than the quicker route which omits this step. One possible implementation of these schemes is to use the quick route in an initial approximate optimisation process, after which any minima found can be investigated further with the higher fidelity main route.

ACKNOWLEDGMENTS
The authors are grateful for the funding provided for this research by the National Nuclear Laboratory and the UK Engineering and Physical Sciences Research Council through ICO Centre for Doctoral Training, along with the Cambridge Centre for Data-Driven Discovery for access to the CSD3 high performance computing system.