Monte Carlo Investigation of the UK’s First EPR Nuclear Reactor Startup Core Using Serpent

Abstract

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. For the first time, 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 state-of-the-art 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 reaction rate (power) distributions radially and axially from the three dimensional (3D) single assembly level to a 3D full core. Shannon entropy is quantified to characterise the convergence behaviour 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 and informs the decision-making by energy regulatory bodies and global partners.

1. Introduction

Maintaining electricity supplies reliably is vital during the lockdown as well as the long-term combat against the Coronavirus (COVID-19). As a low-carbon energy solution, nuclear power is projected to supply more than one-third of the UK’s electricity by 2035 [1]. Based on this vision and with the design acceptance by the UK Office for Nuclear Regulation (ONR) [2,3], EDF Energy is building up the first Evolutionary Pressurised Water Reactor (EPR) in the UK, supplied by AREVA [4]. The new generation of nuclear power station is projected to supply electricity for six million homes, whilst mitigating carbon emissions by 9 million tonnes per year for its commissioning lifespan of 60 years [5]. As distinct from existing Pressurised Water Reactors (PWR), the 1600 MWe EPR being constructed at Hinkley Point C is claimed to exhibit a large volume core with a smaller surface to volume ratio, a unique heavy reflector, improved neutron economy, lower power density, enhanced safety with a core catcher, multiple enrichments zoning, and superior fuel cycle flexibility with the mixed oxide fuel (MOX) [2,3,4,5,6,7].

Apart from the Pre-Construction Safety Report [2,3,4,5,6,7,8] submitted to the ONR, little research attention to our knowledge has been paid into an independent evaluation of the UK EPR’s full-core power distribution and safety coefficients. To be more specific, limited models are available from up-to-date literature for the whole-core optimisation, which restricts the public understanding of the discipline regarding reactor core designs and nuclear safety. To this end, this work contributes a realistic EPR reference core model using a state-of-the-art Serpent Monte Carlo code for the first time for an independent review of the analyses performed by EDF, whose results rely on a deterministic code APOLLO 2 for the 2-group macroscopic data library generation, and then a nodal code SMART [8] for the whole core calculations. The proposed Serpent code model in this work allows ongoing optimisations of the EPR compared with the recent advances in other PWR designs (e.g., AP1000, Sizewell B) to inform the regulatory authorities.

The main target of modelling a nuclear reactor core is reliably predicting 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 production and loss of neutrons are closely involved in the prediction of the reactivity and power distributions. This reflects the significance of the Boltzmann transport equation [9], an integro-differential equation in seven 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 [10] and the stochastic method [11].

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 [10]. 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 advanced Monte Carlo codes. As a 3D continuous-energy Monte Carlo burnup code, Serpent is specialised in lattice physics calculations for generating homogenised group constants [11]. Code-to-code benchmarks between Serpent 2 and other more “classic” Monte Carlo codes (e.g., MCNP 6) have been applied for the licensing and verification of other thermal reactors (e.g., the VVER-1000 mock-up [12]), while no benchmarking study has delved into the UK’s first EPR. 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 [13], especially for 3D full-core burnup calculations. Nevertheless, it is critically important 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 because of the statistical nature of the Monte Carlo approach, which can be further mitigated by simulating more neutron histories and implementing the source biasing.

2. 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.

2.1. 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,14,15,16,17,18]. The whole-core configuration with 241 assemblies is precisely represented in Figure 2 (radial view) and Figure 3 (axial view). Each assembly is represented by 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 because of 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) [14], 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.2 wt%) are loaded in the peripheral region, while multiple relatively low enrichment assemblies (2.1 wt%) are interspersed in the core interior in a checkerboard fashion.

2.2. Multi-Dimensional Branches of Diverse Operating Conditions

The calculation is firstly tasked with deriving the multiplication factor versus 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 [14,15,16,17,18].

2.3. Modelling Gadolinia Burnable Poison Pins and Heavy Reflectors

