On the Feasibility of Neutron Noise Diagnostics in Fast Gen-IV Reactors Comparisons of One-Dimensional Neutron Noise Characteristics in Fast and Thermal Systems using Linear Multi-Group Diffusion Theory Master’s thesis in Physics HAMPUS HANSEN THEODOR HANSSON DEPARTMENT OP PHYSICS CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2025 www.chalmers.se www.chalmers.se Master’s thesis 2025 On the Feasibility of Neutron Noise Diagnostics in Fast Gen-IV Reactors Comparisons of One-Dimensional Neutron Noise Characteristics in Fast and Thermal Systems using Linear Multi-Group Diffusion Theory HAMPUS HANSEN THEODOR HANSSON Reactor Physics, Modelling and Safety Group Division of Subatomic, High-Energy and Plasma Physics Department of Physics Chalmers University of Technology Gothenburg, Sweden 2025 On the Feasibility of Neutron Noise Diagnostics in Fast Gen-IV Reactors Comparisons of One-Dimensional Neutron Noise Characteristics in Fast and Thermal Systems using Linear Multi-Group Diffusion Theory HAMPUS HANSEN, THEODOR HANSSON Acknowledgements, dedications, and similar personal statements in this thesis, reflect the authors’ own views. © HAMPUS HANSEN, THEODOR HANSSON, 2025. Supervisor: Christophe Demazière Examiner: Christophe Demazière Master’s thesis 2025 Department of Physics Division of Subatomic, High-Energy and Plasma Physics Reactor Physics, Modelling and Safety Group TIFX05 Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Spatial noise profiles for a BWR (top) and LFR (bottom). Each coloured line represents a different frequency, where the horizontal direction indicates the axial direction and the vertical direction indicates noise strength. Frequencies increase from right to left. For the BWR, the frequencies span the range 4.23 mHz-15.6 kHz, and the LFR covers 4.23 mHz-159 kHz. Typeset in LATEX Printed by Chalmers digitaltryck Gothenburg, Sweden 2025 ii On the Feasibility of Neutron Noise Diagnostics in Fast Gen-IV Reactors Comparisons of One-Dimensional Neutron Noise Characteristics in Fast and Thermal Systems using Linear Multi-Group Diffusion Theory HAMPUS HANSEN & THEODOR HANSSON Department of Physics Chalmers University of Technology Abstract This master’s thesis studies neutron noise, i.e. the temporal fluctuations from the mean neutron flux, in nuclear reactors. Neutron noise can be used to identify and localise faults and anomalies in such systems through a technique called neutron noise diagnostics, potentially providing early warning and leading to increased safety and operational availability. Due to the limited amount of neutron flux instrumentation available in a typical power reactor, the spatial shape of the noise is the key quantity to characterise the type of noise present. This study differentiates itself from prior research on neutron noise diagnostics, which mostly focused on Light-Water Reactors (LWR), whereas this work investigates the properties of the neutron noise in both LWRs and in Liquid Metal-cooled Fast Reactors (LMFR). More specifically, we compare the axial shape of the neutron noise of: Pressurised Water Reactors (PWR), Boiling Water Reactors (BWR), Sodium-cooled Fast Reactor (SFR) and Lead-Cooled Fast Reactors (LFR). The former two are LWRs, and the latter two are LMFRs. Each system’s deviation from point-kinetics is analysed, which is essential to evaluate the capability of noise diagnostics. To this end, a 1D numerical neutron transport solver was developed using linear multi-group diffusion theory with a variable number of delayed neutron families. The solver uses group constants generated by the Monte Carlo code Serpent-2 for each modelled system. The static neutron flux is described by an eigenvalue problem, which is solved by the power iteration method, accelerated using Wieland’s technique. The induced noise arises from an absorber of variable strength-type perturbation in the macroscopic absorption cross-section, and the resulting source-type problem is solved by matrix inversion. Our computations show that the induced neutron noise of the LWRs deviates more from point-kinetics, and their axial shape is thus more informative than that of the LMFRs. Also, perturbations closer to the reactor’s periphery induce a neutron noise that is more informative. The point-kinetic behaviour of the LFR is strongly dictated by the active height of the core, more so than by the fast spectrum. We also show that reactors fuelled by plutonium are more point-kinetic than those fuelled by uranium. For the results to become more conclusive, more work is needed regarding the types of reactors studied, with more analysis aimed towards what parameters affect the point-kinetic nature of the noise. Studies can also be made in higher spatial dimensions that more accurately take into account the complex geometries of nuclear reactors and incorporate burn-up calculations when generating the group constants. Keywords: Neutron noise, Gen-IV, Fast reactor, LMFR, Multi-group diffusion theory, Numerical meth- ods, Point-kinetics iii Acknowledgements Without the guidance and helpful advice of our supervisor, Professor Christophe Demazière, this thesis would not have been possible. Always ”leaving his door open” shows his openness to discussing both physics and any issues we might have been encountering. Even though Associate Professor Paolo Vinai has not been directly involved in this project, his company and encouragement has been appreciated. We would also like to thank Håkan T. Johansson for helping us with all things Linux, and Fredrik Öhrlund for giving us a head start with Serpent-2. To our friends, classmates, and colleagues who have made our time at Chalmers memorable. To our families, who love and support us. Hampus Hansen & Theodor Hansson, Gothenburg, June 2025 v Abbreviations and Notation List of Abbreviations BWR Boiling Water Reactor GCR Gas-Cooled Reactor Gen-IV Generation IV nuclear reactor LFR Lead-cooled Fast Reactor LMFR Liquid Metal Fast Reactor LWR Light-Water Reactor MOX Mixed Oxide fuel PK Point-Kinetic PWR Pressurised Water Reactor SFR Sodium-cooled Fast Reactor UN Uranium Nitride UO2 Uranium Dioxide ZPTF Zero-Power Transfer Function Notation Guide ω Frequency ψ Angular neutron flux ϕ Scalar neutron flux ϕ+ Adjoint scalar neutron flux δϕ Neutron noise ϕ Neutron flux vector J Neutron current ν Average neutron yield from fission v Average neutron speed D Diffusion constant Σα Macroscopic cross-section of type α δΣα Perturbation of cross-section of type α Ci Delayed precursor concentration for precursor group i λi Decay constant for precursor group i β̃i Effective delayed neutron fraction in precursor group i Λ0 Mean neutron generation time χp g Prompt fission spectrum χd g Delayed fission spectrum κf Average energy release per fission event ρ Reactivity keff Effective multiplication factor Ψ Shape function P Amplification factor δϕpk PK-component dev PK-deviation measure Gsta 0 Analytical ZPTF Gdyn 0 ZPTF through system response vii Contents Abbreviations and Notation vii List of Figures xi List of Tables xiii 1 Introduction 1 1.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Aim . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Theory 5 2.1 Neutron balance equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.1 Neutron transport theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.2 Multi-group diffusion theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.1.3 Steady-state neutron balance equations . . . . . . . . . . . . . . . . . . . . . . . . 10 2.1.4 First-order neutron noise balance equations . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Point-kinetics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.1 Shape and point-kinetic components . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2.2 Zero-power transfer function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 3 Methodology 17 3.1 Stochastic methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.1.1 Reactor Specifications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.1.2 Monte-carlo parameter estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.2 Deterministic methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2.1 Spatial discretisation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 3.2.2 Steady-state solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 3.2.3 Dynamic noise solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4 Results and Discussion 25 4.1 Model verification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 4.2 PK-deviation measure comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 4.3 Analysis of representative cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 4.4 Reactor size: differentiate spectra- from size-differences . . . . . . . . . . . . . . . . . . . . 35 5 Conclusion and Outlook 39 A Specifications A-1 A.1 SERPENT-2 micro-group specifications . . . . . . . . . . . . . . . . . . . . . . . . . . . . A-1 ix List of Figures 3.1 A horizontal representation of the BWR fuel pin. The inner most region represents the fuel followed by a helium gap and then the Zircaloy-2 cladding. Each fuel pin is surrounded by coolant. The square symmetry is also shown. . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.2 A horizontal representation of the PWR fuel pin. The innermost region represents the fuel, followed by a helium gap and then the Zircaloy-4 cladding. The square symmetry is also shown. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3 A horizontal representation of the LFR fuel pin. The innermost region represents the fuel, followed by a helium gap. The cladding consists of two different alloys, an outer corrosion-resistant and an inner structural layer. The fuel pin lattice has a hexagonal symmetry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.4 A horizontal representation of the SFR fuel pin. The innermost region represents the central void, followed by the fuel and a helium gap. The cladding consists of a 15-15Ti steel. The fuel pin has a hexagonal symmetry. . . . . . . . . . . . . . . . . . . . . . . . . . 19 3.5 Residual r of the steady-state solver as a function of the iteration number i for the two different static fluxes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.1 Comparison between the frequency-dependence of the analytical expression of the ZPTFs, denoted Gsta 0 (ω) (represented by full lines), and the ZPTF through system response (rep- resented by diamond markers). The dashed lines mark the 1/β̃ magnitude, and the dot- dashed lines mark the bounds of the plateau region, [λ, β̃/Λ0]. . . . . . . . . . . . . . . . 26 4.2 The relative deviation for each energy group, and for each system. The deviation is defined as |ϕg − ϕS2 g |/ϕS2 g , where ϕS2 g is the flux as calculated by Serpent-2; see Section 3.1.2. . . 27 4.3 Comparison between two different PK-deviation measures: the well-behaved measure dev(ω) used in this work, and the faulty groupwise integrated measure ∑Ng g devg(ω). The source term consists of a 0.1 % increase of the absorption cross-section at the lowest energy- group, localised 45 % off-centre. The noise and its components at frequency 2.78 Hz are presented in (a), and the PK-deviation measures are presented in (b). . . . . . . . . . . . 28 4.4 Deviation from point-kinetics as a function of frequency, with each system being repre- sented as: BWR - yellow full line; PWR - blue dashed line; SFR - red dot-dashed line; LFR - green dotted line. The shaded regions represent the energy-group dependence of the source term for each system. Spatially, the left subfigure has a centred source term axially, with zp = 0, and the right subfigure has an off-centred source term, with zp = −0.35H. . . 30 4.5 Normalised static scalar flux, presented both groupwise, ϕg,0(z) in (a), and total, ϕ0(z) in (b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.6 Axial shape of the neutron noise induced by a centred perturbation zp = 0 of energy group g = 7 at frequency ω = 2.78 Hz. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.7 Axial shape of the neutron noise induced by an off-centred perturbation zp = −0.35H of energy group g = 7 at frequency ω = 2.78 Hz. . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.8 Axial shape of the total neutron noise, for four different frequencies, represented (in in- creasing order) as: blue, red, yellow and purple. . . . . . . . . . . . . . . . . . . . . . . . . 34 4.9 Axial shape of the different components of the noise, summed over all energy-groups, at frequency ω = 2.78 Hz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 xi List of Figures List of Figures 4.10 Deviation from point-kinetics as a function of frequency, with same-sized systems, with each system being represented as: BWR - yellow full line; PWR∗ - blue dashed line; SFR∗ - red dot-dashed line; LFR∗ - green dotted line; SFR∗ U - purple dot-dashed line. The shaded regions represent the energy-group dependence of the source-term for each system. Spatially, the left subfigure has a centred source term axially, with zp = 0, and the right subfigure has an off-centred source term, with zp = −0.35H. . . . . . . . . . . . . . . . . . 36 4.11 Axial shape of the neutron noise at frequency ω = 2.78 Hz for the height-modified models, including the fuel-modified model SFR∗ U. The noise is induced by an off-centred perturba- tion at zp = −0.35H in energy group g = 7. . . . . . . . . . . . . . . . . . . . . . . . . . . 37 xii List of Tables 3.1 Boiling water reactor fuel specifications. Based on SVEA Optima 8 × 8 BWR fuel assembly. 18 3.2 PWR specifications. Based on Westinghouse’s 17×17 fuel design. . . . . . . . . . . . . . . 18 3.3 Lead-cooled fuel pin specifications. Based on a theoretical 80 MWt design. . . . . . . . . . 19 3.4 Sodium-cooled fuel pin specifications. Based on the Superphénix reactor. . . . . . . . . . 19 3.5 Isotopic composition of plutonium in spent GCR nuclear fuel [28]. . . . . . . . . . . . . . 20 3.6 Macro-group energy boundaries of the thermal systems. The lower bound of energy group g = 7 is 0 eV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.7 Macro-group energy boundaries of the thermal systems. The lower bound of energy group g = 8 is 0 eV. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.8 Spatially homogenised parameters calculated using Serpent-2. Subscript g indicates that the parameter is dependent on the macro-group and subscript i is the delayed neutron family number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 3.9 Multiplication factor k∞ and effective delayed neutron fraction β̃ for our four reactors as simulated by Serpent-2 based on technical data from Section 3.1.1. . . . . . . . . . . . . 21 4.1 Characteristic amplitudes and frequencies of the ZPTF, where λ corresponds to the upper plateau frequency, β̃/Λ0 to the lower plateau frequency, and 1/β̃ is the plateau amplitude. 25 A.1 Tabulated index-value pairs in compact format for the Serpent-2 energy micro group. The values in the table indicate the group boundaries, which is the reason that there are 71 total entries for 70 total groups. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A-1 xiii 1 Introduction Neutron noise, or the temporal fluctuations in the neutron flux, is a quantity that reflects the dynamics in the core of a fission reactor. Small perturbations, which could be caused by anomalies or mechanical deteriorations, affect the neutron noise through small changes in fission, capture and scattering cross- sections, and consequently propagate throughout the core due to the transport of neutrons. The neutron noise of the system globally reflects local changes of reactor properties, such as temperature, flow rate and pressure. Using this unique property of the neutron noise to monitor the core with neutron detectors is referred to as noise diagnostics. Noise diagnostics is an analysis technique that can be used to detect, classify and sometimes localise anomalies in Light Water nuclear Reactors (LWR). This is a promising technique that can increase the uptime and safety of nuclear reactors [1], both of utmost relevance. It has, for instance, been used to analyse the coolant flow velocity in a reactor [2], and could also be used to measure various mechanical vibrations [3] or boiling patterns [4]. However, the technique has not been extensively studied for Liquid Metal-cooled Fast Reactors (LMFR), which are part of the Gen-IV fleet of nuclear reactors. This is mainly since these latter systems are uncommon, currently at a relatively low technological readiness level (TRL) [5]. It is hypothesised that performing noise diagnostics for these LMFRs is more challenging because of the fast neutron spectra and the longer diffusion lengths associated with these systems. This likely results in a less informative axial shape of the neutron flux, a concern due to the scarce instrumentation available in the core. However, this hypothesis has not been tested before, and the extent of which the fast spectra affects the spatial information carried by the noise is unknown. For this purpose, the authors have developed a numerical solver to analyse the noise properties in both LWRs and LMFRs, with respect to the neutron noise’s frequency and axial distributions. This preliminary study will not only provide an insight to what extent the fast spectra reduces the information carried by the spatial shape of the neutron noise, but also outline measures that need to be taken in order to implement noise diagnostic techniques to the fast reactors of the future, consequently improving their operational availability. 1.1 Background The splitting of heavy nuclei by neutron absorption, fission, is a process that releases large amounts of energy and between 2-3 high energy neutrons. These neutrons can split other fissile atoms, producing more neutrons, which split even more fissile atoms. Nuclear fission reactors work on this principle. The reactor is critical when the chain reaction is controlled such that the neutron population is constant from one generation to the next. This is achieved when the production of neutrons through fission is fully compensated by the loss of neutrons through capture and leakage. In 2023, 4 % of the world’s energy was produced by nuclear power, with the fraction being 20 % in Sweden, 35 % in France, which is the highest in the world [6]. The vast majority of the current operational nuclear reactor fleet consists of LWRs. However, in the future reactor fleet—consisting of Gen-IV reactors—other types of reactors, such as LMFRs, are expected to play a larger role. Reactors of this class generally aim to improve safety, reliability and operating economics. LMFRs make up two out of the three types of fast Gen-IV reactors. More specifically, these are the Lead- cooled Fast Reactor (LFR) and Sodium-cooled Fast Reactor (SFR). LMFRs have been built and tested in the past, examples of which are the Lead-Bismuth eutectic reactors used in Soviet naval reactors [7], and the full-scale sodium reactors built in France [8] and Japan [9]. All of the reactors from the previous 1 1. Introduction 1.2. Aim examples have since been shut down, but development of new reactors is underway. Examples of new LFR projects are the Swedish Blykalla, Italian ALFRED and Newcleo [7], while American TerraPower and Chinese CFR-600 are examples of up-and-coming SFR projects [10]. The most common LWRs are Boiling Water Reactors (BWR) and Pressurised Water Reactors (PWR). These are thermal reactors, since they moderate the neutrons, i.e. slowing them down from high energies after emission, approximately 2 MeV, to thermal energies, approximately 25 meV through scattering re- actions with the hydrogen in the water. The primary purpose of slowing down neutrons is that the ratio of fission to absorption of 235U is greater at lower energies, and thus allows a lower enrichment of fissile materials. In contrast, the aforementioned LMFRs are fast reactors, as they lack significant moderation. Higher levels of enrichment are thus required to attain criticality in such systems. A benefit, however, is that 238U can undergo fission at higher neutron energies, and at higher neutron incident energies, more neutrons are released per fission event [11, Ch. 1]. This increase in available neutrons allows a significantly increased fuel economy from the increased transmutation of 238U into 239Pu, which is fissile. This allows the extension of Earth’s uranium reserve by an approximate factor of 60–80 [11, Ch. 1]. The ”volume, toxicity and lifespan of the longest-living radioactive waste” is also significantly reduced in comparison to waste from thermal reactors [12]. Lower system pressures and higher coolant boiling points increase safety, and higher operating temperatures allow a higher thermodynamic efficiency. However, drawbacks exist as well. The opaque nature of metals hinders visual inspections that are possible in LWRs, and higher temperatures accelerate corrosion-induced embrittlement [13]. Neutrons with higher energies generally have longer diffusion lengths, and since LFMRs are fast systems, their neutrons in such systems will have a longer diffusion length on average. This larger diffusion length of a fast spectrum compared to a thermal spectrum means that the neutrons travel longer distances before interacting significantly with the material. This is complicated by the ever-shifting balance between production of new neutrons through fission events and removal of neutrons through absorption and leakage. A simplified but illustrative picture is now presented to emphasise the importance of the shape of the noise: let us consider the appearance of an anomaly at a point in the reactor, and that this anomaly effectively acts as a perturbation to the absorption cross-section. Longer diffusion lengths will make the neutrons interact less with the perturbed point, thus decreasing the effect the perturbation has on the neutron flux. The neutron flux, specifically its temporal fluctuations around its mean (neutron noise), will be more homogenous in its spatial shape, resembling the static neutron flux with respect to its spatial dependence. So, a fast spectrum leads to a spatial shape of the neutron noise that contains less information about the perturbation compared to a thermal spectrum. This loss of information carried by the neutron noise can be problematic for the detection and classification of the anomaly, and localisation of the anomaly becomes almost impossible with if the spatial shape of the noise is not informative enough [1]. These are paramount goals of core monitoring. The importance of the information carried by the noise is amplified by the sparsity of instrumentation available in a nuclear power reactor. Western LWRs are typically fitted with neutron detectors outside and within the core, as well as temperature detectors [2, p. 1699]. The temperature detectors are, however, not commonly mounted within the core. Even so, a local anomaly that causes a temperature disturbance would only propagate downstream from the location of the disturbance [1]. Neutron noise, on the other hand, does propagate throughout the entire reactor. A localised anomaly can therefore be measured globally by all neutron detectors —allowing its type to be classified and its location to be inferred. 1.2 Aim The promising outlook of noise diagnostics as a core monitoring technique, strengthened by the fact that it has been demonstrated to work in LWRs, makes the study of the feasibility of neutron diagnostics in fast systems a natural and important aim. If noise diagnostics is a possible core monitoring technique in fast Gen-IV systems, its implementation will help maximize the potential of the future’s reactors by increasing their operational availability and safety, while decreasing the cost [1]. This thesis aims to investigate how the spatial shape of the neutron noise differs depending on the type of reactor. This property is quantified by measuring the responses’ deviations from point-kinetics (see Section 2.2). The following questions are of special interest: 1. To which extent does the fast spectra of the LMFRs affect the axial shape of the neutron noise? 2 1.3. Scope 1. Introduction 2. How does the deviation from point-kinetics vary for the different systems, depending on frequency and energy group of the perturbation? 3. How feasible is it to use neutron diagnosis for fast reactors? The lack of prior research conducted for this technique applied to fast systems specifically means that this work acts as a preliminary study. It is preliminary in the sense that the complexity of the problem will be simplified drastically, see Section 1.3. Consequently, the modelling and analysis must be made in a more sophisticated manner and the scope needs to be broadened in order to make more definite statements. However, this study can not only start the development of solvers and methods to be used in future research, but also grant insights into which relevant questions that need to be asked and answered, which effects need to be studied and what tools and method that need to be developed in order to advance—hopefully limiting the scope of future studies. 1.3 Scope The exploratory nature of this work and its role as a preliminary study necessitate a reduction in com- plexity and scope, which in turn limits the depth of the discussion and the conclusions that can be drawn. Perhaps the largest simplification is that a full-core modelling will not be conducted, but instead, a representative unit fuel cell with surrounding coolant is instead modelled. This is then extended to simulate a fuel pin in the axial dimension. The nuclear fuel in each rod is assumed to be fresh, and no burn-up is considered. The reactor is assumed to be homogeneous, so the group constants have no spatial dependence. The neutron noise is only solved for the axial dimension. The goals of this work are, in most regards, achievable for these simplified homogenous and static one-dimensional reactor models. In this work, the thermal reactor systems considered are a PWR and a BWR. They are compared with two types of fast Gen-IV reactor systems, an SFR and an LFR. However, there exists more types of fast and thermal reactor systems which could be considered. Furthermore, only one specific unit cell of each respective type of reactor system is modelled. These unit cell are specifically modelled after a SVEA Optima 8x8 BWR fuel assembly [14, Table A-3], a Westinghouse 17×17 PWR fuel assembly [14, Table A-4], a LFR-system proposed by SUNRISE [13] and the Superphénix SFR [15, Tab. 1-2]. The results for each type of model may therefore not be representative of all reactor systems of that type. The physical anomalies of the system, the inducers of the neutron noise, will be modelled as a localised perturbation of the macroscopic cross-sections, referred to as a absorber of variable strength. Other perturbations are possible, for example axially travelling perturbations, fuel assembly vibrations, control rod vibrations and core barrel vibrations. These are more complex and many require more sophisticated models. Moreover, they are all fundamentally based on absorber of variable strength-type perturbations, making this a natural starting point for the purpose of this study. Furthermore, only perturbations of strength 0.1 % applied to the macroscopic absorption cross-section are considered. 3 2 Theory The theory of nuclear reactors is a complex subject with many different branches within physics and engineering. The curious reader is directed towards introductory books in the field, such as Lewis [16] or books in reactor analysis by Bell and Glasstone [17], or Demazière [18]. For theory specific to fast neutron reactors, see Waltar and Reynolds [11]. From this point on, we assume that the reader is familiar with the subatomic nuclear interactions that occur in nuclear reactors, as well as the definitions of key quantities such as macroscopic cross-sections. This chapter starts from the governing neutron balance equation, namely the neutron transport equation. The transition from transport theory to multi-group diffusion theory is then presented. This reduces the complexity of the equations, and the ability to numerically calculate the flux is greatly improved. The steady-state and dynamical balance equations are then derived from the multi-group diffusion equations. Some methods for evaluating the properties of the neutron noise are then discussed, such as point-kinetics and zero-power transfer functions. 2.1 Neutron balance equations In nuclear reactors, there are several processes that create, remove or displace neutrons, which can be described by balance equations. To understand them, some critical parameters are defined at position r, direction Ω, energy E and time t: Angular neutron flux: ψ(r,Ω, E, t) ≡ v(E)n(r,Ω, E, t) (2.1) Scalar neutron flux: ϕ(r,E, t) ≡ ∫ 4π ψ(r,Ω, E, t)d2Ω (2.2) Neutron current density vector: J(r,E,t) ≡ ∫ 4π Ωψ(r,Ω, E, t)d2Ω (2.3) Effective delayed neutron fraction for precursor family i: β̃i(r) = ∫ V ∫ ∞ 0 βi(r, E)ν(r, E)Σf (r, E)ϕ(r, E)dEd3r∫ V ∫ ∞ 0 ν(r, E)Σf (r, E)ϕ(r, E)dEd3r (2.4) Effective delayed neutron fraction: β̃(r) = Nd∑ i β̃i (2.5) The angular flux ψ(r,Ω, E, t) represents the neutron flux at a specific energy E and direction Ω. It is defined as the speed of the neutrons v(E) at energy E multiplied with the neutron density n(r,Ω,E,t). Meanwhile, the scalar flux is the angular neutron flux integrated over all directions. Both of these quantities are space, time and energy dependent. The net neutron flux at a point is the neutron current 5 2. Theory 2.1. Neutron balance equations J(r,E,t). The delayed neutron fraction β(r,E) represents the total fraction of neutrons that are emitted from fission fragments after they decay above the virtual energy, summed over all precursor families i and averaged over its energy dependence. Averaging β(r) over its energy dependence yields the effective delayed neutron fraction, denoted by β̃(r), which is used in this work. The quantity ν(r, E)Σf (r,E,t) is the space- and time-dependent product between the average number of neutrons released by fission event and the macroscopic fission cross-section at energy E. Additional variables are then introduced, beginning with those that depend on time and space: Σs(r,Ω′ → Ω, E′ → E, t) is macroscopic scattering cross-section from energy E′ and solid angle Ω′ to E and Ω; Σt(r, E, t) is the macroscopic total cross-section at energy E for all nuclear reactions, specifically scattering, fission and capture reactions; Ci(r, t) represents concentration of the precursors of delayed neutrons belonging to the group of precursors i. The following variables are assumed to be time independent: λi(r) is the decay constant that corresponds to each delayed neutron precursor group; the two fission spectra, prompt χp(r, E) and delayed χd i (r,E), represent what fraction of neutrons are created at energy E as prompt and delayed neutrons respectively, and from which precursor group i it originated from for the latter; v(E) is the speed of a neutron at energy E [18], [19], [17]. 2.1.1 Neutron transport theory The foundation on which every subsequent section of the theory is build upon is the balance equation, sometimes called the Boltzmann or neutron transport equation. It is defined as: 1 v(E) ∂ψ ∂t (r,Ω, E, t) + Ω · ∇ψ(r,Ω, E, t) + Σt(r, E, t)ψ(r,Ω, E, t) = ∫ 4π ∫ ∞ 0 Σs(r,Ω′ → Ω, E′ → E, t)ψ(r,Ω′, E′, t) d2Ω′dE′ + χp(r,E) 4π [ 1 − β̃(r) ] ∫ 4π ∫ ∞ 0 ν(r,E′)Σf (r, E′, t)ψ(r,Ω′, E′,t) d2Ω′dE′ + Nd∑ i=1 χd i (r,E) 4π λi(r)Ci(r, t), (2.6) where ∂Ci ∂t (r, t) = β̃i(r) ∫ ∞ 0 ν(r,E)Σf (r, E, t) ∫ 4π ψ(r,Ω′, E, t) d2Ω′dE − λi(r)Ci(r, t), i = 1, . . . , Nd, . (2.7) Each term of Eqs. (2.6) and (2.7) describe either the consumption of neutrons, the production of neutrons or the time variation of the number of neutrons. The equations can thus be regarded as balance equations, where the physical meaning of the terms are as follows: • The first term on the left-hand side of Eq. (2.6) 1 v(E) ∂ψ ∂t (r,Ω, E, t) (2.8) corresponds to the time variation of the number of neutrons at position r, direction Ω, energy E and time t. It represents the imbalance of production and consumption of neutrons at the specific configuration. • The second term on the left-hand side of Eq. (2.6) Ω · ∇ψ(r,Ω, E, t) (2.9) is the disappearance of neutrons from position r with regards to the spatial variation, referred to as the leakage term. • The third term on the left-hand side of Eq. (2.6) Σt(r, E, t)ψ(r,Ω, E, t) (2.10) is the consumption of neutrons at position r, direction Ω, energy E and time t with regard to the total nuclear reactions, including scattering, fission and capture reactions. 6 2.1. Neutron balance equations 2. Theory • The first term on the right-hand side of Eq. (2.6)∫ 4π ∫ ∞ 0 Σs(r,Ω′ → Ω, E′ → E, t)ψ(r,Ω′, E′, t) d2Ω′dE′ (2.11) is the production of neutrons with direction Ω and energy E from all other directions and energies through scattering reactions at position r and time t. • The second term on the right-hand side of Eq. (2.6) χp(r,E) 4π [ 1 − β̃(r) ] ∫ ∞ 0 ν(r,E′)Σf (r, E′, t)ϕ(r, E′,t) dE′ (2.12) is the production of neutrons at position r, direction Ω, energy E and time t due to fission reac- tions from prompt neutrons. The factor 1 4π accounts for the isotropic emission of neutrons in the laboratory reference system from fission reactions, Σf (r,Ω′ → Ω, E′ → E, t) = Σf (r, E′ → E, t) 4π . (2.13) The prompt fission spectrum factor χp(r,E) account for the distribution of the energy of the neu- trons released by (prompt) fission events, so Σf (r, E′ → E, t) = χp(r,E)Σf (r, E′, t). (2.14) • The third term on the right-hand side of Eq. (2.6) Nd∑ i=1 χd i (r,E) 4π λi(r)Ci(r, t) (2.15) is the decay of the precursors of delayed neutrons, which corresponds to the production of neutrons at position r, solid angle Ω, energy E and time t due to fission reactions from the delayed neutrons. These neutrons are also released isotropically, and their energies are distributed according to the delayed neutron spectrum factor χd i (r,E). • The left-hand side of Eq. (2.7) ∂Ci ∂t (r, t) (2.16) is the time variation of the concentration of precursors, which corresponds to the imbalance between production and decay of precursors. • The first term on the right-hand side of Eq. (2.7) β̃i(r) ∫ ∞ 0 ν(r,E)Σf (r, E, t)ϕ(r, E, t) dE (2.17) is the production of precursors of delayed neutrons at position r and time t through fission reactions. • The second term on the right-hand side of Eq. (2.7) −λi(r)Ci(r, t) (2.18) is the loss of precursors at position r and time t through radioactive decay. These terms make up the integro-differential neutron transport equation, which is very challenging to solve. Some simplifications are made to reduce the complexity and degrees of freedom of the problem. This will reduce the problem from solving the neutron transport equation to solving the far simpler multi- group diffusion equation. These simplifications are: only resolving lower moments of angular neutron flux, applying the diffusion approximation, and finally employing a multi-group formalism. These steps are explained in more depth in the following section. 7 2. Theory 2.1. Neutron balance equations 2.1.2 Multi-group diffusion theory Lower order moments of scalar flux To change the framework of the neutron balance equations from transport theory to multi-group diffusion theory we start with resolving lower moments of the angular neutron flux ψ(r,Ω,E,t), specifically the scalar neutron flux ϕ(r,E,t) and the neutron current density vector J(r,E,t). This is done by integrating Eq. (2.6) over the unit sphere. Performing this integration yields, 1 v(E) ∂ϕ ∂t (r, E, t) + ∫ 4π Ω · ∇ψ(r,Ω, E, t) d2Ω + Σt(r, E, t)ϕ(r, E, t) = ∫ 4π ∫ 4π ∫ ∞ 0 Σs(r,Ω′ → Ω, E′ → E, t)ψ(r,Ω′, E′, t) dE′d2Ω′d2Ω + χp(r,E) [ 1 − β̃(r) ] ∫ ∞ 0 ν(r,E′)Σf (r, E′, t)ϕ(r, E′,t) dE′ + Nd∑ i=1 χd i (r,E)λi(r)Ci(r, t), (2.19) where, ∂Ci ∂t (r, t) = β̃i(r) ∫ ∞ 0 ν(r,E)Σf (r, E, t)ϕ(r, E, t) dE − λi(r)Ci(r, t), i = 1, . . . , Nd. (2.20) The factors 1 4π are cancelled by integrating over the unit sphere. The leakage term in Eq. (2.19) can be rewritten as ∫ 4π Ω · ∇ψ(r,Ω, E, t) d2Ω = ∇ · ∫ 4π Ωψ(r,Ω, E, t) d2Ω = ∇ · J(r, E, t), (2.21) where Eq. (2.3) was used in the last equality. Typically, reactor simulation tools and programs only resolve the lower moments of the angular neutron flux. This is in part due to the large decrease in computational cost, but also because neutron detectors cannot measure the angular neutron flux, only the scalar neutron flux, so resolving the former in simulations or computations carries little value [18]. Multi-group formalism The multi-group formalism is a treatment of the energy dependence of the system. It consists of separating the neutron energy continuum into discrete energy bins, referred to as (macro-) energy groups, ordered in decreasing neutron energy [Emin,Emax] = 1⋃ g=Ng [Eg,Eg−1], (2.22) then integrating Eqs. (2.19) and (2.20) over each energy group {g}NG g=1, where Ng is the total number of energy groups. The bounds of each energy group for both the LWRs and the LMFRs are presented in Table 3.6-3.7 in Section 3.1.2. This further reduces the complexity of the transport equation. To maintain the reaction rates when the transport equation is integrated over each energy bin, the following methodology is applied to all energy dependent quantities: ϕg(r,t), ψg(r,Ω,t), χp g(r) and χd i,g(r) are simply integrated over each energy group as presented here for the neutron current: Jg(r,t) = ∫ Eg−1 Eg J(r,E,t) dE. (2.23) For the cross-sections, the metrics are weighted according to the flux. This is the case for 1 vg , Σt,g(r,t), Σs,g′→g(r,Ω′ → Ω,t), all of which are normalized as presented here for the macroscopic fission cross- section: (νΣf )g (r,t) = ∫ Eg−1 Eg ν(r,E)Σf (r,E,t)ϕw(E) dE∫ Eg−1 Eg ϕw(E) dE , (2.24) where ϕw(E) is a typical flux. The specifics of this flux-weighting is described in [20]. 8 2.1. Neutron balance equations 2. Theory Performing the integrations in Eqs. (3.1)-(3.3) on all energy-dependent quantities in Eqs. (2.19) and (2.20) now yields the multi-group transport equations integrated on all solid angles, 1 vg ∂ϕg ∂t (r, t) + ∇ · Jg(r, t) + Σt,g(r, t)ϕg(r, t) = Ng∑ g′=1 ∫ 4π Σs,g′→g(r,Ω′ → Ω, t)ψg′(r,Ω′, t) d2Ω′ + χp g(r) [ 1 − β̃(r) ] Ng∑ g′=1 (νΣf )g′ (r, t)ϕg′(r, t) + Nd∑ i=1 χd i,g(r)λi(r)Ci(r, t), (2.25) and ∂Ci ∂t (r, t) = β̃i(r) Ng∑ g=1 (νΣf )g (r, t)ϕg(r, t) − λi(r)Ci(r, t), i = 1, . . . , Nd. (2.26) Diffusion approximation The last simplification applied to Eq. (2.25) and (2.26) is the transition to diffusion theory through the diffusion approximation, which resolves the remaining angular dependence. This transition is formally derived through spherical harmonics and the associated Legendre polynomials, through the so-called PN - method, the details of which are addressed in [18, Ch. 4]. The idea is to consider the time-independent multi-group transport equation, expand the angular dependence of the angular neutron flux and the angular dependence of the scattering matrix onto spherical harmonics and the Legendre polynomials, then truncating at first order. After some more work, the P1-equations are acquired ∇ · Jg(r) + Σt,gϕg(r) = Ng∑ g′=1 Σs0,g′→g(r)ϕg′(r) + χg(r) k Ng∑ g′=1 (νΣf )g′(r)ϕg′(r), (2.27) 1 3∇ϕg(r) + Σt,g(r)Jg(r) = Ng∑ g′=1 Σs1,g′→g(r)Jg′(r), (2.28) where Σs0,g′→g(r) and Σs1,g′→g(r) are the first two Legendre moments of the scattering cross-section defined as Σs0,g′→g(r) = 2π ∫ 1 −1 Σs,g′→g(r,µ) dµ (2.29) Σs1,g′→g(r) = 2π ∫ 1 −1 Σs,g′→g(r,µ)µdµ, (2.30) where µ ≡ Ω · Ω′ (2.31) that represents the cosine of the deviation angle between incoming and outgoing neutrons. Lastly, we define the following group-wise fission spectra, χg(r) ≡ χp g(r)[1 − β̃] + Nd∑ i=1 χi,g(r)β̃i(r). (2.32) Note that the zeroth order Legendre moment Σs0,g′→g(r) is isotropic in the laboratory reference system, but the first order Legendre moment Σs1,g′→g(r) is anisotropic. The central assumption in the diffusion approximation is that the anisotropic in-scatter into group g is exactly equal to the anisotropic out-scatter from g, represented by Ng∑ g′=1 Σs1,g′→g(r)Jg′(r) = Ng∑ g′=1 Σs1,g→g′(r)Jg(r). (2.33) This means that the anisotropic in-scatter to some energy group is countered by anisotropic out-scattering from that group to any other, so the neutron current will be driven solely by neutron flux gradients. As a 9 2. Theory 2.1. Neutron balance equations consequence of this assumption, the neutron current vector is proportional to the divergence of the scalar neutron flux, seen by inserting Eq. (2.33) into the Eq. (2.28). This proportionality relation is known as Fick’s Law, 1 3∇ϕg(r) + Σt,g(r)Jg(r) = Ng∑ g′=1 Σs1,g→g′(r)Jg(r) ⇐⇒ Jg(r) = −Dg(r)∇ϕg(r), (2.34) where the proportionality constants define the diffusion coefficients, Dg(r) = 1 3 [ Σt,g(r) − ∑Ng g′=1 Σs1,g→g′(r) ] ≡ 1 3Σ0 t,g(r) , (2.35) where Σ0 t,g(r) is known as the transport corrected total cross-section. This quantity is considered trans- port corrected because it includes the Σs1-term, and thus includes most of the anisotropic behaviour associated with the angular dependence. However, since the derivation above is made through the P1- method, where transport corrected cross-sections cancel out in Eq. (2.27), a discussion about transport correction has little importance. One must only ensure that both Σt,g(r) and Σs0,g′→g(r) in Eq. (2.27), which are generated through external simulations, are either transport corrected or not. Now, generalizing this to the time dependent problem, using the diffusion approximation, expressing the leakage term with the diffusion coefficients and expanding χg(r), we end up with the multi-group diffusion equation 1 vg ∂ϕg ∂t (r, t) =∇ · [Dg(r, t)∇ϕg(r, t)] + Ng∑ g′=1 Σs0,g′→g(r, t)ϕg′(r, t) + [ 1 − β̃(r) ] χp g(r) Ng∑ g′=1 νg′(r)Σf,g′(r, t)ϕg′(r, t) + Nd∑ i=1 χd g,i(r)λi(r)Ci(r, t) − Σt,g(r, t)ϕg(r, t), (2.36) for the neutron flux and the balance equation for each of the groups of precursors of delayed neutrons ∂Ci ∂t (r, t) = β̃i(r) Ng∑ g′=1 νg′(r)Σf,g′(r, t)ϕg′(r, t) − λi(r)Ci(r, t), i = 1, ..., Nd. (2.37) 2.1.3 Steady-state neutron balance equations All time-dependent quantities X(r,t) can be separated into X(r,t) = X0(r) + δX(r,t), (2.38) where X0(r) is the time average of the corresponding quantity and δX(r,t) is the time-dependent fluc- tuations around the time-average. If the scalar neutron flux is considered as the general quantity X(r,t), then the fluctuations are referred to as the neutron noise. The fluctuations are assumed to be small in comparison to the mean valued quantities, |δX(r,t)| << |X0(r)| ∀(r,t), (2.39) and static, such that the time average is zero,〈 δX(r,t) 〉 = 0 ∀ (r,t) =⇒ 〈 ∂δX ∂t (r,t) 〉 = 0 ∀ (r,t). (2.40) The neutron noise is characterised by the small, time-dependent fluctuations of the scalar neutron flux. In the following two sections, balance equations of the neutron noise are derived by separating the multi- group diffusion equations, Eqs. (2.36) and (2.37), into one set of balance equations for the steady-state (presented in this section) and one set of balance equations for the dynamical, time-dependent case 10 2.1. Neutron balance equations 2. Theory (presented in Section 2.1.4). The solution to the former set of equations yields the static flux ϕg,0(r) and the latter yields the neutron noise, δϕg(r,ω). The steady-state set of equations must be solved to compute the neutron noise. The solutions to these sets of equations are described in Section 3. Separating the instantaneous quantities X(r,t) in Eqs. (2.36) and (2.37) into their time-averaged mean X0(r) and their time-dependent fluctuations δX(r,t) yields 1 vg ∂ ∂t [ ����ϕg,0(r) + δϕg(r, t) ] =∇ · [ Dg(r, t)∇ [ ϕg,0(r) + δϕg(r, t) ]] + Ng∑ g′=1 [ Σs0,g′→g,0(r) + δΣs0,g′→g(r, t) ][ ϕg′,0(r) + δϕg′(r, t) ] + [ 1 − β̃(r) ] χp g(r) Ng∑ g′=1 νg′(r) [ Σf,g′,0(r) + δΣf,g′(r, t) ][ ϕg′,0(r) + δϕg′(r, t) ] + Nd∑ i=1 χd g,i(r)λi(r) [ Ci,0(r) + δCi(r, t) ] − [ Σt,g,0(r) + δΣt,g(r, t) ][ ϕg,0(r) + δϕg(r, t) ] , (2.41) and ∂ ∂t [ ����Ci,0(r) + δCi(r, t) ] =β̃i(r) Ng∑ g′=1 νg′(r) [ Σf,g′,0(r) + δΣf,g′(r, t) ][ ϕg′,0(r) + δϕg′(r, t) ] − λi(r) [ Ci,0(r) + δCi(r, t) ] , (2.42) where we have assumed that the fluctuations of the diffusion coefficient are negligible, such that Dg(r, t) ≈ Dg,0(r) ∀ (r,t), (2.43) which is a valid assumption for LWR systems [21], but also assumed here to be valid for fast systems. Whether this assumption holds is not certain and should preferably be subject to further studies. The steady-state balance equations are time-independent by definition. Assuming no time dependence of the quantities in Eqs. (2.41) and (2.42), or equivalently only regarding the mean values and using Eq. (2.40), we get the steady-state balance equations, ∇ · [Dg,0(r)∇ϕg,0(r)] + Σt,g,0(r)ϕg,0(r) = Ng∑ g′=1 Σs0,g′→g,0(r)ϕg′,0(r) + [ 1 − β̃(r) ] χp g(r) Ng∑ g′=1 νg′(r)Σf,g′,0(r)ϕg′,0(r) + Nd∑ i=1 χd g,i(r)λi(r)Ci,0(r), (2.44) and β̃i(r) Ng∑ g′=1 νg′(r)Σf,g′,0(r)ϕg′,0(r) − λi(r)Ci,0(r) = 0, i = 1, ..., Nd, (2.45) which are satisfied by definition. The equations above are slightly idealised, since we have assumed a steady-state. It is highly unlikely that the system will be exactly critical for a given set of macroscopic cross-sections. Further deviation from exact criticality are expected from the discretisation of the system, which is required to solve the equations numerically. However, since there exists a configuration at each time step where the system is exactly in steady-state conditions, it is possible to maintain the steady- state by renormalizing the macroscopic fission cross-section with the effective multiplication factor keff [18]. This quantity is defined as the neutron production from fission in one generation divided by the 11 2. Theory 2.1. Neutron balance equations neutron absorption and leakage in the preceding generation, and it is equal to one only when the system is critical. This results in the final time-independent steady-state equations: ∇ · [Dg,0(r)∇ϕg,0(r)] + Σt,g,0(r)ϕg,0(r) = Ng∑ g′=1 Σs0,g′→g,0(r)ϕg′,0(r) + χp g(r) keff [ 1 − β̃(r) ] Ng∑ g′=1 νg′(r)Σf,g′,0(r)ϕg′,0(r) + Nd∑ i=1 χd g,i(r)λi(r)Ci,0(r), (2.46) and β̃i(r) keff Ng∑ g′=1 νg′(r)Σf,g′,0(r)ϕg′,0(r) − λi(r)Ci,0(r) = 0, i = 1, ..., Nd. (2.47) This steady-state set of equations is interpretable as an eigenvalue problem, where we identify ϕg,0(r) as the eigenvector and keff as its corresponding eigenvalue. The method of solving this problem is described in Section 3.2.2. 2.1.4 First-order neutron noise balance equations To retrieve the dynamical equations we expand the binomial products in Eqs. (2.41) and (2.42), then renormalize the macroscopic fission cross-section with keff. We then subtract the steady-state Eqs. (2.46) and (2.47) from Eqs. (2.41) and (2.42) respectively, yielding 1 vg ∂ ∂t δϕg(r, t) =∇ · [Dg,0(r)∇δϕg(r, t)] + Ng∑ g′=1 Σs0,g′→g,0(r)δϕg′(r, t) + Ng∑ g′=1 δΣs0,g′→g(r, t)ϕg′,0(r) + Ng∑ g′=1�����������: O(δ2) δΣs0,g′→g(r, t)δϕg′(r,t) + χp g(r) keff [ 1 − β̃(r) ] Ng∑ g′=1 νg′(r)Σf,g′,0(r)δϕg′(r, t) + χp g(r) keff [ 1 − β̃(r) ] Ng∑ g′=1 νg′(r)δΣf,g′(r, t)ϕg′,0(r) + χp g(r) keff [ 1 − β̃(r) ] Ng∑ g′=1������������:O(δ2) νg′(r)δΣf,g′(r, t)δϕg′(r, t) + Nd∑ i=1 χd g,i(r)λi(r)δCi(r, t) − Σt,g,0(r)δϕg(r, t) − δΣt,g(r, t)ϕg,0(r) − ���������: O(δ2) δΣt,g(r, t)δϕg(r, t) , (2.48) and ∂ ∂t δCi(r, t) = β̃i(r) keff Ng∑ g′=1 νg′(r)Σf,g′,0(r)δϕg′(r, t) + β̃i(r) keff Ng∑ g′=1 νg′(r)δΣf,g′(r, t)ϕg′,0(r) + β̃i(r) keff Ng∑ g′=1������������:O(δ2) νg′(r)δΣf,g′(r, t)δϕg′(r, t) − λi(r)δCi(r, t). (2.49) The terms second-order in fluctuations, denoted by O(δ2), are negligible due to the small magnitude of the fluctuations, based on the assumption in Eq. (2.39). This leaves us with the linearised first-order 12 2.1. Neutron balance equations 2. Theory noise. From here, a Fourier transform from the time domain to the frequency domain is applied, F { X(r,t) } = ∫ ∞ −∞ X(r,t)e−iωtdt ≡ X(r,ω), (2.50) such that F { ∂ ∂t X(r,t) } = iωX(r,ω) (2.51) for a quantity X(r,t) with corresponding Fourier-pair function X(r,ω). This results in iω vg δϕg(r, ω) =∇ · [Dg,0(r)∇δϕg(r, ω)] + Ng∑ g′=1 Σs0,g′→g,0(r)δϕg′(r, ω) + Ng∑ g′=1 δΣs0,g′→g(r, ω)ϕg′,0(r) + χp g(r) keff [ 1 − β̃(r) ] Ng∑ g′=1 νg′(r)Σf,g′,0(r)δϕg′(r, ω) + χp g(r) keff [ 1 − β̃(r) ] Ng∑ g′=1 νg′(r)δΣf,g′(r, ω)ϕg′,0(r) + Nd∑ i=1 χd g,i(r)λi(r)δCi(r, ω) − Σt,g,0(r)δϕg(r, ω) − δΣt,g(r, ω)ϕg,0(r), (2.52) and iωδCi(r, ω) = β̃i(r) keff Ng∑ g′=1 νg′(r)Σf,g′,0(r)δϕg′(r, ω) + β̃i(r) keff Ng∑ g′=1 νg′(r)δΣf,g′(r, ω)ϕg′,0(r) − λi(r)δCi(r, ω). (2.53) The goal is to isolate the unknown quantity δϕg(r,ω). Notice that from Eq. (2.53) we can express the tidal fluctuations of the concentration of precursors of delayed neutrons as a function of ϕg,0(r) and δϕg(r,ω) as δCi(r, ω) = 1 iω + λi(r) [ β̃i(r) keff Ng∑ g′=1 νg′(r)Σf,g′,0(r)δϕg′(r, ω) + β̃i(r) keff Ng∑ g′=1 νg′(r)δΣf,g′(r, ω)ϕg′,0(r) ] = 1 keff β̃i(r) iω + λi(r) [ Ng∑ g′=1 νg′(r)Σf,g′,0(r)δϕg′(r, ω) + Ng∑ g′=1 νg′(r)δΣf,g′(r, ω)ϕg′,0(r) ] . (2.54) Now substitute δCi(r, ω) with Eq. (2.54) in Eq. (2.52), and put the δϕg-dependencies on the left-hand- side and the ϕg,0-dependencies on the right-hand side, yielding iω vg δϕg(r, ω) − ∇ · [Dg,0(r)∇δϕg(r, ω)] − Ng∑ g′=1 Σs0,g′→g,0(r)δϕg′(r, ω) − χp g(r) keff [ 1 − β̃(r) ] Ng∑ g′=1 νg′(r)Σf,g′,0(r)δϕg′(r, ω) − Nd∑ i=1 χd g,i(r)λi(r) { 1 keff β̃i(r) iω + λi(r) [ Ng∑ g′=1 νg′(r)Σf,g′,0(r)δϕg′(r, ω) ]} + Σt,g,0(r)δϕg(r, ω) = Ng∑ g′=1 δΣs0,g′→g(r, ω)ϕg′,0(r) + χp g(r) keff [ 1 − β̃(r) ] Ng∑ g′=1 νg′(r)δΣf,g′(r, ω)ϕg′,0(r) + Nd∑ i=1 χd g,i(r)λi(r) { 1 keff β̃i(r) iω + λi(r) [ Ng∑ g′=1 νg′(r)δΣf,g′(r,ω)ϕg′,0(r) ]} − δΣt,g(r, ω)ϕg,0(r). (2.55) 13 2. Theory 2.2. Point-kinetics After factorizing the four terms including the macroscopic fission cross-section we get iω vg δϕg(r, ω) − ∇ · [Dg,0(r)∇δϕg(r, ω)] − Ng∑ g′=1 Σs0,g′→g,0(r)δϕg′(r, ω) + Σt,g,0(r)δϕg(r,ω) − 1 keff { χp g(r) [ 1 − β̃(r) ] + Nd∑ i=1 χd g,i(r)λi(r)β̃i(r) iω + λi(r) } Ng∑ g′=1 νg′(r)Σf,g′,0(r)δϕg′(r, ω) = Ng∑ g′=1 δΣs0,g′→g(r, ω)ϕg′,0(r) − δΣt,g(r, ω)ϕg,0(r) + 1 keff { χp g(r) [ 1 − β̃(r) ] + Nd∑ i=1 χd g,i(r)λi(r)β̃i(r) iω + λi(r) } Ng∑ g′=1 νg′(r)δΣf,g′(r,ω)ϕg′,0(r). (2.56) Notice that the right-hand side of Eq. (2.56) is dependent on the static solution of the neutron flux ϕg,0(r) alone, while the fluctuations are isolated to the macroscopic cross-sections. Given that the group constants and cross-sections are known (see Section 3.1.2), the steady-state solution calculated (see Section 3.2.2) and the fluctuations of the macroscopic cross-sections are specified (see Section 3.2.3), the right-hand side of Eq. (2.56) is independent of the induced neutron noise. The dynamical equations of the linearised neutron noise in multi-group diffusion theory thus result in the following source-type problem iω vg δϕg(r, ω) − ∇ · [Dg,0(r)∇δϕg(r, ω)] − Ng∑ g′=1 Σs0,g′→g,0(r)δϕg′(r, ω) + Σt,g,0(r)δϕg(r,ω) − 1 keff { χp g(r) [ 1 − β̃(r) ] + Nd∑ i=1 χd g,i(r)λi(r)β̃i(r) iω + λi(r) } Ng∑ g′=1 νg′(r)Σf,g′,0(r)δϕg′(r, ω) = S(r,ω), (2.57) with the noise source S(r,ω) ≡ Ng∑ g′=1 δΣs0,g′→g(r, ω)ϕg′,0(r) − δΣt,g(r, ω)ϕg,0(r) + 1 keff { χp g(r) [ 1 − β̃(r) ] + Nd∑ i=1 χd g,i(r)λi(r)β̃i(r) iω + λi(r) } Ng∑ g′=1 νg′(r)δΣf,g′(r,ω)ϕg′,0(r), (2.58) To define the noise source S(r,ω), we need to specify the fluctuations of the cross-sections in Eq. (2.58), which induces the neutron noise. These are modelled as a point-like perturbations of the macroscopic cross-sections, referred to as a absorber of variable strength, δΣα,g(r, ω) = A(ω) δ(r − rp), (2.59) where δΣα,g(r, ω) is a general fluctuation of a macroscopic cross-section of type α, A(ω) is the perturbation strength and rp is the location of the perturbation. Motivated in Section 1.3, only frequency-independent perturbations of strength 0.1 % are applied to the macroscopic abortion cross-section in this work, and Eq. (2.59) reduces to δΣa,g(r) = 0.001 Σa,g(rp) δ(r − rp). (2.60) Given that the group constants, static flux and source term are known, the induced neutron noise can then be calculated by solving Eq. (2.57). This is described in Section 3.2.3. 2.2 Point-kinetics A common method of analysing the kinetic response of a reactor is to employ point-kinetics (PK). The PK-framework is based on separating the dependencies of the time dependent neutron flux, which results in a separation of the neutron noise into two different components. This allows a means of measuring the amount of information that is carried by the neutron noise about the perturbation. The derivation of the 14 2.2. Point-kinetics 2. Theory Zero-power transfer function (ZPTF), used to validate the developed solvers, is also derived within this framework. This framework is referred to as point-kinetics since one of the components will be spatially distributed according to the static flux, so it will be spatially independent of the perturbation. The perturbation will only affect the amplitude of this quantity as a time-dependent amplitude factor. This component is hence referred to as the point-kinetic component, since the reactor is reduced to a single spatial point with respect to the response of the perturbation [17, Ch. 9]. 2.2.1 Shape and point-kinetic components The point-kinetic framework is based on factorising the scalar neutron flux into a solely time-dependent amplitude factor P (t) and a shape function Ψg(r,t) such that ϕg(r,t) = P (t)Ψg(r,t). (2.61) Following the same procedure as before, the time-dependent quantities in Eq. (2.61) are separated into their time-average and tidal fluctuations, as in Eq. (2.38). Enforcing the initial condition of the shape function ϕg,0(r) = P0Ψg,0(r), (2.62) leads to the following expression of the linearised neutron noise: δϕg(r,t) ≈ δϕpk g (r,t) + P0δΨg(r,t), (2.63) where δϕpk g (r,t) = δP (t) P0 ϕg,0(r) (2.64) is referred to as the point-kinetic component of the neutron noise, and the term P0δΨg(r,t) is the so- called shape component or space-dependent term. Note how the spatial dependence of the latter term is expressed by the fluctuation of the shape function, while the spatial dependence of the point-kinetic component is solely described by the static flux. Introducing the amplitude factor and shape function requires some constraints through normalisation. In this work, the shape function is normalised according to [17, p. 469] ∂ ∂t ∫ G∑ g=1 1 vg ϕ+ g,0(r)Ψg(r,t) dr = 0, (2.65) and the time-averaged amplitude factor is set equal to the time average reactor power [17] P0 = ∫ G∑ g=1 κf,gΣf,g,0(r)ϕg,0(r) dr, (2.66) where κf,g is the average energy release per fission event for energy group g. In Eq. (2.65), ϕ+ g,0(r) is the adjoint function of the static flux. Analogously to Eq. (2.46), the adjoint static equation can be identified as an eigenvalue problem, and can be solved in the same way, only with transposed matrices. The details of the numerical solution are presented in Section 3. To quantify to which extent the total noise is comprised of the PK-component, a measure referred to as the PK-deviation measure, is introduced as [22] dev(ω) = 1 H ∫ H/2 −H/2 ∣∣∣∣∣∣∣∣ δϕ(z,ω) δϕpk(z,ω) − 1 ∣∣∣∣∣∣∣∣ dz, (2.67) here written for the axial dimension r = z, where H is the height of the active part of the reactor, δϕ(ω) is the renormalised total noise after being summed over all energy groups and δϕpk(ω) is the renormalised PK-component of the neutron noise after being summed over all energy groups. This integral measure is equal to zero if the total noise is fully point-kinetic, and grows as the shape component increases. In other words, it measures deviation from point-kinetics of the total noise, and thus the amount of information about the perturbation carried by the neutron noise. 15 2. Theory 2.2. Point-kinetics 2.2.2 Zero-power transfer function Following the procedure described in detail in [17], [19] and [22], coupling the diffusion equation of the neutron noise, Eq. (2.36), with the factorised neutron flux, Eq. (2.61), and enforcing Eq. (2.65) yields the point-kinetic equations. From these, the so called point-reactor model can be derived [17], formulated in frequency-domain as [22] δP (ω) P0 = G0(ω)δρ(ω), (2.68) with an analytical expression for the point-kinetic zero-power transfer function G0(ω) = 1 iω ( Λ0 + ∑Nd i=1 βi iω+λi ) (2.69) for Nd precursors of delayed neutrons, and Λ0 = 1 F0 ∫ G∑ g=1 1 vg ϕ+ g,0(r)Ψg,0(r) dr (2.70) F0 = ∫ G∑ g=1 (νΣf )g,0(r)ϕ+ g,0(r)Ψg,0(r) dr. (2.71) The amplitude of the ZPTF is expected to be constant at roughly 1 β̃ within the frequencies ω ∈ [λ, β̃/Λ0], where λ is the total decay constant [23] and Λ0 is the mean neutron generation time. Within this range, the phase of the ZPTF is expected to be close to zero. Because of these characteristics, this frequency range is referred to as plateau frequencies. In general, the analytical expression for the zero-power transfer function exists for all systems [17], so the point-kinetic component of the neutron noise can be expressed analytically through the zero-power transfer function by inserting Eq. (2.68) in (2.64), yielding δϕpk g (r,ω) = G0(ω)δρ(ω)ϕg,0(r). (2.72) Note that the point-kinetic component, formulated through the zero-power transfer function, is solely determined by the static solution, its adjoint, and the change in reactivity. Also note that the analytical expression of the ZPTF is independent of any neutron noise in the system. As described in [22], the ZPTF can also be computed through the dynamic solution. Utilising the system response when a source term is introduced allows one to express the fraction on the left-hand side of Eq. (2.68) in terms of the neutron noise, and consequently the ZPTF as well. This is done by combining Eqs. (2.63) and (2.64), performing some algebraic manipulations, utilizing the initial condition Eq. (2.62) and the normalisation Eq. (2.65) with the time-dependence separated according to Eq. (2.38). After performing a Fourier transform as in Eq. (2.50), one gets δP (ω) P0 = ∫ ∑G g=1 1 vg ϕ+ g,0(r)δϕg(r,ω) dr∫ ∑G g=1 1 vg ϕ+ g,0(r)ϕg,0(r) dr . (2.73) This method yields an expression of the fraction δP (ω) P0 , and now, according to Eq. (2.68), the ZPTF can be expressed by dividing Eq. (2.73) with the change in reactivity from the source term δρ(ω), yielding Gdyn 0 (ω) = 1 δρ(ω) ∫ ∑G g=1 1 vg ϕ+ g,0(r)δϕg(r,ω) dr∫ ∑G g=1 1 vg ϕ+ g,0(r)ϕg,0(r) dr ≈ 1 δρ ∫ ∑G g=1 1 vg ϕ+ g,0(r)δϕg(r,ω) dr∫ ∑G g=1 1 vg ϕ+ g,0(r)ϕg,0(r) dr , (2.74) where we have assumed that the reactivity change is independent of the frequency. Using Eq. (2.68), this assumption allows us to express the reactivity change as δρ(ω) ≈ δρ = δP (ω1) P0 Gdyn 0 (ω1) , ∀ω (2.75) where ω1 is an arbitrarily fixed frequency, as long as it is within the range of physically relevant frequencies. This version of the ZPTF is dependent on the neutron noise, and thus dependent on the dynamical solution and corresponding source term. This is in contrast with the ZPTF in Eq. (2.69), which is solely dependent on the static solution. By comparing these two functions with each other, one can verify that the solver works as intended. 16 3 Methodology A nuclear reactor consists of different geometric structures with varying length scales. Starting from the largest, the core. It consists of a large pressure vessel with the fissile material contained in the centre, often with some shielding surrounding it. The fissile material is then arranged in several hundred fuel assemblies. The fuel assemblies all have different levels of enrichment and burn-up. Interlaced are control rods with a high boron concentration. The fuel elements themselves also contain tens to hundreds of fuel rods. The fuel rod is in turn made of a central uranium pellet encased in a strong alloy, with a small gap between filled with helium. From this, one can see that the components of a nuclear reactor span several orders of magnitude. The pressure vessel measures several meters across [24], while a fuel assembly is roughly a decimetre or two in width [14]. A fuel pin narrows down to about a centimetre, and the helium gap within the fuel pin is even smaller, measuring less than a millimetre [14]. To accurately simulate every physical phenomena across the entire reactor with the resolution of the smallest component quickly becomes unfeasible, especially if one simulates every dynamic property of the reactor. For this reason, our system of interest is reduced to a one-dimensional system along the height of a fuel pin, where the properties of the system are based on a horizontal slice of a unit cell. From this, average system properties are used to calculate the steady-state, or static, neutron flux distribution along the fuel pin using deterministic methods. The deterministic method requires both spatial and energy discretisation of the scalar neutron flux. The static flux is then used to calculate the neutron noise’s spatial distribution throughout the system, which is calculated by introducing a perturbation in the frequency domain. For each reactor system, it is assumed that the material properties of the reactor are static. This means no burn-up of nuclear material, thermal material expansion or temperature changes. The properties are also spatially homogenised in Section 3.1.2 in the axial plane with respect to the spatial neutron flux distribution. We also assume that the properties do not change in the axial direction, effectively homogenising the entire reactor core. 3.1 Stochastic methods The first section of the methodology details the generation of some macro-parameters, referred to as group constants, using stochastic Monte-Carlo methods, that are capable of simulating the neutronics of an arbitrary geometry. For this purpose, unit cells, consisting of a fuel pin and surrounding coolant are defined, and used as input in the stochastic calculations. 3.1.1 Reactor Specifications The Monte Carlo simulations in the next section require physical dimensions and material properties of each of the four reactor types. The information was gathered from a variety of academic and industrial sources. Some specifications were calculated by the authors, such as material densities and temperatures. The temperatures and densities were assumed to be constant in each material. The filler gas in each fuel pin was assumed to be helium for all different reactor types [25, Sec. 9.4.5]. No information regarding the helium density of the LMFRs was found by the authors, instead, it was assumed to be the average of the LWRs’ helium density. The temperatures of the LWRs’ unit cells were calculatedusing the method described in [18, Ch. 5]. 17 3. Methodology 3.1. Stochastic methods Boiling Water Reactor This design is based on the design of a SVEA Optima 8x8 BWR fuel assembly [14, Table A-3]. The specifications are collected in Table 3.1 and the geometry is shown in Figure 3.1. The coolant has a significant steam void fraction during nominal operations, which we have assumed to be 40%. The fuel consists of uranium oxide enriched to 3.5% 235U. The helium pressure was assumed to be 7.5 bar [25, Sec. 9.4.5]. Table 3.1: Boiling water reactor fuel specifications. Based on SVEA Optima 8×8 BWR fuel assembly. Pin pitch 12.8 mm Pin active fuel length 3712 mm Pin assembly symmetry Square Fuel pellet diameter 10.4 mm Fuel material UO2 Fuel density 10.3 g cm−3 Fuel 235U enrichment 3.5% Fuel temperature 1293 K Gap width 0.075 mm Gap temperature 632 K Cladding material Zircaloy-2 Cladding thickness 0.605 mm Cladding temperature 604 K Coolant material H2O Coolant density 0.436 g cm−3 Coolant void fraction 40% Coolant temperature 559 K Figure 3.1: A horizontal representation of the BWR fuel pin. The inner most region represents the fuel followed by a helium gap and then the Zircaloy-2 cladding. Each fuel pin is surrounded by coolant. The square symmetry is also shown. Pressurised Water Reactor The representative PWR is based on a Westinghouse 17×17 fuel assembly [14, Table A-4]. The specifica- tions are collected in Table 3.2 and the geometry is shown in Figure 3.2. The water coolant was assumed to be free of vapour. The gap helium pressure was assumed to be 15 bar [25, Sec. 9.4.5]. Table 3.2: PWR specifications. Based on Westing- house’s 17×17 fuel design. Pin pitch 12.6 mm Pin active fuel length 3658 mm Pin assembly symmetry Square Fuel pellet diameter 8.19 mm Fuel material UO2 Fuel density 10.3 g cm−3 Fuel 235U enrichment 4.0% Fuel temperature 1184.9 K Gap width 0.085 mm Gap temperature 633 K Cladding thickness 0.57 mm Cladding material Zircaloy-4 Cladding temperature 608 K Coolant material H2O Coolant density 0.709 g cm−3 Coolant temperature 570 K Figure 3.2: A horizontal representation of the PWR fuel pin. The innermost region represents the fuel, followed by a helium gap and then the Zircaloy-4 cladding. The square symmetry is also shown. 18 3.1. Stochastic methods 3. Methodology Lead-Cooled Fast Reactor Our reference LFR design is based on a theoretical 80 MWt design which is cooled by pure elemental lead [13]. The design employs two different layers of cladding on the fuel pins, with the outer layer protecting against corrosion and the inner layer for structural purposes. The material temperatures were calculated from the average temperature distributions from [13, Fig. 9]. Figure 3.3: A horizontal representation of the LFR fuel pin. The innermost region represents the fuel, followed by a helium gap. The cladding consists of two different alloys, an outer corrosion-resistant and an inner structural layer. The fuel pin lattice has a hexagonal symmetry. Table 3.3: Lead-cooled fuel pin specifications. Based on a theoretical 80 MWt design. Pin pitch 13.0 mm Pin active fuel length 1010 mm Pin assembly symmetry Hexagonal Fuel pellet diameter 10.0 mm Fuel material UN Fuel density 13.78 g cm−3 Fuel 235U enrichment 12.0% Fuel temperature 1163 K Gap width 0.2 mm Gap temperature 973 K Inner cladding thickness 0.50 mm Inner cladding material 15-15Ti-AIM1 Inner cladding temperature 783 K Outer cladding thickness 0.20 mm Outer cladding material Fe-10Cr-4Al-RE Outer cladding temperature 765 K Coolant material Pb Coolant density 10.47 g cm−3 Coolant temperature 759 K Sodium-Cooled Fast Reactor Our reference SFR design is based on the Superphénix reactor [15, Tab. 1-2]. The reactor is fuel by MOX-fuel, which is a mix of uranium- and plutonium oxide, see Table 3.5. We have assumed that the uranium component is natural uranium. There exists a central helium void in each fuel pin. Figure 3.4: A horizontal representation of the SFR fuel pin. The innermost region represents the cen- tral void, followed by the fuel and a helium gap. The cladding consists of a 15-15Ti steel. The fuel pin has a hexagonal symmetry. Table 3.4: Sodium-cooled fuel pin specifications. Based on the Superphénix reactor. Pin pitch 9.8 mm Pin active fuel length 1000 mm Pin assembly symmetry Hexagonal Centre void diameter 2.0 mm Fuel pellet diameter 7.14 mm Fuel materiala UnatO2 + PuO2 Fuel density 10.54 g cm−3 Fuel PuO2 fraction 15.0% Fuel temperature 1500 K Gap width 0.115 mm Gap temperature 1132 K Cladding thickness 0.565 mm Cladding material 15-15Ti [26] Cladding temperature 783 K Coolant material Na Coolant density 10.47 g cm−3 Coolant temperature 745 K aPuO2 is replaced by 235UO2 for some specific calculations. 19 3. Methodology 3.1. Stochastic methods The plutonium in the SFR’s MOX was primarily sourced from reprocessed Gas Cooled Reactor (GCR) fuel [27, p. 96]. The isotopic composition of such plutonium is shown in Table 3.5 [28]. For calculations related to point-kinetic properties, the plutonium in the SFR was replaced with an equal mass of 235UO2, in order to compare the systems at different points in the fuel cycle. Table 3.5: Isotopic composition of plutonium in spent GCR nuclear fuel [28]. Isotope Mass fraction Pu238 0.00 wt% Pu239 76.70 wt% Pu240 20.10 wt% Pu241 2.80 wt% Pu242 0.40 wt% 3.1.2 Monte-carlo parameter estimation The geometries described in the previous section contain physical structures that are composed of ma- terials at different temperatures, elements and densities. These all affect the neutronics of the reactor in different ways. The goal of the following section is to preserve reaction rates for each reaction type when the dimensionality of the reactor is reduced to one. For this purpose, the neutron transport code Serpent-2 was utilised [20]. It uses a continuous-energy Monte-Carlo neutron transport method to simulate various spatially and energy dependent variables, such as cross-sections and neutron flux dis- tributions for arbitrary geometries. Our calculations used the JEF-3.1.1 Nuclear Data Library [29], which provides cross-section data for different reaction types at different incident neutron energies. Serpent-2 was used to calculate the neutron flux in a horizontal slice of each unit fuel cell, as shown by the geometries in Figures 3.1-3.4. For the LWRs, a thermal scattering correction was applied to the hydrogen in the water. The boundary conditions of each fuel pin are periodic, as to simulate an infinite system in the horizontal plane. The nature of the 2D-slice implies that the system effectively becomes infinite in the axial direction as well. A total of 1000 neutron generations, each consisting of 100 000 neutrons, was simulated, leading to a total of 100 000 000 neutron generations for each reactor system. The individual neutrons are binned to a discretised neutron flux into so-called micro-energy groups. The micro group is user-defined, with the default being {h}NH h=1, where NH = 70. See Appendix A.1 for details about the micro-group energy boundaries. The flux for each micro group is calculated as, ϕh = ∫ Eh−1 Eh ∫ S ϕ(r,E) dEd2r, (3.1) with the associated cross-sections for each group Σh = ∫ Eh−1 Eh ∫ S Σ(r, E)ϕ(r, E) dEd2r∫ Eh−1 Eh ∫ S ϕ(r, E) dEd2r . (3.2) which is used to calculate the spatially dependent parameters for each macro-energy group Σg = ∑ h∈g Σhϕh∑ h∈g ϕh . (3.3) This preserves the reaction rates when transforming from a system that is continuous in both energy and space, to one with a scalar flux with discretised energy. The selection of the macro-energy groups is a non-trivial task, which is described in detail by Moghrabi [30] for the curious reader. We have selected to use a standard 7 group structure (CASMO-code) for the LWRs [31] and a different 8 group structure for the LMFRs [32, Tab. 4.2]. The energy boundaries of the two groups are shown in Tables 3.6-3.7. Note the shift towards lower neutron energies for the thermal group as compared with the fast group. All parameters that were calculated through Serpent-2 and that are utilised in the next sections are shown in Table 3.8. These parameters are all spatially homogenised through the previously mentioned method. Due to the sheer number of numerical values associated with the parameters in Table 3.8, the 20 3.1. Stochastic methods 3. Methodology Table 3.6: Macro-group energy boundaries of the thermal systems. The lower bound of energy group g = 7 is 0 eV. g Upper bound [eV] 1 2.00 · 107 2 8.21 · 105 3 5.53 · 103 4 4.00 · 100 5 6.25 · 10−1 6 1.40 · 10−1 7 5.80 · 10−2 Table 3.7: Macro-group energy boundaries of the thermal systems. The lower bound of energy group g = 8 is 0 eV. g Upper bound [eV] 1 2.00 · 107 2 2.23 · 106 3 8.21 · 105 4 3.02 · 105 5 1.11 · 105 6 4.09 · 104 7 1.50 · 104 8 7.49 · 102 numerical values are not shown, except for a small selection in Table 3.9 that are used for validation purposes. The effective multiplication factor k∞ is larger than unity for each reactor, which can be expected since fresh fuel is assumed, without burnable absorbers or inserted control rods, and for an infinite system. The effective delayed neutron fraction β̃ varies slightly between the reactors, due to material and neutron spectrum differences, as expected. Table 3.8: Spatially homogenised parameters calculated using Serpent-2. Subscript g indicates that the parameter is dependent on the macro-group and subscript i is the delayed neutron family number. Name Symbol Fission cross-section Σf,g Isotropic scattering cross-section Σs0,g′→g Anisotropic (P1) scattering cross-section Σs1,g′→g Absorption cross-section Σa,g Prompt fission neutron spectrum χp g Average number of neutrons released per fission event νg Neutron speed vg Delayed fission neutron spectrum χd g Delayed neutron decay constant λi Delayed neutron fraction βi Table 3.9: Multiplication factor k∞ and effective delayed neutron fraction β̃ for our four reactors as simulated by Serpent-2 based on technical data from Section 3.1.1. BWR PWR LFR SFR SFR (235U) k∞ 1.22 1.36 1.08 1.32 1.26 β̃ 0.709% 0.691% 0.746% 0.352% 0.712% 21 3. Methodology 3.2. Deterministic methods 3.2 Deterministic methods The latter part of the methodology utilises the parameters that are calculated in the previous sections. These methods are entirely deterministic and thus quite computationally performant. Two different solvers, one for the static flux and one for the noise, were developed. These solvers require further discretisation, specifically in space, where the axial fuel pin is divided into vertical nodes. 3.2.1 Spatial discretisation So far, the scalar neutron flux has been described as ϕ(r, E) or with discretised energy as ϕg(r). Here, the flux is discretised in the axial direction of the fuel pin and denoted as ϕn g , where n ∈ [1, NZ ] is the position index and g be the macro-energy group. In this context, the flux can be interpreted as volume-averaged quantities, where the fuel pin has been divided into vertical slices of equal height. The number of spatial points NZ is chosen such that the distance between nodes is smaller than the accepted minimum accurate distance of 1 − 2 cm for LWRs, when using a discretisation based on finite differences. We have chosen the smaller spatial resolution of 1 mm, as to reduce the errors from the finite differences significantly. The discretised form of the balance equations requires minimal modifications from the continuous version. The interaction rate preservation from Section 3.1.2 implies that a direct substitution can be made for all terms except the diffusion term from Eq. (2.8). ∂ ∂zϕg(zn + ∆z/2) represents the central finite-difference at position zn + ∆z/2, which lies between our discretised nodes n and n+ 1. It is defined as, ∂ϕg ∂z (zn + ∆z 2 ) ≈ ϕn+1 g − ϕn g ∆z , (3.4) where ∆z is the distance between nodes. Using such derivatives, the second-order derivative can be calculated at each node n as ∂2ϕg ∂z2 (zn) ≈ ∂ϕg ∂z (zn + ∆z 2 ) − ∂ϕg ∂z (zn − ∆z 2 ) ∆z ≈ ϕn−1 g − 2ϕn g + ϕn+1 g ∆z2 . (3.5) The discretised multi-group transport equation can thus be written as −Dg ∂2ϕg ∂z2 (zn) + Σt,gϕ n g − Ng∑ g′=1 Σs,g′→gϕ n g′ = 1 keff (1 − β̃)χp g Ng∑ g′=1 Σf,g′ϕn g′ + 1 keff β̃χd g Ng∑ g′=1 Σf,g′ϕn g′ . (3.6) Marshak boundary conditions are assumed by stipulating a gradient in Eq. (3.5), based on the flux at the boundaries. At the bottom of the fuel pin, the flux is ϕg(0), which lies ∆z/2 below our lowest discretised node ϕ1 g. The Marshak conditions relate the neutron current Jg(0) with the neutron flux according to, Jg(0) = −1 2ϕg(0). (3.7) The neutron current is defined as, Jg(0) = −Dg ∂ϕg ∂z (0) ≈ Dg ϕg(0) − ϕ1 g ∆z/2 . (3.8) Setting Eq. (3.7) equal to Eq. (3.8) yields, ϕg(0) = ϕ1 g 1 + ∆z 4Dg , (3.9) which gives the neutron current at z = 0 using Eq. (3.8) as, Jg(0) = 1 2 · ϕ1 g 1 + ∆z 4Dg = 2ϕ1 g 4 + ∆z Dg . (3.10) 22 3.2. Deterministic methods 3. Methodology This is directly proportional to the gradient, according to, ∂ϕg ∂z (0) = 2ϕ1 g 4Dg + ∆z , (3.11) which allows the second-order derivative to be defined at z1 as, while remembering that 0 = z1 − ∆z/2, ∂2ϕg ∂z2 (z1) ≈ ∂ϕg ∂z (z1 + ∆z 2 ) − ∂ϕg ∂z ϕg(0) ∆z ≈ − [ 2 ∆z(4Dg + ∆z) + 1 ∆z2 ] ϕ1 g + 1 ∆z2ϕ 2 g. (3.12) The relation in Eq. (3.12) supersedes the one defined in Eq. (3.5) at z = z1. A similar expression exists for the upper border, but with ϕg(H) and Jg(H) instead, where H is the fuel pin’s height. At the top boundary, the Marshak condition given is as: Jg(H) = 1 2ϕg(H), (3.13) which can be compared with Eq. (3.7). Note that Eqs. (3.5), (3.6), and (3.12) form a linear system of equations that can be represented using a vector and two matrices. Let the neutron flux vector be ϕ = [ϕ1,ϕ2, ...,ϕNG ], where ϕg = [ϕ1 g, ϕ 2 g,..., ϕ NZ g ], and M and F be matrices. This allows us to rewrite Eq. (3.6) as an eigenvalue problem, Mϕ = 1 keff Fϕ =⇒ F −1Mϕ = 1 keff ϕ. (3.14) 3.2.2 Steady-state solver After performing the discretisations described in Section 3.1.2 and 3.2.1, and representing the terms in Eqs. (2.46) and (2.47) as matrices, the problem can be formulated as finding the eigenvector ϕ that corresponds to Eq. (3.14) such that 1/keff is the largest eigenvalue. The largest eigenvalue can be found using an iterative power iteration method [33], or by using the accelerated Wielandt shift method [34]. The procedure for the power iteration for solving Eq. (3.14) requires an initial guess ϕ0 and is performed iteratively for each iteration i, ϕi+1 = 1 ki M−1F ϕi, (3.15) ki+1 eff = ki eff ϕi · ϕi+1 ϕi · ϕi , (3.16) ϕi+1 = ϕi+1 max g,z ϕi+1 . (3.17) The initial flux ϕ0 was set to unity and k0 eff was also set to unity. This can be further accelerated using Wielandt’s method, which reduces the magnitude of the dominance ratio d = |λ1/λ0|, where λ0 and λ1 are the largest and second largest eigenvalues, respectively. The convergence rate for the power iteration is largely determined by the starting guess and the dominance ratio [35]. By subtracting 1/kestϕ i from both sides of Eq. (3.14), one can rewrite the equation using a new matrix W and Λ, as: Mϕ − 1 kest F ϕ = 1 keff Fϕ − 1 kest F ϕ, (3.18) W ϕ = ΛFϕ. (3.19) This equation is also solvable using power iteration. Selecting the estimated kest is done according to kn est = kn − 0.001, which ensures that Λ remains positive. For computational efficiency, W was factorised using an LUPQ-decomposition, such that, P W Q = LU . (3.20) where P and Q are permutation matrices, and L and U are lower and upper triangular matrices, respectively. The triangular matrices allow for efficient solution via forward and backwards substitution. Thus, the procedure becomes, ϕi+1 = ΛiQU−1L−1P F ϕi, (3.21) 23 3. Methodology 3.2. Deterministic methods where U−1 and L−1 denote solutions to triangular systems, not explicit matrix inverses. The adjoint flux vector ϕ+ is calculated using the same methods as the static flux, where the matrices M and F are transposed. Determining whether the iterative method has converged can be performed by calculating the residual r, ri ≡ ∥Mϕi − 1 ki F ϕi∥, (3.22) the results of which are shown in Figure 3.5. 0 5 10 15 i [1] 10!12 10!9 10!6 10!3 r [1 ] BWR PWR SFR LFR (a) Static flux ϕi. 0 5 10 15 i [1] 10!12 10!9 10!6 10!3 r [1 ] BWR PWR SFR LFR (b) Static adjoint flux ϕ+,i. Figure 3.5: Residual r of the steady-state solver as a function of the iteration number i for the two different static fluxes. 3.2.3 Dynamic noise solver In a similar manner as previously, the source-type problem, Eq. (2.57), is discretised as S(ϕ0, gp, np) = iω vg δϕn g +Dg −δϕn−1 g + 2δϕn g − δϕn+1 g ∆z2 + Σt,gδϕ n g − Ng∑ g′=1 Σs,g′→gδϕ n+1 g − 1 keff [ χp g(1 − β̃) + Nd∑ i=1 χd g iω + λi ] Ng∑ g′=1 νg′δϕn g , (3.23) where S(ϕ0, gp, np) is the perturbation source term. Note that the boundary conditions from Eq. (3.12) still apply. ϕ0 is the static flux as calculated in the previous section, gp and np are indices for where the disturbance is located in space and energy group and are defined as: S(ϕ0, gp, np) ≡ { 0.001 · Σa,g · ϕn g , if n = np and g = gp 0, otherwise. (3.24) This relation is rewritten using matrices, Mδϕ = S, (3.25) where S is the source vector from Eq. (3.24) and Mδϕ is the right hand side of Eq. (3.23). This is an ordinary linear system of equations and is solved using standard numerical methods. The noise flux δϕ is calculated for each ω of interest. In our case, 100 logarithmically spaced values between 1 × 10−4 and 1 × 104 Hz. 24 4 Results and Discussion In this section, the solver is first validated by comparing the steady-state ZPTF with the noise-dependent ZPTF. The PK-deviation measure introduced in Section 2.2 is then presented for the different systems, along with a motivation for the choice of measurement. Once this overview of the complete results for all systems are presented, some representative cases are analysed in more detail. Following this, dev(ω) is computed for systems of the same size. This is done in order to distinguish spectra-related effects from effects related to the reactor size. Finally, a brief discussion about the spatial relaxation lengths concludes the results. 4.1 Model verification To verify that the numerical solver works as intended, the ZPTF function calculated by the two methods outlined in Section 2.2 are compared. This results in the ”analytical solution” given by Eq. (2.69), which is solely dependent of the steady-state neutron flux and denoted Gsta 0 (ω), and the ”numerical ZPTF” given by Eq. (2.74), which is dependent on the neutron noise and denoted Gdyn 0 (ω). The two quantities are normalised to the same magnitude for frequency ω1 = 14.5 Hz, as described by Eq. (2.75), which is well within the plateau region as described by the bounds in Table 4.1. The magnitudes of these two functions are compared in Figure 4.1a, and the phases of the ZPTFs are compared in Figure 4.1b. The most apparent result is that both the phase and the magnitude of the Gdyn 0 (ω) show the characteristic shape of the ZPTF, with the plateau region within the bounds λ and β/Λ0, where the magnitude is constant at 1/β and the phase is close to zero. The values of the characteristic amplitudes and frequencies are generated by the stochastic simulation in Serpent- 2 (except Λ0, which is computed as Eq. (2.70)) and benchmarked against the shape of ZPTF. The computed values of the plateau amplitudes and frequencies for all systems are presented in Table 4.1. They align as expected, thus validating the solver. There is, however, a small discrepancy between the lower bound of the plateau region and the upper bound. This is most clearly seen in Figure 4.1b, where the ZPTFs do not appear to be centred around the plateau region bounds. Table 4.1: Characteristic amplitudes and frequencies of the ZPTF, where λ corresponds to the upper plateau frequency, β̃/Λ0 to the lower plateau frequency, and 1/β̃ is the plateau amplitude. λ [Hz] β̃/Λ0 [Hz] 1/β̃ [1] BWR 0.0816 94.0 141 PWR 0.0795 93.7 145 SFR 0.0881 1860 284 LFR 0.0846 5010 134 A comparison between Gdyn 0 (ω) and Gsta 0 (ω) provides a validation that the solver seems to work as in- tended. It is clear that the dynamic ZPTFs align very well with their static counterparts, especially for the fast systems. The thermal systems show a slightly worse alignment, which is indicative that some- thing in the solver might not work perfectly. However, this only affects the solver for frequencies larger than approximately 50 Hz, and even though the exact impact on the result is unknown, the qualitative similarities between Gdyn 0 (ω) and Gsta 0 (ω) means that the solver likely captures the general behaviour of the system, even for larger frequencies. Furthermore, in real applications of noise diagnostics, the 25 4. Results and Discussion 4.1. Model verification classified noise sources that are physically relevant are below 50 Hz, examples of which are: boiling in BWR at 0.1-10 Hz [2, p. 1700], control rod oscillation 5-6 Hz [3], fuel rod oscillations 0.6-10 Hz [36] and PWR core barrel oscillations 8-15 Hz [37]. Given these two points, and the purpose of this study, the accuracy of the solvers are deemed sufficient. 10!4 10!2 100 102 104 ! [Hz] 100 101 102 103 104 105 jjG 0 (! )jj [1 ] jjGsta0 (!)jj jjGdyn0 (!)jj 1=~-h 6; ~-=$0 i ZPTFs and benchmarks BWR PWR SFR LFR (a) Magnitude of the ZPTFs as a function of frequency for all systems. -80 -60 -40 -20 0 P h as e [d eg ] BWR PWR 10!3 100 103 ! [Hz] -80 -60 -40 -20 0 P h as e [d eg ] SFR 10!3 100 103 ! [Hz] LFR arg ! Gsta 0 (!) " arg ! Gdyn 0 (!) " [6;-=$0] ZPTFs and Benchmarks (b) Phase of the ZPTFs as a function of frequency for all systems. Figure 4.1: Comparison between the frequency-dependence of the analytical expression of the ZPTFs, denoted Gsta 0 (ω) (represented by full lines), and the ZPTF through system response (represented by diamond markers). The dashed lines mark the 1/β̃ magnitude, and the dot-dashed lines mark the bounds of the plateau region, [λ, β̃/Λ0]. 26 4.1. Model verification 4. Results and Discussion To evaluate whether the static flux is accurate, the flux calculated in the infinite system (see Section 3.1.2) is compared with the flux calculated by our static solver. A deviation is defined as |ϕg − ϕS2 g |/ϕS2 g , where ϕS2 g is the flux as calculated by Serpent-2 in the infinite system, and ϕg is the static flux as calculated in Section 3.2.2. Both fluxes are normalised such that ∑ g ϕg = 1. The magnitude of the deviation falls within the range 15−0.02 %, where the LMRFs have a larger deviation on average. When the axial length of the reactor is increased, all average deviations decrease, thus indicating that the flux differences are caused by physical factors such as neutron leakage. This inspires confidence that our static solver is accurate. 10!1 101 103 105 107 Neutron energy [eV] 10!1 100 R el at iv e n eu tr on . u x d ev ia ti o n [% ] BWR PWR SFR LFR Figure 4.2: The relative deviation for each energy group, and for each system. The deviation is defined as |ϕg − ϕS2 g |/ϕS2 g , where ϕS2 g is the flux as calculated by Serpent-2; see Section 3.1.2. 27 4. Results and Discussion 4.2. PK-deviation measure comparison 4.2 PK-deviation measure comparison In order to get a comprehensive overview of how different source terms affect the spatial shape of the neutron noise, a measure is introduced that encapsulates the deviation from point-kinetics given a source term, thus reflecting the feasibility of performing noise-diagnostic techniques on the system at hand. The specific choice of the measure is somewhat arbitrary, as long as it properly represents the system at hand while accurately reflecting the difference between systems. In this report, the measure is defined as Eq. (2.67). This integral measure has been used in previous works [22], and seems to correctly represent the system responses for all source terms implemented in this work. 5000 10000 15000 20000 25000 0.2 0.4 0.6 0.8 1 N o rm al iz ed n oi se [1 ] PWR 22000 23000 24000 25000 0.1 0.2 0.3 0.4 PWR, g = 7 2000 4000 6000 8000 Vector index [1] 0.2 0.4 0.6 0.8 1 N o rm al iz ed n oi se [1 ] LFR 7200 7400 7600 7800 8000 Vector index [1] 0.05 0.1 0.15 LFR, g = 8 Total noise Point-kinetic Shape Components (a) Normalised total noise vector (in blue), with its PK-component (in red) and shape-component (in dashed yellow) for the PWR and LFR systems. The left-hand plots show all energy groups, while the right-hand plots only show the lowest energy group. 10!4 10!2 100 102 104 ! [Hz] 10!2 10!1 100 101 102 P K -d ev ia ti on [1 ] PWR: dev(!) LFR: dev(!) PWR: PNg g devg(!) LFR: PNg g devg(!) (b) Two different integral measures showing the deviation from point-kinetics for the PWR and LFR. Figure 4.3: Comparison between two different PK-deviation measures: the well-behaved measure dev(ω) used in this work, and the faulty groupwise integrated measure ∑Ng g devg(ω). The source term consists of a 0.1 % increase of the absorption cross-section at the lowest energy-group, localised 45 % off-centre. The noise and its components at frequency 2.78 Hz are presented in (a), and the PK-deviation measures are presented in (b). To strengthen the claim that this measure is sufficient, Figure 4.3 shows the robustness of the chosen measure compared to an alternative measure. Compared to the PK-deviation measure presented in Eq. (2.67), the integration of this alternative measure is performed groupwise before averaging over groups, and is therefore denoted ∑Ng g devg(ω). The source term consists of a perturbation of 0.1 % increase of the macroscopic absorption cross-section of the lowest energy group, localised 45 % off-centre. The alternative measure ∑Ng g devg(ω) incorrectly reflects a larger deviation from point-kinetics for the fast systems than what can be seen in Figure 4.3a (only PWR and LFR presented here as thermal and fast systems respectively). This is due to the extremely low static flux in the lowest energy group in the fast 28 4.2. PK-deviation measure comparison 4. Results and Discussion systems (seen in Section 4.3, Figure 4.5a). From Eq. (2.64), it follows that the corresponding point- kinetic component is extremely small. The total noise in this energy group consists almost fully of its shape component, and the corresponding integrand becomes disproportionately large in this group since δϕg(ω) δϕpk g (ω) >> 1. This does not happen if the noise is summed over groups before integrating, as in Eq. (2.67). The alternative measure exaggerates the deviation from point-kinetics to the point where, for this specific source term, the LFR appears to deviate more than the PWR, as seen in Figure 4.3b. This is clearly a faulty representation when compared to the shape of the noise, Figure 4.3a. On the other hand, the chosen integral measure dev(ω) does not exhibit this faulty behaviour, and correctly reflects the systems at hand. Only this measure is retained in the following analysis. Now that the integral measure of the deviation from point-kinetics has been chosen, it can be computed for a plethora of cases to get an overview of how the different systems respond to different source terms. The number of possible source terms in Eq. (2.58) are plenty: the perturbations can be applied to all cross-sections in the noise source, or even any combination of them; they can be applied to different energy groups, and any combination of them; the strength of the perturbation can be frequency-dependent; the location of the perturbation, zp, can be changed. Considering all of them would exceed the scope of this work, and not all of these noise sources are physically relevant. As described in Section 1.3, the noise source has been limited to time-independent perturbations of the macroscopic absorption cross-section with strength 0.1 %. No combination of energy groups will be perturbed simultaneously. The spatial location of the perturbation will be either centred at zp = 0, or 35 % off-centre at zp = −0.35 ·H. The deviation from point-kinetics of the neutron noise induced by these noise sources is presented in Figure 4.4. The first general observation is that the LWRs deviate far more from point-kinetics than the LMFRs, with the PWR deviating most, closely followed by the BWR. The LFR deviates far less than these two, and the SFR deviates the least. This is in part due to the longer diffusion lengths of fast neutrons, which are more abundant in the LMFRs. This lack of thermal neutrons leads to more homogeneous shape of the noise, and therefore a smaller deviation from point-kinetics. Due to the lower coolant density in BWR, as compared with PWR, the neutron spectrum will also be slightly harder than that of the PWR, from the reduced moderation. Another reason for the difference in point-kinetics is the physical size differences between the systems. The higher density of the LMFRs’ coolant, as compared with that of LWRs, constrains the size of the reactor, due to safety concerns regarding construction and maintenance [13]. A smaller reactor can generally be assumed to behave more like a point-reactor. The size difference will be discussed further in Section 4.4. The shaded bounds in Figure 4.4 constitute the variation of dev(ω) within the system depending on the energy group of the perturbation. At plateau frequencies, this difference in deviation does not exceed approximately 0.1 for any system, i.e.∣∣∣∣max g {dev(ω)}g − min g {dev(ω)}g ∣∣∣∣ ≲ 0.1. (4.1) Since these differences within systems are small absolute, while the differences between systems can differ by orders of magnitude, the bound for systems appears negligible in systems with greater PK-deviation. Typically, the deviation from point-kinetics increases when the energy group of the source term is lower. So perturbations in the absorption cross-section for lower energy groups typically tend to cause a greater deviation from point-kinetics than if the change in cross-section occurs in a higher energy group. This result is consistent for almost all systems and all perturbations seen in the figure. The only exceptions to this general behaviour are the off-centred source terms of the LWRs, whose PK-deviation increases for higher energy group perturbations. The final general observation from Figure 4.4 is that when the perturbation is axially located closer to the periphery of the core, the PK-deviation increases. An interesting result is that this effect is more drastic for systems with larger PK-deviation. This also explains the apparent decrease in the difference in dev(ω) with respect to the energy group of the source term. This is attributable to the increased PK- deviation, which downplays the absolute energy group-related differences of around 0.05 when plotted on a logarithmic scale. The main takeaway from these results is that, contrary to the LMFRs, the LWRs deviate strongly from point-kinetics already at plateau frequencies. The smaller deviation from point- kinetics exhibited by the LMFRs at plateau frequencies will complicate the task of performing noise diagnostics drastically, especially for the SFR system. However, as discussed above, if the perturbation is axially located further towards the periphery of the core, the PK-deviation increases. This is especially 29 4. Results and Discussion 4.3. Analysis of representative cases interesting for the case regarding the LFR, which might deviate sufficiently from point-kinetics to allow for some noise diagnostic applications. But due to the complexities associated with noise diagnostics, the neutron detector configuration, and the difference between the modelled LFR and other LFR-type systems, the specifics of which noise term provides the best outlook for noise diagnostics will not be specified further. 10!4 10!2 100 102 104 ! [Hz] 10!2 10!1 100 101 102 d ev (! ) [1 ] Centered perturbation, zp = 0 10!4 10!2 100 102 104 ! [Hz] 10!2 10!1 100 101 102 O,-centered perturbation, zp = !0:35H BWR PWR SFR LFR Figure 4.4: Deviation from point-kinetics as a function of frequency, with each system being represented as: BWR - yellow full line; PWR - blue dashed line; SFR - red dot-dashed line; LFR - green dotted line. The shaded regions represent the energy-group dependence of the source term for each system. Spatially, the left subfigure has a centred source term axially, with zp = 0, and the right subfigure has an off-centred source term, with zp = −0.35H. 4.3 Analysis of representative cases To better understand the general results presented in Figure 4.4, it is valuable to delve deeper into some representative cases. In this subsection, a deeper analysis is performed for all reactor systems, with both a centred (zp = 0) and an off-centred (zp = −0.35H) perturbation of the same energy group, specifically energy group g = 7. This energy group is chosen for two reasons. Firstly, it will provide the greatest PK- deviation according to Figure 4.4 (with the exception of g = 8 for the LMFRs). Secondly, the magnitude of the static flux in this energy group is comparable to the other energy groups (in contrast to g = 8 for the LMFRs). Since the neutron noise is computed from the steady-state solution of the flux, a natural starting point is to consider the static flux. Both the groupwise static flux ϕg,0(z) and the total static flux ϕ0(z) are presented in Figure 4.5a and 4.5b respectively. The latter quantity is normalised and summed over all energy groups according to ϕ0(z) = ∑Ng g ϕg,0(z) max z (∑Ng g ϕg,0(z) ) . (4.2) The static flux is expected to be cosine-shaped due to the homogeneity of the systems, which is confirmed in the figures and strengthens the validation of the steady-state solver. The groupwise static flux for both LWRs shows a similar distribution of static flux throughout their energy groups. The largest differences are in their thermal energy groups, for g > 3. The BWR exhibits a small but steady decrease of the static flux as energy decreases, whereas the thermal flux of the PWR seems constant for all of these groups. For the LMFRs, the general character is also similar for both systems. They differ slightly in energy groups g = 1,5,6,7,8, where the LFR has a smaller flux. 30 4.3. Analysis of representative cases 4. Results and Discussion 5000 10000 15000 20000 25000 0 0.5 1 ? g ;0 (z ) [1 ] BWR 5000 10000 15000 20000 25000 0 0.5 1 ? g ;0 (z ) [1 ] PWR 1000 2000 3000 4000 5000 6000 7000 8000 0 0.5 1 ? g ;0 (z ) [1 ] SFR 1000 2000 3000 4000 5000 6000 7000 8000 Vector index [1] 0 0.5 1 ? g ;0 (z ) [1 ] LFR (a) Groupwise static flux, where the vector index consists of Ng · NZ discrete point running from the highest energy group g = 1 and one border z = −H/2 to the lowest energy group g = Ng and the other border z = H/2. -185 0 185 0 0.5 1 ? 0 (z ) [1 ] BWR -182 0 182 0 0.5 1 ? 0 (z ) [1 ] PWR -50 0 50 0 0.5 1 ?