Monte Carlo Simulation of the UK’s First EPR Nuclear Reactor Startup Core using Serpent

Computationally modelling a nuclear reactor startup core for a benchmark against the existing models is highly desirable for an independent assessment informing nuclear engineers and energy policymakers. This work presents a startup core model of the UK’s first Evolutionary Pressurised Water Reactor (EPR) based on Monte Carlo simulations of particle collisions using Serpent 2, a continuous-energy Monte Carlo reactor physics burnup code. Coupling between neutronics and thermal-hydraulic conditions with the fuel depletion is incorporated into the multi-dimensional branches, obtaining the thermal flux and fission rate (power) distributions radially and axially from the three dimensional (3D) single assembly level to a 3D full core. Shannon entropy is employed to characterise the convergence of the fission source distribution, with 3 billion neutron histories tracked by parallel computing. Source biasing is applied for the variance reduction. Benchmarking the proposed Monte Carlo 3D full-core model against the traditional deterministic transport computation suite used by the UK Office for Nuclear Regulation (ONR), a reasonably good agreement within statistics is demonstrated for the safety-related reactivity coefficients, which creates trust in the EPR safety report.

optimisations of the EPR compared with the recent advances in other PWR designs (e.g. AP1000, Sizewell B) to inform the regulatory authorities.
The primary goal of modelling a nuclear reactor core is to reliably predict the reactivity and power distribution accurately and efficiently. The mechanism of the full-core computational modelling is depicted in Figure 1. First, few-group homogenised cross sections for 2D individual assemblies are generated based on solving neutron transport equations. Subsequently, multi-parameter cross sections dependencies are utilised for nodal diffusion calculations and the coupling between neutronics and thermal-hydraulic feedbacks with the evolution of fuel depletion in the 3D full-core model. Rates of the neutron production and loss are closely involved in the prediction of the reactivity and power distributions. This reflects the significance of the Boltzmann transport equation [8], an integro-differential equation in 7 dimensions, representing the neutron balance, i.e. the number of neutrons generated, equals to those consumed by fission, capture and other losses. It can be solved by both the deterministic approach [9] and the stochastic method [10]. The deterministic approach copes with it in a numerically approximated manner, such that the continuous energy dependence of cross sections is condensed into discrete energy groups (i.e. multi-group approximation), with the angular dependency handled by discretisation or functional expansions [9]. How precise nowadays the neutronics and via this the power distribution of the large core of the 1600 MWe EPR can be simulated is of academic research interest. It is important that the neutronics of the EPR is simulated by the state-of-the-art Monte Carlo codes. As a 3D continuous-energy Monte Carlo burnup code, Serpent is specialised in lattice physics calculations for generating homogenised group constants [10]. As compared with traditional deterministic methods, the state-of-the-art Monte Carlo method relies on fewer approximations, e.g. no discretisation of space or energy required. However, this results in a higher computational cost and a longer running time [11], especially for 3D full-core burnup calculations. Nevertheless, it is of extreme importance to provide the nuclear regulatory body with an independent consideration of the neutronics and thereby the nuclear safety of the EPR, as an independent evaluation and benchmark can create trust in the safety reports. Thereby, using Serpent 2 as a reference in this work to benchmark against conventional deterministic codes could gain confidence if the cycle length claimed by the designers can be achieved while satisfying the safety-related and operational constraints on power peaking, reactivity coefficients, and maximum boron concentrations. Section 2 of the paper presents the key elements of the Serpent-based Monte Carlo model setup, including the core geometry descriptions, the multi-dimensional branches covering the burnup and thermal-hydraulic operating conditions, the modelling of self-shielded burnable poison pins and heavy reflectors, the convergence of Shannon entropy, and the source weighting. Section 3 reports the Serpent simulation results from the 2D assembly, 3D assembly to the 3D full core. Relative fission rate and thermal flux distributions are visualised. Safety-related reactivity coefficients and power peaking factors are computed and validated with those obtained by deterministic codes submitted to the ONR. Negligible deviation does exist due to the statistical nature of the Monte Carlo method, which is further mitigated by simulating more neutron histories and implementing the source biasing.