Precisely modelling the burnup characteristics of gadolinia (Gd) burnable poison pins in a reactor core 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, the highest neutron absorption cross sections are exhibited by Gd-157 and Gd-155; 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.

Based on this, Figure 7 depicts our Serpent results of burning an EPR 3.2 wt% fuel pin with 8 wt% Gd (modelled in 10 annular rings), with a special focus on the atomic density of Gd-157 versus the pin 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 behaviour 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 assemblies interfacing 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).

2.4. Convergence and Fission Source Entropy

Concerning the Monte Carlo criticality problems, the effective multiplication factor might converge sooner than the fission source distribution [19]. The convergence behaviour of Shannon entropy [19] 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. Running options setting in Serpent for convergence study.

Models	Convergence in Monte Carlo Iterations
	Neutron Histories	Number of Neutrons Per Cycle	Active Cycles	Inactive Cycles
2D assembly	10 million	20,000	500	50
3D assembly	30 million	30,000	1000	100
3D full core	3 billion	2,000,000	1500	200

according to the convergence behaviour of the multiplication factor as well as the quantified Shannon entropy as illustrated in Figure 9 for the 3D whole core at the BOL under HFP. Source convergence takes a longer full-core computational time. 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 within the first 200 inactive neutron cycles. The distribution is subsequently employed to adapt the emitted fission neutrons’ number during the active cycles. A 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.

2.5. 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 aggravated computation burden (space and time). There remains a huge potential for the variance reduction [20] at a reasonable cost of the computational burden. In our EPR model, the upper 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 [20], 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 core centre as illustrated in Figure 10b (with source biasing) to reduce the regional uncertainties and statistical noises occurring in Figure 10a (without 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.

3. Results and Analysis

3.1. 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. Modern lattice codes are reasonably expected to predict a single assembly’s multiplication factor within a few hundred percent milli-rho (pcm) of accuracy. To avoid inconsistencies, same compositions of the fuel and identical operating settings are employed in the two-codes benchmark models. Both the assemblies with and without Gd are examined with respect to burnup. The BOL calculation is featured with using smaller sub-steps of irradiation to accurately characterise the material compositions’ variations in a faster rate. Figure 11 reports the results for the 2.1 wt% assembly (8 Gd pins admixed), illustrating the larger number of discrete time steps employed at the BOC to precisely predict the Gd depletion. The benchmark here shows a good agreement.

Table 2. Serpent vs. WIMS results benchmark for EPR 2D assemblies (BOL 0.5 GWd/t, nominal power, equilibrium Xenon, 0 ppm boron, Rod Cluster Control Assembles (RCCA) out, cycle 1) in 5 uranium enrichment zones: zone 1 (4.2% + 16 Gd), zone 2 (3.2% + 20 Gd), zone 3 (3.2% + 16 Gd), zone 4 (2.1% + 8 Gd), zone 5 (2.1%, no Gd). All Gd in 8 wt%. 2-sigma uncertainties for Serpent results.

Zones	Multiplication Factor
	Serpent 2D	WIMS 2D	Deviation (pcm)
1	1.15425 (±0.00101)	1.15448	−17
2	1.03430 (±0.00121)	1.03363	62
3	1.07991 (±0.00114)	1.07942	42
4	1.06194 (±0.00109)	1.06158	31
5	1.18747 (±0.00104)	1.18748	−0.7

below quantifies the deviation of the multiplication factors between WIMS and Serpent regarding all the five types of EPR fuel assemblies burning to 0.5 GWd/t at the BOL (Xenon building up to the 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.

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. Serpent vs. PANTHER results benchmark for the EPR 3D assembly (2.1 wt%, RCCA out, boron 1000 ppm, cycle 1, nominal power). 2-sigma uncertainties for Serpent results.

Burnup Steps (GWd/t)	Multiplication Factor
	Serpent 3D	PANTHER 3D	Deviation (pcm)
0	1.08255 (±0.00014)	1.08172	71
0.1	1.07937 (±0.00014)	1.07832	90
1	1.07120 (±0.00015)	1.07007	99
10.5	0.99507 (±0.00016)	0.99410	98
21	0.92778 (±0.00017)	0.92685	108

below), hence gaining confidence for the Serpent 3D single assembly model setup. Due to the increased amount and variety of fission products produced with the fuel depletion, the statistical uncertainties in Serpent as well as the results deviation with PANTHER are largely growing.

To calculate the power distribution axially over the first cycle, we firstly perform 3D single-channel modelling of the EPR fresh core based on 21 axial meshes (layers), with thermal-hydraulic conditions incorporated. Each mesh (layer) is characterised with 0.2 m in height. Figure 12 shows a cosine-like shape of the flux (axial power profile) at the hot zero power (HZP) state due to the highly uniform axial distribution of the water density, i.e., the core’s mid-plane exhibits the flux peak, which is within expectation. Progressing towards the hot full power (HFP) condition at the BOC, the axial temperature distribution of the coolant is changing, with the outlet temperature higher than that for the inlet, the resulting thermal expansion effect of which modifies the moderator-to-fuel ratio (in the atomic number). This ratio is dropping axially from the bottom core to the top, i.e., density in the coolant is reducing (less efficient moderation) with the core height. Therefore, a hardened neutron spectrum and higher leakage occur at the region above the core’s mid-plane. The peak position of the flux (power profile) is accordingly shifting from the mid-plane downwards by 10% (in height). Due to this shifting of the higher thermal flux towards the bottom core, the fuel below the mid-plane is burning down much quicker (with more Xenon building up), giving rise to another skewing (re-distribution) of the axial power from the middle of cycle (MOC) to the end of cycle (EOC) in Figure 12, i.e., the top region of the core height exhibits the flux peaking.

3.2. 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 13a for the symmetric radial distributions of the normalised thermal flux as well as the fission rate (power). The convergence characteristics for 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 13b 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 13b. 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.

3.3. 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. EPR full-core reactivity coefficients and hot channel factors benchmark (cycle 1, HFP, RCCA out): Serpent Monte Carlo (this work) vs. deterministic codes sequence by the UK Office for Nuclear Regulation (ONR).

Safety-Related Reactivity Coefficients		Serpent in this Work	ONR [4]	Difference
Total enthalpy rise hot channel factor		1.63 (±0.01)	1.61	1.24%
Total heat flux hot channel factor		2.69 (±0.01)	2.82	−4.61%
Critical boron concentration (ppm)		1029 (±25)	1026	3 ppm
Boron worth (pcm/ppm)	BOC	−9.39 (±0.12)	−9.70	0.31 pcm/ppm
	EOC	−9.11 (±0.16)	−9.30	0.19 pcm/ppm
Doppler coefficient (pcm/K)	BOC	−2.90 (±0.05)	−2.50	−0.4 pcm/K
	EOC	−3.17 (±0.13)	−3.05	−0.12 pcm/K
Moderator temperature coefficient (pcm/K)	BOC 1	−13.71 (±0.48)	−13.00	−0.71 pcm/K
	EOC	−64.24 (±2.13)	−60.60	−3.64 pcm/K

¹ Results comparison based on HZP (for both Serpent and ONR).

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 [21], while up to 10 MeV for APOLLO based on JEF-2.2 [22]). 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 [23] and the zoning of poisoned rods [24]. Low leakage [25] in/out fuel management regime can be examined. Furthermore, the optimisation between the modelling accuracy and the computational burden can be attempted by strategically adjusting the branch points. While the thermal-hydraulic conditions have been incorporated by varying the axial temperature based on the power distribution, future work may explore different heat exchange regimes by means of coupled multi-physics simulations with computational fluid dynamics (CFD) to further justify the thermal-hydraulic feedback.

4. 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 for the fission source distribution is demonstrated 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 [8] 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 [26]. 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 whole-core Monte Carlo burnup. In response to this, hybrid deterministic/Monte Carlo and domain-decomposition methods [27] could be incorporated for future development. 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 [28], low-carbon, cost-competitive, and innovative nuclear energy system [29,30] to approach the future energy challenges [31] in the UK.