Serpent Monte Carlo Model Setup
To obtain the flux and power distributions of the startup core, the modelling follows the steps from the lattice transport to the diffusion core calculations. The principal route is the homogenisation by preserving the reaction rates, tracking the isotopic changes in materials with burnup, parameterisation in branches, and multi-dimensional table interpolations. Proper model parameters are strategically chosen, allowing a high modelling accuracy without degrading the computational efficiency.

Full-core Geometry Specification of EPR
The fresh core model developed in this work adheres to the design specifications as listed in the Pre-Construction Safety Report [4,[12][13][14][15][16]. The whole-core configuration with 241 assemblies is precisely represented in Figure 2 (radial view) and Figure 3 (axial view). Each assembly is a 17×17 square lattice composed of 264 fuel rods with highly corrosion-resistant M5 alloy claddings, as well as 25 guide thimbles for the Rod Cluster Control Assembles (RCCA) and in-core instrumentations. The top and bottom of the fuel assemblies are smeared regions of the construction materials.  To partially represent the thermal-hydraulic conditions using Serpent, relatively hotter water (density of 0.6456 g/cm³ at the outlet temperature of 605 K at hot full power) is deployed at the upper part of the core as shown in Figure 3, while relatively colder water (density of 0.7524 g/cm³ at the inlet temperature of 563 K at cold zero power) is deployed at the lower part of the core. Note that the coolant temperature in the core's mid-plane at the beginning of cycle (BOC) is lower than the average coolant temperature, as the axial power peak is biased towards the bottom half of the core due to the more neutron thermalisation (the analysis is detailed later in section 3).
Note that the traditional baffle assembly is replaced with a heavy radial reflector [4] made of thick steel slabs (95.6% by volume) and water (4.4% by volume), which are modelled as full-size assemblies and homogenised with the barrel and downcomer water in this work. Surrounding the fresh core with a heavy reflector can boost the neutron economy by minimising the leakage of neutrons out of the core [4], increase the burnup with a given enrichment (or reduce the enrichment requirement) [12], secure a flatter flux profile, and protect the reactor pressure vessel against the fast neutron fluence-induced aging and embrittlement [4].
The fresh core's fuel loading in this work is based on a 5-regions out-in checkerboard loading pattern as specified in Figure 4. To reduce the power peaking at the beginning of life (BOL), the highest enriched fuel assemblies (4.2wt%) are loaded in the peripheral region, while multiple relatively low enrichment assemblies (2.1wt%) are interspersed in the core interior in a checkerboard fashion.

Multi-dimensional Branches of Diverse Operating Conditions
The calculation is firstly tasked with deriving the multiplication factor as a function of burnup for 18 months cycle. Each assembly is irradiated to 80 GWd/t (discharge burnup). Since the Xenon transient happens relatively fast, a small burnup step at the beginning is selected, i.e. 0.1 GWd/t to get the Xenon to equilibrium. For the criticality calculations and obtaining the safety-related reactivity coefficients, multi-dimensional branches of different operating states are implemented.
First, nodes at different positions of the core are distinct in the operating conditions. First, the moderator density ρ(mod) varies with the core height, i.e. a higher moderator temperature T(mod) near the top as opposed to that at the bottom. Thus, the multiplication factor has a functional dependence on T(mod) and ρ(mod). With the reference coolant temperature at the core average value of 580 K, branch points from 0.7524 g/cm³ at the inlet temperature of 563 K (cold zero power) to 0.6456 g/cm³ at the outlet temperature of 605 K (hot full power) are employed in this work.
Second, the criticality changes with the fission products building up and the fissile material depletion. To compensate for this change, boron is diluted progressively with fresh water, thus the multiplication factor must reflect this reduction in the boron concentration (BC) with burnup. Branch points of the boron concentration from 0 to 2000 ppm are employed.
Third, the fuel temperature T(fuel) varies with the core height (hotter in the middle as per the power profile). Thus, the multiplication factor also exhibits a dependence on T(fuel). With the reference operating fuel temperature of 900 K (hot full power), the branch points from 340 K (cold zero power) to 1500 K are chosen in the multi-parameter tabulations. The insertion and withdrawal of the Rod Cluster Control Assembles (RCCA) influence the multiplication factor as well. To obtain the safety-related reactivity coefficients over cycle 1, this work starts from parameterising a single dependence on burnup, and subsequently branching out in diverse directions of the operating conditions, i.e. T(mod), ρ(mod), T(fuel), BC, and RCCA in/out, obtaining multi-dimensional dependencies upon the operation states as per the regulations [12][13][14][15][16].

Modelling Gadolinia Burnable Poison Pins and Heavy Reflectors
Precisely modelling the burnup behavior of gadolinia (Gd) burnable poison pins in a nuclear reactor is tricky, as it is a very strong absorber of thermal neutrons. The highly self-shielded burnable poison depletes largely from the outermost zones inwards, presenting strong flux gradients around the pin. Note that self-shielding in energy reduces the neutron capture in the resonance of fertile material, while spatial self-shielding is caused by the fact that some of the thermal neutrons entering the fuel from the moderator are absorbed near the surface of the fuel. They do not survive to contribute to the thermal flux in the interior. The outer portion of the fuel thus shields the interior, indicating that modelling the whole Gd-bearing pin in only one zone would potentially overestimate the poison's burning rate at the BOL.
Note that among the isotopes of Gd, Gd-157 and Gd-155 exhibit the highest neutron absorption cross sections, hence contributing to the most of the thermal neutron absorptions. To characterise this self-shielding effect in the EPR Gd-bearing pins, the fuel pin admixed with Gd is divided into 10 annular segments of an equal area as depicted in Figure 5. With this method, the density variation of Gd-157 and Gd-155 versus the assembly burnup is investigated using Serpent. Figure 6 sketches the results of the atomic densities of Gd versus burnup, where Gd-155 and 157 burn out from about 10 GWd/t. By capturing a neutron, Gd-155 transmutes to Gd-156. Albeit with small capture cross sections, Gd-156 serves as a constant source to feed Gd-157 with very high absorption cross sections, i.e. Gd never burns down completely, resulting in a residual poison penalty.  radial distance at different burnup stages. The observed burning down layer by layer from the outermost rings inwards quantifies the rate of the spatial self-shielding behavior for the EPR's Gd-bearing fuel pins. Precise modelling the heavy reflector requires additional consideration, as errors can arise from the high flux gradients between nodes. A single-assembly infinite lattice model cannot produce accurate prediction because it requires the incoming current of neutrons with energy distributions representative of the actual interface between the core and the reflector. Thereby, characteristic interaction parameters for peripheral fuel assemblies facing the radial reflector and the radial reflector itself are produced using a fuel-reflector super cell model including the interfacial interactions. The simulated thermal flux and fission rate distribution of different super cell combinations are presented in Figure 8 at the BOL under hot full power (HFP), illustrating the effect of neutron thermalisation on fission rates. The relative fission power distributions are distinguished in "hot" shades of red and yellow, while the relative thermal flux distributions are represented by "cold" shades of blue (energy < 0.625 eV). The magnitudes are visually linked to the luminance, namely, brighter areas suggest a better moderation and higher power, while darker spots imply a lower local flux and power (e.g. Gd-bearing rods in dark red). The intra barrel-vessel area is occupied by water at the system pressure and temperature. The reflector, which is non-fissile, is heated by the absorption of gamma radiation. Therefore, the water temperature and boron branches are implemented by using the three different water densities, i.e. 0.7524 g/cm³ at the inlet temperature of 563 K, 0.6456 g/cm³ at the outlet temperature of 605 K, and 0.09841 g/cm³ at 620 K (steam only at the operating pressure).

Convergence and Fission Source Entropy
For Monte Carlo criticality problems, the effective multiplication factor might converge sooner than the fission source distribution [17]. The convergence behavior of Shannon entropy [17] sheds light on the suitable number of inactive cycles that should be discarded before starting the Monte Carlo tallies. Neutron histories are summarised in Table 1 based on the convergence behavior of the multiplication factor as well as the Shannon entropy as illustrated in Figure 9 for the 3D full core at the BOL under HFP. Source convergence takes a longer time for the full-core calculations. The full-core geometry is covered with a cartesian mesh (17×17×21). In this mesh, the distribution of collision points and fissions are collected during the first 200 inactive neutron cycles. The distribution is subsequently used to adjust the number of emitted fission neutrons during the active cycles. The jump of the fission source entropy at the cycle 200 is illustrated in the cumulative values and errors. This is the point where active cycles start, and the statistics are reset.  Figure 9. Convergence of fission source entropy for the EPR 3D full-core simulation using Serpent.

Variance Reduction by Source Biasing
The modelling accuracy of Serpent can be improved by running more neutron histories and cycles but at the cost of a considerably increased computation time. There is a huge potential for the variance reduction [18] at a reasonable cost of the computational burden. In our EPR model, the upper Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 5 August 2020 doi:10.20944/preprints202008.0111.v1 corners of the core exhibit larger statistical errors because of fewer reactions, with a longer computation time accordingly for these regions to converge during the power iterations. The statistical accuracy could be cleverly improved by a uniform fission source method [18], i.e. source biasing more neutron samplings in the less important regions with low fission power, and adjusting the collision weights for a statistically equivalent collision density distribution. As sketched in Figure 10 using the Serpent's mesh plotter for the EPR 3D full-core normalised collision density, higher statistical weights are assigned to the centre of the core as shown in Figure 10 (b) (with source biasing) to reduce the regional uncertainties and statistical noises occurring in Figure 10 (a) (without source biasing).
(a) (b) Figure 10. Serpent mesh plots (radial view) of the EPR full-core relative collision density: (a) no source biasing; (b) with source biasing.
When simulating collisions with an analogue Monte Carlo method, all particles share an equal statistical weight, i.e. the weight turns from 1 to 0 when absorbed. Whereas for the implicit capture that Serpent is currently using, the particles gradually lose their statistical weight as they transport through the system. There is a cut off below which the particle will be killed. Multiplying the number of collisions by the statistical weight, the collision density distribution is obtained accordingly as shown in Figure 10 above. The implicit capture extends the length of history, allowing particles to have more collisions at these corners, increasing the number of events, and thereby reducing the variance.

2D and 3D Assemblies Results
WIMS (a deterministic lattice transport code) is employed to validate the reliability of the 2D assembly results obtained by Serpent. It is expected from the modern lattice codes to predict a single assembly's multiplication factor within a few hundred percent milli-rho (pcm) of accuracy. Identical fuel compositions (in terms of atomic number fractions) and operating conditions are set for the two codes to avoid inconsistencies. Both the assemblies with and without Gd are examined with respect to burnup. Finer irradiation sub-steps are chosen at the BOL to precisely capture the relatively fast changes in the material compositions. The results for the 2.1wt% assembly (with 8 Gd of burnable poisons) are shown in Figure 11, characterised by using a higher number of discrete time steps at the BOL to accurately characterise the depletion of Gd. The comparison shows a good agreement. Table  2 below quantifies the deviation of the multiplication factors between WIMS and Serpent regarding all the 5 types of EPR fuel assemblies burning to 0.5 GWd/t at the BOL (Xenon building up to the Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 5 August 2020 doi:10.20944/preprints202008.0111.v1 equilibrium level, RCCA out, and 0 ppm boron). The deviation ranges from -0.7 pcm to 62 pcm, which is largely within statistics between the two codes. Figure 11. Serpent vs. WIMS results benchmark for the EPR 2D assembly (2.1wt%+8Gd), depletion analysis up to 80 GWd/t. Applying an axially and radially reflective boundary condition, the validity of the Serpent 3D assembly model is checked against PANTHER (a 3D nodal diffusion solver whose multi-dimensional cross sections are produced by WIMS). Following the same calculation settings, the difference between PANTHER and Serpent in predicting the multiplication factor versus burnup is largely within 100 pcm over the first cycle (as evidenced in Table 3 below), hence gaining confidence for the Serpent 3D single assembly model setup. Thermal-hydraulic conditions are incorporated into the EPR 3D single-channel modelling. 21 axial layers (with the height of 0.2 m for each mesh) are employed to simulate the axial power distribution over the core life as reported in Figure 12. At hot zero power (HZP), the axial distribution of the water density is nearly uniform, resulting in the flux peak at the centre of the core height (a cosine-like shape as expected). When it comes to hot full power (HZP) at the BOC, coolant temperature varies with the core height, i.e. the outlet is higher than the inlet. Due to the thermal expansion, the atomic number ratio of the moderator to fuel decreases with the core height. The decreasing coolant density from the bottom to the top results in a less efficient moderation above the core mid-plane, and hence the neutron spectrum is hardened and the leakage is increased, shifting the flux peak from the centre towards the lower part (at about 40% of the core height). With more fission reactions occurring due to the higher thermal flux, the fuel at the bottom achieves a higher burnup. The axial power re-distribution occurs with the bottom fuel burning down faster and a higher concentration of Xenon building up at the bottom compared with the fuel at the top. This progressively leads to another skewing of the flux profile, i.e. peaking towards the top as evidenced in the results for the middle of cycle (MOC) and the end of cycle (EOC) in Figure 12.

3D Full-core Normalised Fission Rate and Thermal Flux Distributions
The whole-core depletion using Serpent is based on 10 radial rings for each pin with Gd, and 21 axial layers for 3D analysis. Tracking 3 billion neutron histories in the Serpent full-core simulation run on a high-performance computing multi-core parallel system, a reasonable source convergence is achieved in the power iteration procedure, as evidenced in Figure 13 (a) for the symmetry of the normalised thermal flux and fission rate (power) radial distributions. The convergence behavior of the 3D fission source distribution is evidenced in the Shannon entropy already shown and explained in Figure 9. For the whole-core axial flux profile, the variance is considerably mitigated for the power distribution at the core edge, with Figure 13 (b) presenting a bottom-peaked skewed power shape at the BOC under hot full power (HFP), which agrees with the 3D assembly axial power profile at the BOC as shown in Figure 12. The variation of the moderator density with the core height together with the non-uniform axial burnup effect contribute to this asymmetric (bottom-peaked) axial flux distributions. Under HFP at the BOC, the moderator temperature is higher near the top as opposed to that at the bottom. The lower water density and the resulting less efficient moderation along with the higher leakage reduce the thermal flux at the higher core positions, thus skewing the peak of the power shape towards the bottom as evidenced in Figure 13 (b). Due to the limitation of the Serpent's mesh plotter (i.e. reaction rates are normalised to the maximum values but without presenting scales), detectors are defined to obtain the numerical data. Figure 14 presents the results of channel form factors varying with time in the first core cycle. The power peaking factors exhibit fluctuations with the burnup, as it is a complicated process encompassing the depletion of Gd, buildup of fission products and plutonium at different regions of the core. Burnable poisons have burnt out in all types of assemblies since about 15 GWd/t, and hence mitigating the form factor thereafter.

Benchmark and Independent Evaluation
Reactivity coefficients constitute an important part of the safety and licensing case. The Doppler temperature coefficient (DTC) in this work is calculated by a uniform perturbation of the fuel temperature (varying from 340 K to 1500 K) independently, whilst freezing the coolant temperature (and density) to separate the effects of the fuel heat-up and the coolant heat-up. Uncertainty propagation is employed for the Serpent results. A similar perturbation mechanism is applied for obtaining the moderator temperature coefficient (MTC). The initial boron content (zero burnup, no Xenon) search to maintain the core's criticality is performed for the first core cycle at the nominal HFP state (RCCA out). To benchmark and independently review the full-core modelling results submitted to the Office for Nuclear Regulation (ONR), Table 4 presents a quantitative comparison of the simulation results obtained by the single-stage Monte Carlo Serpent code (this work) versus those calculated by multi-stage deterministic codes APOLLO-SMART-FLICA III-F [3]. The comparison results show that the total enthalpy rise hot channel factors and the critical boron search exhibit negligible difference between the Monte Carlo and the deterministic predictions, while the total heat flux hot channel factor exhibits a more significant difference (up to -4.61%). The higher deviations in the doppler temperature coefficient as well as the moderator temperature coefficient could be due to the different perturbation methods, diverse treatments of the scattering kernel S (α, β), as well as the different upper limits of the highest energy boundary based on the different nuclear data library employed (i.e. up to 25 MeV for Serpent based on ENDF/B-VII [19], while up to 10 MeV for APOLLO based on JEF-2.2 [20]). Arguably, it would be more informative to employ the same nuclear data library in both methods for future comparison work, so that the discrepancy due to the data source could be factored out. Another implication from the core modelling study is associated with the core dimensions. Although the larger the surface area, the larger leakage, it is the surface-to-volume ratio that matters. The EPR large core offers a low surface-to-volume ratio of 1.79 and a relatively low leakage, hence contributing to a better neutron economy in contrast to the conventional PWRs. However, to propagate a disturbance across a large core requires many neutron generations, i.e. the core might become neutronically decoupled. It takes some time to propagate perturbations in one part of the core, which might exhibit Xenon-induced power oscillations. Power shape is targeted to be as flat as possible for the sake of the thermal performance optimisation. Progressive enhancement work might follow concerning the optimisation of the loading patterns [21] and the zoning of poisoned rods [22]. Low leakage [23] in/out fuel management regime could be examined. Furthermore, the optimisation between the modelling accuracy and the computational burden can be attempted by strategically adjusting the branch points.

Conclusions
For the first time to the best of our knowledge, this paper develops and verifies a startup core model of the UK's first EPR nuclear reactor using a Serpent Monte Carlo method as distinct from the conventional deterministic computing suite. Running 3 billion neutron histories in Serpent and quantifying the Shannon entropy, the convergence of the fission source distribution is achieved in the proposed Serpent 3D full-core model. To represent the thermal-hydraulic conditions, the full-core model incorporates the water density variations with the core height as well as the other branches of the operating conditions. A source biasing method is applied for the variance reduction efficiently, with less statistical noises and a more reliable prediction of the relative core power distribution observed especially for the periphery of the core, i.e. the highest uncertainty has gone down. However, the lowest uncertainty is likely to go up with this method, which merits a further investigation for the variance optimisation. Figure-of-merit can be quantified to characterise the accuracy per unit computer time for the future efficiency study. The obtained full-core safety coefficients and power peaking factors in this work agree well in statistics with the Preliminary Safety Report [4] based on the conventional deterministic code packages [7] employed by AREVA/EDF.
Based on the Monte Carlo simulation of particle collisions, Serpent's 2D and 3D single-stage precise modelling capability without resorting to approximations is arguably an advantage over the multi-stage deterministic transport methods. Nevertheless, obtaining a reliable 3D power density distribution within an acceptable computation time remains challenging, especially for obtaining the 3D pin-wise results for the full-core Monte Carlo burnup. Sensitivity calculations and validation against the experimental data (after the core starts commissioning in the future) could be performed with ease based on the whole-core model developed in this work for multi-objective optimisations. The knowledge obtained from the EPR Monte Carlo modelling is also transferrable for simulating an AP1000 startup core, so that a comparison with the EPR's core performance and neutron economy could be performed, targeting a more reliable, low-carbon and cost-competitive nuclear power solution to approach the future energy challenges [24][25][26] in the UK.
Funding: This research received no external funding.