Validation of the Westinghouse BWR nodal core simulator POLCA8 against Serpent2 reference results

(1)

INOM

EXAMENSARBETE TEKNISK FYSIK, AVANCERAD NIVÅ, 30 HP

STOCKHOLM SVERIGE 2021,

Validation of the Westinghouse BWR nodal core simulator

POLCA8 against Serpent2 reference results

MATHILDE GAILLARD

KTH

SKOLAN FÖR TEKNIKVETENSKAP

(2)

nodal core simulator POLCA8 against Serpent2 reference results

Mathilde Gaillard

Date: February, 2021

Degree Project in Engineering Physics (30 ECTS credits)

Double engineering degree between Grenoble INP Phelma and KTH Royal Institute of Technology

Supervisor at KTH: Vasily Arzhanov

Examiner at KTH: Prof. Jan Dufek

(3)

Mathilde Gaillard, MSc Master Thesis Report

TRITA-SCI-GRU 2021:018

Royal Institute of Technology School of Engineering Sciences KTH SCI

SE-100 44 Stockholm, Sweden URL: www.kth.se/sci

Westinghouse Eletric Sweden AB, Västerås, Sweden, February 2021 ii

(4)

Aknowledgement

First of all, I would like to thank Mr Jan Dufek, one of my teacher at KTH, for helping me to find this master thesis thanks to his contacts, as well as I would like to thank Mr Jean-Marie Le Corre, the first person I met from Westinghouse, for helping me to find an interesting topic in the field I am most interested in within the company, and Mr Magnus Limbäck, for helping me with the technical details once I was accepted for this master thesis.

I would also like to thank and express my gratitude to Mr. Petri Forslund Guimarães for accepting me into the project and taking me on as a student. He spent a lot of time teaching me all the important aspects of the subject, and we had many interesting conversations about how to solve the problems we were facing. He was also very helpful with his comments on how to write the thesis correctly, and I know that all the lessons he taught me will be useful to me throughout my career.

I would like to thank Mr Vasily Arzhanov for advice on how to write this master thesis.

I would also like to thank my friends and coworkers from KTH with whom we helped each other during this particular period of teleworking and kept ourselves motivated during the whole thesis period. I am thinking in particular of one of my closest friends, Dr Benjamin Portal, who supported me both morally and technically with his advice on graphics and writing.

Finally, I would like to thank all the people who have helped me, or who have helped to make the project a success, in one way or another.

A special thought belongs to my father, who, I know, from where he is, would have been really proud of me.

(5)

Abstract

When a new nodal core simulator is developed, like all other simulators, it must go through an extensive verification and validation effort where, in the first stage, it will be tested against appropriate reference tools in various theoretical benchmark problems. The series of tests consist of comparing several geometries, from the simplest to the most complex, by simulating them with the nodal core simulator developed and with some higher order solver representing the reference solution, in this case on the Serpent2 Monte Carlo transport code. The aim of this master’s thesis is to carry out one part of these tests. It consisted in simulating a three-dimensional (3D) 2x2 mini boiling water reactor (BWR) core with the latest version of the Westinghouse BWR nodal core simulator POLCA8, and in comparing the outcome of these simulations against Serpent2 reference results.

Prior to this work, POLCA8 was successfully tested on a 3D single-channel benchmark problem using the same Serpent2/POLCA8 methodology. However, this benchmark problem considered in this work is challenging in several aspects. Indeed, the nodal core simulator should accurately predict the eigenvalues and power distributions against reference results, and this by taking into account axial leakage, resulting from the passage from two-dimensional (2D) infinite lattice physics calculations to 3D simulations, or strong axial flux gradients due to the insertion or withdrawal of the control rods after a certain depletion. This last effect is known as the Control Blade History (CBH) effect and will be the main focus of this study. In addition to the development of a new version of the nodal core simulator, a new version of the Westinghouse deterministic transport code PHOENIX5 is also under development. The accuracy of PHOENIX5 was indirectly tested through this benchmark by providing the cross sections for the POLCA8 simulations. In addition, Serpent2 based nodal cross sections were generated to POLCA8 to provide means of comparing these two sets of nodal cross section data. The results obtained lead to the conclusion that the CBH model gives very good results, especially with regard to all power distributions, and especially those after the removal of the control bars when needed most.

keywords: Nodal Core Analysis, Monte Carlo Methods, CBH Effects

Westinghouse Eletric Sweden AB, Västerås, Sweden, February 2021 iii

(6)

Sammanfattning

När en ny nodal-kärnsimulator utvecklas, som alla andra simulatorer, måste den genomgå en omfattande verifierings och valideringsinsats där den i det första steget kommer att testas mot lämpliga referensverktyg i olika teoretiska riktmärkesproblem. Testserien består av att jämföra flera geometrier, från den enklaste till den mest komplexa, genom att simulera dem med den utvecklade nodkärnsimulatorn och med någon högre ord- ningslösning som representerar referenslösningen, i detta fall på Serpent2 Monte Carlo-transportkoden. Syftet med detta examensarbete är att genomföra en del av dessa tester. Den bestod av att simulera en tredimensionell (3D) 2x2 mini-kokande vattenreaktor (BWR) -kärna med den senaste versionen av Westinghouse BWR- nodalkärnasimulator POLCA8, och att jämföra resultatet av dessa simuleringar mot Serpent2-referensresultat.

Före detta arbete testades POLCA8 framgångsrikt på ett 3D-enkanaligt riktmärkesproblem med samma Serpent2 / POLCA8-metodik. Detta riktmärkesproblem som beaktas i detta arbete är dock utmanande i flera aspekter. I själva verket bör nodkärnsimulatorn noggrant förutsäga egenvärdena och kraftfördelningarna mot referensresultat, och detta genom att ta hänsyn till axiellt läckage, resulterande från övergången från tvådimensionella (2D) oändliga gitterfysikberäkningar till 3D-simuleringar eller starkt axiellt flöde gradienter på grund av att styrstavarna sätts in eller dras ut efter en viss utarmning. Denna sista effekt är känd som CBH-effekten (Control Blade History) och kommer att vara huvudfokus för denna studie. Förutom utvecklingen av en ny version av nodal core-simulatorn är också en ny version av Westinghouse deterministiska transportkod PHOENIX5 under utveckling. PHOENIX5: s noggrannhet testades indirekt genom detta riktmärke genom att tillhandahålla tvärsnitt för POLCA8-simuleringar. Dessutom genererades Serpent2-baserade nodtvärsnitt till POLCA8 för att tillhandahålla medel för att jämföra dessa två uppsättningar av nodtvärsnittsdata. De erhållna resultaten leder till slutsatsen att CBH-modellen ger mycket bra resultat, särskilt med avseende på alla effektfördelningar, och särskilt de som har tagits bort när man behöver mest.

Nyckelord: Nodal kärnanalys, Monte Carlo-metoder, CBH-effekter

(7)

Abreviations

2D: two-dimensional 3D: three-dimensional BA: Burnable Absorbers BWR: Boiling Water Reactor CBH: Control Blade History

CCCP: Current Coupling Collision Probability CC: Current Coupling

CP: Collision Probability CR Control Rods

CRAM: Chebyshev Rational Approximation Method LWR: Light Water Reactor

MC: Monte Carlo

PHX5 xs data: PHOENIX5 cross section data Serp xs data: Serpent2 cross section data

Westinghouse Eletric Sweden AB, Västerås, Sweden, February 2021 v

(8)

Aknowledgement ii

Abstract iii

Sammanfattning iv

Abreviations v

List of Figures x

List of Tables xi

Presentation of the company xii

1 Introduction 1

1.1 Nodal core analysis methodology . . . 1

1.2 Mission for the project . . . 3

2 Codes in use and cross-section generation 4 2.1 The stochastic Monte Carlo transport code Serpent2 . . . 4

2.2 The deterministic transport code PHOENIX5 . . . 5

2.3 The nodal core simulator POLCA8 . . . 6

2.3.1 Cross-section model . . . 6

2.3.2 Pin power reconstruction methodology . . . 8

2.4 Cross-section generation methodology . . . 9

3 Specification of benchmark problems 11 3.1 Description of the configurations . . . 11

3.1.1 Benchmark configuration . . . 11

3.1.2 The Control Blade History (CBH) effect . . . 13

3.2 Setting up the reference solution in Serpent2 . . . 14

3.2.1 Bias in the multiplication factor . . . 14

3.2.2 Convergence of the fission source . . . 15

3.2.3 Void Profile . . . 17

3.3 Quantities of interest and statistics . . . 17

3.4 Normalization of the results . . . 17

(9)

4 Results 19

4.1 Analysis at the fresh core conditions . . . 19

4.1.1 Core multiplication factor predictions . . . 19

4.1.2 Nodal power distributions . . . 20

4.1.3 Pin power distributions . . . 21

4.2 Analysis of the depleted core conditions . . . 23

4.2.1 Core multiplication factor predictions . . . 24

4.2.2 Nodal power distributions . . . 26

4.2.3 Pin power distributions . . . 30

5 Conclusions 38 Bibliography 40 Appendix 41 A Development of results analysis tools 42 A.1 Resuts format . . . 42

A.2 The Matlab analysis tool . . . 43

B Gantt chart and diagram of tasks performed 45 C Results at fresh core conditions 47 C.1 Axial pin power distribution . . . 47

C.2 Radial pin power distribution . . . 48

C.2.1 Axial node 5, PHOENIX5 cross section data . . . 49

C.2.2 Axial node 5, Serpent2 cross section data . . . 50

D Results with PHOENIX5 cross section data, depleted core conditions 57 D.1 Axial pin power distribution . . . 57

D.1.1 Control rods corner . . . 57

D.1.2 Detector corner . . . 58

D.2 Radial pin power distribution . . . 59

D.2.1 Axial node 13, control rods half inserted . . . 60

D.2.2 Axial node 5, CBH model activated . . . 61

D.2.3 Axial node 5, CBH model inactivated . . . 62

Westinghouse Eletric Sweden AB, Västerås, Sweden, February 2021 vii

(10)

E Results with Serpent2 cross section data, depleted core conditions 67

E.1 Axial pin power distribution . . . 67

E.1.1 Control rods corner . . . 67

E.1.2 Detector corner . . . 68

E.2 Radial pin power distribution . . . 69

E.2.1 Axial node 13, control rods half-inserted . . . 70

E.2.2 Axial node 5, CBH model activated . . . 71

E.2.3 Axial node 5, CBH model inactivated . . . 72

(11)

List of Figures

1.1 Lattice cell used in this project . . . 2

3.1 Fuel assemblies at the bottom of the core . . . 12

3.2 Axial subdivision of the core . . . 12

3.3 The imposed axial coolant density profile in POLCA8 and Serpent2 calculations . . . 13

3.4 Stabilisation of the Shannon entropy as function of the number of inactive cycles . . . 15

3.5 Convergence of the eigenvalue as function of the number of inactive cycles . . . 16

4.1 Comparison between PHX5/POLCA8 and Serp2/POLCA8 against 3D Serpent2 reference simulation for the nodal power in the bundle A/01, CBH model inactivated . . . 21

4.2 Differences in the axial pin power distribution between POLCA8 simulations and the reference simulation inside bundle A/01, pin next to CR corner . . . 22

4.3 Differences in the axial pin power distribution between POLCA8 simulations and the reference simulation inside bundle A/01, pin next to detector corner . . . 22

4.4 Eigenvalue evolution and differences between PHOENIX5/POLCA8 and Serpent2 reference simulations . . . 24

4.5 Eigenvalue evolution and differences between Serpent2/POLCA8 and Serpent2 reference simulations . . . 25

4.6 Differences in the axial power distribution between PHOENIX5/POLCA8 and the Serpent2 reference simulations inside bundle A/01 during burnup with CBH model activated . . . 27

4.7 Differences in the axial power distribution between PHOENIX5/POLCA8 and the Serpent2 reference simulations inside bundle A/01 during burnup with CBH model inactivated . . . 28

4.8 Differences in the axial power distribution between Serpent2/POLCA8 and the Serpent2 reference simulations inside bundle A/01 during burnup with CBH model activated . . . 28

4.9 Differences in the axial power distribution between Serpent2/POLCA8 and the Serpent2 reference simulations inside bundle A/01 during burnup with CBH model inactivated . . . 29

4.10 PHOENIX5/POLCA8 axial power distributions at the withdrawal of the control rod . . . 29

4.11 Serpent2/POLCA8 axial power distributions at the withdrawal of the control rod . . . 30

4.12 Differences in the axial pin power distributions between POLCA8 and Serpent2 reference simulations inside bundle A/01, the pin next to the control rod corner pin with CR inserted, PHX5 xs data . . . 31

4.13 Differences in the axial pin power distributions between POLCA8 and Serpent2 reference simulations inside bundle A/01, the pin next to the control rod corner pin with CR withdrawn, PHX5 xs data . . . 31

ix

(12)

4.14 Differences in the axial pin power distributions between POLCA8 and Serpent2 reference simulations inside bundle A/01, the pin next to the control rod corner pin with CR inserted, Serp2 xs data . . . 32 4.15 Differences in the axial pin power distributions between POLCA8 and Serpent2 reference

simulations inside bundle A/01, the pin next to the control rod corner pin with CR withdrawn, Serp2 xs data . . . 32 4.16 Differences in the axial pin power distributions between POLCA8 and Serpent2 reference

simulations inside bundle A/01, the pin next to the detector corner pin with CR inserted, PHX5 xs data . . . 33 4.17 Differences in the axial pin power distributions between POLCA8 and Serpent2 reference

simulations inside bundle A/01, the pin next to the detector corner pin with CR withdrawn, PHX5 xs data . . . 33 4.18 Differences in the axial pin power distributions between POLCA8 and Serpent2 reference

simulations inside bundle A/01, the pin next to the detector corner pin with CR inserted, Serp2 xs data . . . 34 4.19 Differences in the axial pin power distributions between POLCA8 and Serpent2 reference

simulations inside bundle A/01, the pin next to the detector corner pin with CR withdrawn, Serp2 xs data . . . 34 4.20 Differences in the radial pin power distribution between PHOENIX5/POLCA8 and Serpent2

reference simulations inside axial node 13 of bundle A/01 at the burnup point of CR withdrawal (CRout) and with the CBH model activated . . . 35 4.21 Differences in the radial pin power distribution between PHOENIX5/POLCA8 and Serpent2

reference simulations inside axial node 13 of bundle A/01 at the burnup point of CR withdrawal (CRout) and with the CBH model inactivated . . . 36 4.22 Differences in the radial pin power distribution between Serpent2/POLCA8 and Serpent2

reference simulations inside axial node 13 of bundle A/01 at the burnup point of CR withdrawal (CRout) and with the CBH model activated . . . 37 4.23 Differences in the radial pin power distribution between Serpent2/POLCA8 and Serpent2

reference simulations inside axial node 13 of bundle A/01 at the burnup point of CR withdrawal (CRout) and with the CBH model activated . . . 37

(13)

List of Tables

3.1 Bias in the eigenvalue for the evaluated 3D mini-core . . . 15

4.1 Core eigenvalues . . . 20

4.2 Reactivity components . . . 20

4.3 Summary of radial pin power distributions . . . 23

4.4 Core eigenvalues at the withdrawal of the control rod . . . 25

4.5 Reactivity components at the withdrawal of the control rod for PHOENIX5 cross section data . 26 4.6 Reactivity components at the withdrawal of the control rod for Serpent2 cross section data . . 26

4.7 Differences in the radial power distributions between POLCA8 simulations and the Serpent2 reference solution at the burnup point of the withdrawal of the control rod . . . 26

4.8 Summary of radial pin power errors at the control rod withdrawal point (PHX5 xs data) . . . . 35

4.9 Summary of radial pin power distributions at the withdrawal point. Serp2 xs data . . . 36

xi

(14)

Presentation of the company

Westinghouse Electric Compagny LLC is an American nuclear power company created in 1999. It was the former nuclear power division of the original Westinghouse Electric Corporation. Westinghouse Electric Corporation was founded in 1886 and helped in developing electric infrastructure throughout the United States thanks to the first industrial AC system for generating power. The company helped the US government’s military program for nuclear energy applications by building the reactor on the world’s first nuclear submarine, and was also instrumental in the development and commercialization of nuclear energy systems for electric power generation. The world’s first pressurised water reactor (PWR) was designed and built by Westinghouse in 1957 in Shippingport, Pennsylvania, U.S. Nowadays, Westinghouse nuclear technologies are used in several countries such as France, EDF is a long-time licensee of them, or Sweden with the Ringhals Nuclear Power Plant.

In Sweden, Westinghouse designed and built 12 nuclear power reactors between 1966 and 1985. As part of the Swedish nuclear program, a nuclear fuel manufacturing plant was established in 1966 in Västerås.

Westinghouse bought it in 2000. This is one of the most modern nuclear fuel manufacturing plants in the world with a capacity of 900 tonnes of uranium per year and a production of nuclear fuel, fuel components and control rods. Within this manufacturing plant, an engineering service centre was developed to participate in the research, development, production, testing and licensing of nuclear fuel. This engineering service centre also acts for the maintenance of nuclear reactors such as boiling water reactor (BWR), PWR and water water energetic reactor (VVER). Part of their work is to develop a new version of reactor simulator such as nodal core simulator.

(15)

Chapter 1

Introduction

1.1 Nodal core analysis methodology

Steady-state neutronics analysis of nuclear reactors relies on accurate computational methods for obtaining the steady-state distribution of the free neutron population everywhere within the core. In Light Water Reactor (LWR) analysis, computational methods for small-scale calculations are usually based on the neutron transport theory, i.e. solving the Boltzmann transport equation. On the other hand, large-scale calculations are usually based on a simplified version of the transport equation, namely, the diffusion equation [1]. Although the transport equation is the most fundamental and accurate description of the spatial, energy and angular distribution of neutrons, performing this type of calculation on an entire reactor core requires extensive computational resources rendering such an approach rather impractical. Moreover, to analyse the behavior of a reactor core, a lot of different simulations have to be considered, such as performing core follow calculations, evaluating thermal margins during core surveillance or performing nuclear design calculations and fuel loading pattern optimization. Consequently, a two-step methodology is adapted in nodal core analysis to predict the neutronic behavior in the reactor core which comprises of two-dimensional (2D), multi-group lattice transport calculations generating nodal cross section data to the nodal core simulator for various domains of the core, and subsequent three-dimensional (3D) few-group nodal calculations considering the full reactor core.

The first step in this two-step nodal core analysis methodology is to perform calculations with a 2D neutron transport code based on the more accurate neutron transport theory. These lattice calculations are not performed on the whole core, but on a sub-region of the core called lattice cell with a fine spatial and energy mesh.

Figure 1.1 represents the lattice cell used for this project. A lattice cell typically represents a 1 cm thick slice that contains a single fuel assembly plus half of the surrounding coolant gap in the radial direction. It is accurately modeled in 2D geometry with materials characterized by fine-group cross-sections. Using the reflective boundary condition, this unique lattice cell becomes an infinite large core, with only one type of fuel assembly [2]. Thanks to generalized equivalence theory, preserving reaction rates and net leakages in an average sense, the resulting neutron flux can be used for spatial homogenisation and condensation of cross-sections (with respect to energy), as well as for computing discontinuity factors, pin power factors and other physics data [3]. These homogenised nodal data will then be used in the second stage of the methodology by parameterizing them as functions of important local state parameters and subsequently extracted by table interpolations or polynomial fitting.

In the second step, a coarse spatial mesh is applied on the entire reactor core domain using diffusion theory and pretabulated homogenized data obtained by lattice physics (previous step). The reactor core is

1

(16)

Figure 1.1: Lattice cell used in this project

subdivided into axial nodes of around 20 cm high with a base corresponding to the lattice cell i.e., one such node corresponds to a single coarse-mesh point. In the nodal diffusion approach [4], the nodal balance equation, derived from the 3D steady-state few group neutron diffusion equation, is solved within each node. Rather than the conventional approach of finite differences to discrete the spatial variable, nodal methods use a high order or analytical expansion of the intra-nodal flux shape to achieve a higher degree of accuracy for a given node size. As a set of multi-dimensional equations is obtained, a transverse integration is often employed to get a set of coupled one-dimensional equations. The resulting system is then solved on a 3D core composed by homogeneous nodes characterized by the data previously generated.

To ensure that operating design and safety parameters, such as the minimum critical power ratio, stay below the specified limits, the power generated in individual fuel rods is then estimated. If one wants to estimate the fluxes in localized regions within the nodes, the nodal solution does not contain enough detail. For that, a more detailed flux distribution is approximated from the nodal solution by superimposing the fine spatial mesh transport solution upon the homogeneous intra-nodal flux solution, a procedure denoted flux reconstruction.

Using a two-step methodology instead of performing full core transport calculations directly on the problem domain will introduce errors in the final results because of some additional assumptions inherent in this procedure, such as:

• Application of reflective boundary conditions in the 2D lattice calculations thereby neglecting ant leakage and spectrum interactions between fuel assemblies at the stage of generating cross section data.

• Assuming some representative core conditions during depletion neglecting the detailed local history that is only known during the full core calculation.

These assumptions are addressed in the nodal core simulator by dedicated but approximative cross section and pin power form function correction terms trying to compensate for these effects at the core level. Consequently, assessing the capability of the nodal core simulator to handle and model these challenging real core conditions using the two-step methodology is of crucial importance.

(17)

1.2 Mission for the project

This M.SC. project focuses on the second part of the nodal core analysis and its predictive capability. The goal is to assess the accuracy of the PHOENIX5 [5] (the deterministic transport code of Westinghouse) and POLCA8 [5] (the nodal core simulator of Westinghouse) code package in simulating challenging benchmark problems anticipated from real core operation, this by comparing POLCA8 nodal and pin power distributions against corresponding 3D full core Serpent2 [6] Monte Carlo transport results. For this purpose, a 3D boiling water reactor (BWR) mini-core control blade history (CBH) benchmark has been constructed and evaluated in this work to assess the accuracy of the dedicated cross section and pin power cross sections implemented in POLCA8 to account for leakage effects and CBH.

In this evaluation, two sets of nodal cross section data (first step of the nodal core analysis) are used and compared:

• Cross section data generated by the deterministic transport code PHOENIX5.

• Cross section data generated by the stochastic Monte Carlo code Serpent2.

By applying nodal cross section data generated by Serpent2 in the POLCA8 calculations, an unambiguous comparison between POLCA8 and Serpent2 is obtained as the same transport solution methodology and cross section library is used both for cross section generation and for obtaining the 3D reference solution.

This thesis is organized as follows. In Chapter 2, a description of the three codes used to carry out this project is given, as well as an explanation of the methodology for generating multi-group cross-sections. In Chapter 3, the specification of benchmark problems will be given through the description of the configuration, how the reference simulation was set-up, and a presentation of the quantities of interest and statistics used will be given too. Before the conclusion, the Chapter 4 will present the numerical results and their analysis.

Westinghouse Eletric Sweden AB, Västerås, Sweden, February 2021 3

(18)

Codes in use and cross-section generation

The purpose of this section is to provide a basic description of the methodology of the three codes used for this project, which are the Serpent2 Monte Carlo stochastic transport code, the PHOENIX5 deterministic transport code and the POLCA8 nodal core simulator. A description of how the cross-sections are generated during lattice physics calculations will also be given, to facilitate the understanding of the work performed and to highlight the challenges faced during this project.

2.1 The stochastic Monte Carlo transport code Serpent2

The Monte Carlo (MC) method is a way of solving complex problems such as neutron transport by a stochastic approach using random numbers. This method is usually used on transport problems of high geometrical complexity which are difficult to solve by a deterministic approach.

One of the major advantages of the MC method is its simplicity: the transport equation does not have to be formulated to find the neutron flux inside the reactor as the a large batch of neutron histories is simulated and the results are combined. It also allows the method to be able to simulate all types of reactor geometries (either the simplest ones or the most complex ones). Another advantage of the MC method is its accuracy in complex problem as the statistical errors decay as the square root of the batch size: the larger the sample, the more accurate the results will be.

However, as previously mentioned, the MC method uses a stochastic approach, which will always induce some statistical errors in the results. Therefore, even with a large batch size, the final result will still come with some uncertainties. The other main disadvantage of the MC method is the computational resources needed to resolve complex problems. The more complex the problem, the greater the resources required will be, especially with a large neutron batch size, and especially for burnup simulations that requires many flux solutions.

The MC code used for this project is Serpent2. Serpent2 was first developed as a MC lattice transport code, and is now developed as multipurpose MC reactor physics burnup calculation code capable to simulate various complex systems, to generate group constants to fuel cycle analysis, or to calculate coupled multi-physics problems [6].

Serpent2, as other MC codes, uses a three-dimensional constructive solid geometry model. This model is built from elementary quadratic and derived surface types. Subsequently, these different surface types are used

(19)

to form two or three dimensional cells. It is also possible to create different levels inside the geometry by using universes which will describe those levels.

Serpent2 uses a combination of two particle tracking methods:

• the ray-tracing based surface tracking method

• the rejection sampling based delta-tracking method

This last method allows neutron tracks to be continued over several material boundaries without calculating the distances to the boundary surface, which leads to a speed-up in the transport simulation.

All neutron interaction data used in Serpent2 are coming from continuous-energy ACE format cross-section libraries. Then, all cross-sections can be reconstructed by using a single unionized energy grid. With this approach, it is possible to pre-calculate material-wise macroscopic cross-sections, and only one energy grid is used to interpolate microscopic cross-sections between two tabulated values. Like that, the global calculation time is significantly reduced. For the burnup simulations, all the data needed for isotopic transmutation (ra- dioactive decay, energy-dependent fission yields and isomeric branching ratios) are read from ENDF format data libraries. Serpent2 solves the Bateman depletion equations by using the CRAM (Chebyshev Rational Approximation Method) matrix exponential method [7]. This method can handle the entire nuclide system without any approximation for short-lived isotopes and cyclic processes, without step length limitation or numerical accuracy problems. Although more advanced time integration schemes are available in Serpent2, such as the Stochastic Implicit Euler method [8] [9] [10], a standard predictor-corrector scheme was utilized in this work.

To further speedup the calculations, Serpent2 employs a hybrid OpenMP/MPI parallelization scheme. In addition, in order to fit large problems into the computer memory, domain decomposition is also available to subdivide the problem into smaller pieces to be distributed over the used computing nodes.

Besides constituting the reference solution for the evaluated 3D mini-core problem, Serpent2 generates all the nodal cross section data necessary for nodal diffusion calculations utilizing its built-in lattice physics branch capabilities, such as the homogenized few-group macroscopic and microscopic cross-sections, the assembly discontinuity factors and pin-wise power distributions for pin-power reconstruction. The cross section data are all written in a Matlab m-file format that are further processed by Westinghouse in-house parser script

"serp2Latt" to convert these data to a format readable by the POLCA8 code package.

2.2 The deterministic transport code PHOENIX5

PHOENIX5 is an advanced lattice physics code developped by Westinghouse, based on the well-established HELIOS lattice code with the main objective to create a reliable and robust tool for 2D neutron and gamma transport calculations. Its main purpose is to accurately model light water reactor (LWR) fuel design, particularly to cope with all the complexities of modern BWR reactor, and to generate nodal cross-section data for a 3D simulator, POLCA8 as example. It can also be used in criticality analyses, control rod design, isotopic inventory or make 2D reference calculations for others codes. PHOENIX5 mainly uses the ENDF/B-VII.1 cross-section library with complements of JEFF-3, JENDL-4 and BROND-2.2 for gamma transport calculations [5].

Similar to HELIOS, PHOENIX5 uses the Current Coupling Collision Probability (CCCP) method to solve the neutron and gamma transport equations. In this method, the system is divided into space elements, which are

(20)

further divided into flat-flux regions. They are internally treated by collision probabilities (CP). The coupling of these space elements is done by the current coupling (CC) method using interface currents with angular discretization on the hemisphere of the incoming directions [11] [12]. In orderto perform depletion calculations and transmutation of isotopes, PHOENIX5 uses a standard predictor-corrector integration scheme employing linearized nuclide chains.

As already explained in the Chapter 1, deterministic transport codes are used in the first step of nodal core analysis to perform lattice physics calculations. The purpose of these calculations is to generate nodal cross section data which will then be used by a nodal core simulator.

2.3 The nodal core simulator POLCA8

In nodal core analysis, after feeding generated nodal cross section data to the nodal core simulator of use, 3D reactor core simulations are conducted. The objective of this type of code is to simulate various conditions foreseen for the design of a nuclear reactor core as well as perform core foloww calculations for real plant operation. For example, the nodal core simulator of Westinghouse, POLCA (Power On-Line CAlculation), was created in 1968 and has been used for many years to perform various BWR analyses such as core and fuel design calculations, as well as core-tracking and in-core instrumentation evaluations.

As already explained in the Chapter 1, a nodal core simulator will solve the nodal diffusion equation on a system subdivided into homogenized nodes, with each such node associated with the previously generated nodal cross section data. Dedicated physics models are employed in the nodal core simulator to handle different phenomena related to reactor physics, such as the cross section and pin power reconstruction being amongst the most important. Consequently, both these models are discussed in more depth in the following.

2.3.1 Cross-section model

The cross-section model of POLCA8 is based on a combination of two types of terms:

• The ’base’ macroscopic cross-section, Σ^base, obtained when a lattice code depletes the fuel at reference conditions, i.e. depletion calculations where only the fuel exposure (burnup) is allowed to change while all other state parameters are fixed to their so-called base values (nominal design values or typical average values). In this regard, the coolant density and its history constitute important exceptions to the above base condition and are included in the set of dependencies base cross sections posses.

• Deviation terms composed of instantaneous effects, depletion history effects and spatial homogenization effects on cross sections.

The cross section model is mathematically expressed as:

(21)

Σ = Σ^base(E, ρh, ρ) + δSG∆Σ^SG(E, ρh, ρ) + δCR∆Σ^CR(E, ρh, ρ, β ) + δDT∆Σ^DT(E, ρh, ρ)

+ c_B(E, ρ_h, ρ)h

C_B−C_B^basei + c_{X e}(E, ρ_h, ρ)h

NX e− N_{X e}^basei + d_Dop(E, ρ_h, ρ)h

pTf−q T_f^base

i + c_T_m(E, ρ_h, ρ)h

T_m− T_m^basei +

∑

i

σi(E, ρh, ρ, wh, δCR, ...)h

N_i− N_i^base(E, ρh, w_CBH)i

+ w_CBH(1 − δ_CR) ∆Σ^CBH,out(E, ρh, ρ) + δCR∆Σ^CBH,in(E, ρh, ρ) + ∆Σ^spat+ ∆Σ^leak+ ∆Σ^het,byp

(2.1) During reconstruction, cross-sections are evaluated at a detailed axial sub-node level with a dependencies on various state parameters, such as the burnup (expressed here as E, in MWd/kg), the coolant density history (ρhin kg/m³) or the instantaneous coolant density (ρ in kg/m³). To convert the sub-node cross-section into the corresponding nodal mesh used by the flux solver, axial homogenization is performed [13].

Some terms, such as the corrections for control rods and spacer grids (∆Σ^CRand ∆Σ^SG), the off-base fuel Doppler temperature (d_Dop.∆Tf), the xenon concentration (c_{X e}.∆NX e) and the soluble boron concentration (c_B.∆CB) are generated by PHOENIX5 lattice physics calculations (first step of the nodal core analysis methodology). However, as PHOENIX5 (or Serpent2) lattice physics calculations are performed for an infinite fuel lattice, certain assumptions are made during these lattice calculations, such as the absence of leakage. Therefore, correction terms are estimated and computed by POLCA8 itself to account for neutron leakage effects (∆Σ^leak), as well as for intra-nodal depletion effects (∆Σ^spat) [13].

During a burnup simulation, the control rods also experience a depletion of their own active absorber material. This will have an impact on the cross-sections, as well and POLCA8 takes this into account in the cross-section model. The depletion state of the control blade is quantified with the control blade depletion fraction, β , include in the set of dependencies of the correction term for the control rods. Although cross-sections are dependent on the coolant density history, some residual history effects are included by simulating real life operation. Therefore, an isotropic correction term based on explicit microscopic depletion, ∑iσi∆Ni, have been add to the model [13].

The very challenging CBH effect, the main target of analysis of this study, is also handled with dedicated correction terms. During the lattice physics calculations, the two extreme conditions of the CBH effect, when the control rods are fully inserted during fuel depletion, ∆Σ^CBH,in, and when they are fully withdrawn during fuel depletion, ∆Σ^CBH,out, are computed. POLCA8 then performs an interpolation between these two extreme states to obtain the cross-section correction required for a certain control blade history condition [13]. Note that by deactivating the CBH model in POLCA8, one is fully relying on that intra-nodal cross section model corrections will capture the CBH effect, an assumption that has proven to be rather poor one.

Finally, to cope with certain transient conditions in the core, a term with an explicit dependence on the moderator temperature is present, c_T_m.∆T_m, as well as a correction term to account for bypass water density and

(22)

soluble boron heterogeneities ∆Σ^het,byp. If detectors are present inside the core, their geometrical effects on the cross-sections are also taken into account thanks to one special correction term, ∆Σ^DT [13].

2.3.2 Pin power reconstruction methodology

In POLCA8, pin powers are calculated inside a materially homogeneous sub-node s, and are expressed as volumetric power densities. A sub-node is obtained after the breaking down of a node, if this node presents material heterogeneities in the axial direction due to the presence of burnable absorbers (BA)/enrichment zoning, spacer grids or control rod for examples. The calculation of pin powers is based on a combination of lattice code physics data and nodal core simulator results. The detailed relative power density is mathematically expressed as:

p(x, y, z) = c^norm

2 g=1

∑

S^rad_g (x, y) S^ax(z) P_g^hom(x, y, z) (2.2) The summation is over two energy groups and the basic assumption that the radial and axial heterogeneous dependencies, represented by the shape function, are separable is made. The radial fine structure shape function, S^rad_g (x, y), is carried over the lattice code evaluation and is given for the given sub-node (one value per pin cell and group). This transport solution accounts for the heterogeneous pin nature of the assembly and the CBH effect. The axial fine structure shape function, S^ax(z), accounts for the node the axial heterogeneities like spacer grids, control rod tips, and enrichment/BA discontinuities. It is energy group independent, and the calculations of it take place in the core itself. The coefficient c^normis the normalisation factor, determined such that p (x, y, z) averages P^sub, the average relative power density for the given sub-node. The quantity P^subis known from the axial homogeneization calculation. [14].

The last term of Equation 2.2, P_g^hom(x, y, z), is known as the homogeneous relative power density distribution per energy group inside a full node and is obtained by solving the two-group diffusion equation with realistic boundary conditions and with slowly varying (close to uniform) cross section. It accounts for global, smooth power tilts caused by an uneven leakage of neutrons from neighbouring nodes and by the fact that the assembly is depleted spatially differently in real life compared to the lattice calculations feeding the nodal core simulator with cross section data. It can be mathematically expressed by the following equation:

P_g^hom(x, y, z) = κ Σ^node_f,g Φ^node_g (x, y, z) Q_relπ_V

Vn

V_n^f

(2.3)

where κΣ^node_f,g is the intra-nodal cross section profile, Φ^node_g (x, y, z) is the smooth flux solution (i.e., plane wave expansion in each energy group), Q_relis the relative core thermal power (W ), πvis the volumetric core power density (W.m⁻³), V_nis the full node volume (m³) and V_n^f is the node so-called "fissile" volume, i.e. only multiplying axial segments of the node (m³) [14].

In Serpent2, one can access both nodal power distributions and pin power distributions inside the considered system through the use of different options such as the cpd card or by a mesh detector. All these power distributions are expressed in watt [6]. However, although data of same resolution can be accessed in POLCA8, powers are specified in terms of relative values with regard to volumetric power densities. Consequently, and before making comparisons between simulations, it is necessary to express results using the same unit. Since the results of Serpent2 are expressed in watt, the choice was made to convert all powers computed by POLCA8 to the unit of watt.

(23)

In terms of nodal powers as predicted by POLCA8, one way to obtain the absolute nodal power in watt for node n is given by Equation 2.4.







P_n^abs[W ] = P_n^rel× Q_rel× πv×V_n 1

N

N n=1

∑

P_n^rel= 1.0







(2.4)

where P_n^relis the relative nodal power density from POLCA8.

The method for converting pin powers is the same as that for nodal powers, with the exception of one term.

It is given by Equation 2.5 for pin i in node n.











p^abs_i,n [W ] = p^rel_i,n× Q_rel× πv×V_n× 1 N_f 1

N_f

∑

i=1

p^rel_i,n= P_n^rel











(2.5)

where p^rel_i,nis the raw data from POLCA8, and N_f is the number of fuel pins inside the node n.

It is important to mention that Equation 2.4 and Equation 2.5 and their normalisation are only valid if all nodes within the core have the same volume [14].

2.4 Cross-section generation methodology

As earlier described, the cross sections used during nodal simulations must be generated before starting the calculations. They are generated by performing 2D lattice physics calculations employing a transport code, i.e., in this project by using the deterministic transport code PHOENIX5 and the stochastic MC transport code Serpent2.

For these lattice calculations, reflective boundary conditions are applied because it is assumed that the fuel assembly properties depend mainly on the heterogeneous nature of the assembly itself and less on its location in the reactor core. As several fuel assembly types are used within the core geometry, this set of calculations must be applied for each of them. Indeed, each set of 2D lattice physics calculations represents a axial fuel assembly segment type defined by its lattice and nuclear design (i.e. geometry, enrichment, BA loading, etc.) [15].

For each such fuel assembly segment type, a series of independent lattice depletion calculations have to be carried out for given depletion histories, i.e., in this case for given void histories (or active coolant density histories), combined with a large number of momentaneous branch calculations representing perturbations to the conditions prevailing during depletion. Only the burnup is allowed to change during these lattice depletion calculations, with all the other state parameters fixed at their base-values (i.e. nominal design values or typical average values). In contrast, the branch calculations are taking into account momentaneous changes in the operating conditions, such as the coolant density, the control rod insertion status, the xenon concentration and the soluble boron concentration. These branches can be applied at different burnup steps chosen by the user.

They are normally performed independently for each state parameter, but a combination of different parameters can also be done. For example, the momentaneous active coolant density is always varied simultaneously with each one or a combination of the other state parameters. Furthermore, the depletion history cases should cover the range of the expected values of the relevant history parameter, so that any extrapolation can be avoided when reconstructing the cross section for each node in the core [15].

(24)

The base conditions defined for each coolant density history are described by the following state parameter setup:

• Hot full power internal and external bypass moderator density corresponding to a saturated temperature at nominal core pressure with no boiling in these bypasses.

• Hot full power active coolant density corresponding to a selected void condition (including zero void).

This also defines the coolant density history for the selected void content.

• Hot full power nominal power density.

• Hot full power nominal fuel temperature.

• No control rods present.

• No spacer grids present.

• No detector presents.

• Reference soluble boron concentration of 0 ppm.

• Equilibrium xenon at nominal power density.

The Control Blade History (CBH) having such a strong impact on the reactivity and the assembly-internal power distribution, requires a speacial treatment. Consequently, dedicated CBH lattice calculations are performed to generate CBH corrections to the base values of cross sections, discontinuity factors and pin power form functions. These CBH correction tables are then tabulated as function of the momentaneous coolant density, coolant density history and fuel assembly exposure [15].

The CBH lattice calculations are performed as follows:

• Control rods always inserted (i.e. rodded) base depletion calculation for each coolant density history at hot full power both voided and non-voided conditions.

• Rodded coolant density branch calculations from rodded base depletion cases for each momentaneous coolant density condition.

• Control rods withdrawal branch calculations from rodded base depletion cases also for each momentaneous coolant density condition.

Finally, it is worthwhile to mention that these lattice calculations can be very time demanding when considering such calculations based on Serpent2. Consequently, only a subset of all possible and supported branch cases were performed in this work, including:

• Coolant density branches.

• Zero xenon branches.

• Control rod insertion branches.

• Zero xenon and control rod insertion branches

Also, for Serpent2 cross section data simulations, PHOENIX5-based CBH-tables were applied, as these lattice calculations are not currently supported by the Serpent2/POLCA8 methodology.

(25)

Chapter 3

Specification of benchmark problems

This chapter will discuss the specification and implementation of the benchmark problem considered in this work. Besides, a complete description of the simulations carried out will be provided, as well as how some important parameters, such as the number of neutrons per cycle or the number of inactive cycles, were chosen before running the reference simulation. A short discussion about the most important physical implications of the CBH effect will also be given including the definitions of the parameters of interest for this project.

3.1 Description of the configurations

3.1.1 Benchmark configuration

The configuration used in this benchmark is a 3D 2x2 mini core, with three different fuel lattices in the axial direction. Each lattice has its enrichment of fissile materials and percentage of burnable poison. The fuel assembly lattices are based on a SVEA-09 Optima3 model for the Olkiluoto nuclear power plant, unit 2, residing in Finland but with some simplifications. Consequently, and to facilitate the Serpent2 calculations, the following assumptions are made:

• Water diamond and water cross are replaced with water pins

• A regular fuel pin pitch is used to facilitate the calculation of pin power

A radial representation of the mini-core system at the bottom of the core is shown by the Figure 3.1. In this geometrical representation, the cruciform control rod is situated in the northwest corner of the problem domain.

As previously said, each fuel assembly lattice type in the axial direction has its own material composition, described here from the bottom to the top of the core:

• Lattice 217: 4.15 w/o²³⁵U, 8x6.00 & 2x3.00 w/o Gd₂O₃

• Lattice 219: 4.19 w/o²³⁵U, 8x6.00 & 4x3.00 w/o Gd2O3

• Lattice 221: 4.15 w/o²³⁵U, 8x6.00 & 4x3.00 w/o Gd₂O₃

Lattice 217 has a fewer amount of burnable poisons rods compared to the other two lattices. Therefore it is expected to get a higher power production in the bottom part of the core than in the upper part.

To follow the conventional apporoach of nodalization normally employed in nodal core analysis, and also to facilitate data processing after simulations, the core is subdivided into 25 axial nodes, as presented in Figure 3.2,

11

(26)

Figure 3.1: Fuel assemblies at the bottom of the core

where the ++ present the position of the control rod inside the core. This is also the spatial mesh used axially by the flux solver of POLCA8. This figure also demonstrates the axial composition of the core, designed by the numbers 217, 219, and 221, and the nodes where the control rod is present when they are half-inserted into the core.

Figure 3.2: Axial subdivision of the core

This mini-core is depleted up to 10 MWd/kg starting at fresh conditions, with steps of 0.250 MWd/kg. The control rod is half inserted up to 5.000 MWd/kg, and then withdrawn from the core. To obtain data at the burnup

(27)

point of 5.000 MWd/kg for both the rodded and unrodded condition, a small step at 4.999 MWd/kg was added of to the POLCA8 and Serpent2 simulations.

Neutrons are kept inside the core in the radial direction thanks to the reflective boundary conditions. How- ever, due to vacuum boundary conditions applied at the top and bottom edge of the core, a strong neutron flux gradient and leakage will be induced in the axial direction of the core.

Inside a BWR core, coolant steam will be created inside the core. This leads to the creation of a void profile between the bottom and the top of the core, which manifests itself as a change in coolant density, as shown in Figure 3.3. Therefore, and to resemble real operation conditions as close as possible, a fixed axial coolant density profile was imposed over the core. Also, fixed values of the moderator temperature, 559 K, and the fuel Doppler temperature, 900 K, were applied in this benchmark problem. The xenon was set to equilibrium conditions.

Figure 3.3: The imposed axial coolant density profile in POLCA8 and Serpent2 calculations

In order to avoid any thermal-hydraulic feedback in the POLCA8 calculations (i.e., enabling the use of the above-fixed values of relevant state parameters), the thermal-hydraulics calculations in POLCA8 were deactivated. In this way, unambiguous comparisons of results between Serpent2 and POLCA8 can be established.

3.1.2 The Control Blade History (CBH) effect

Inside a BWR, the control rod is inserted into the inter-assembly gaps that are normally filled with water.

They are also often inserted during a long irradiation period. These two conditions combined lead to tilted distributions of power or flux compared to the unrodded situation. Subsequent withdrawal of the control rod after such a rodded irradiation period will in turn cause a strong power peaking in the neighborhood of the control rod blades, an effect called the CBH effect.

(28)

In the vicinity of the control rod, two physical phenomena create an accumulation of²³⁵U and²³⁹Pu thanks to a harder neutron spectrum. Due to the large thermal absorption cross-section of the control rod absorber material in a combination of water displacement decreasing the moderation properties locally close to the control rod, fewer²³⁵U nuclides can undergo fission, and therefore retardation in the depletion of this fissile material appears. Also, due to the harder spectrum locally, a larger number of neutrons are captured in

238U resonances, thereby enhancing the breeding of²³⁹Pu in this region. Therefore, when the control rod is withdrawn, a strong spectrum softening will occur due to the removal of the control rod absorber material and due to the increased amount of water appearing in this area. This will in turn lead to a severe pin power peaking in the vicinity of the control rod [16].

Due to this effect, a special correction term is needed for the calculation of the macroscopic cross-sections, the flux discontinuity factors, the detector constants, and the pin power form function in POLCA8. Capturing this CBH power peaking is essential from a safety point of view when evaluating thermal loads of fuel pins against prescribed design limits.

3.2 Setting up the reference solution in Serpent2

To obtain confidence in the Serpent2 reference solution, several issues must be addressed before running the final simulation. This includes the inter-cycle bias of the eigenvalue due to the use of the Power Method to obtain the eigenvalue solution, the convergence of the fission source, and how to implement an axial void profile over the reactor core domain in Serpent2. These items will now be explained in more detail.

3.2.1 Bias in the multiplication factor

In Monte Carlo simulations, the power method iteration process is employed to obtain the eigenvalue solution if the problem considered is a transport problem. With the use of the power method comes also an inter-cycle bias in the computed eigenvalue. It appears that this bias does not depend on the number of cycles used during the simulation but on the inverse of the number of neutrons used in each cycle. Therefore, this bias can be decreased by increasing the number of neutrons per cycle [18]. On the other hand, this process requires a lot of computational time. Before running the reference simulation, the optimal number of neutrons per cycle was sought: a trade-off between accuracy and computing time.

Table 3.1 shows the evolution of the eigenvalue bias directly given by Serpent2 Monte Carlo code, with the number of neutrons per cycle. The five simulations shown were carried out by using the same Serpent2 model of the 3D mini-core, with the same number of active and inactive cycles, 300 and 300, respectively. The use of a large number of inactive cycles ensures convergence of the fission source distribution before performing the calculation.

By also comparing the results of these simulations to previously evaluated configurations considering an unrodded single-channel system (i.e., with the size of one fouth of the current system) where only 1 million neutrons per cycle were used [19], using 4 million neutrons per cycle was found appropriate for this analysis.

(29)

Table 3.1: Bias in the eigenvalue for the evaluated 3D mini-core Neutrons per cycle Eigenvalue Bias (pcm)

1.10⁶ 1.01407 7.5

2.10⁶ 1.03919 4.7

3.10⁶ 1.04360 4.2

4.10⁶ 1.03316 4.0

3.2.2 Convergence of the fission source

In parallel with the determination of the number of neutrons per cycle, the convergence rate of the fission source (i.e., an indirect measure of the eigenvector) was studied. The convergence of the fission source distribution is determined by the power iteration process. In this process, both the eigenvalue and the fission source distribution converge to a certain value. The rate of these convergences depends on the dominance ratio of the system, which is the ratio between the second and the first eigenvalue of the system. In general, the eigenvalue converges faster than the corresponding eigenvector and in systems with a high dominance ratio, this effect is more pronounced. Therefore, to ensure the convergence of the fission source and the corresponding eigenvector (i.e., the neutron flux), enough number of inactive cycles need to be deployed before tallying any desired statistics [18] [20]. One way to obtain the convergence of the neutron flux is by looking at the Shannon entropy, which is a well-known concept from information theory. The Shannon entropy converges to a single steady-state value as the source distribution approaches stationarity [20]. Consequently, in order to find the optimum number of inactive cycles, different simulations (using the same number of active cycles) were performed. These results are presented in Figure 3.5 and in Figure 3.4.

Figure 3.4: Stabilisation of the Shannon entropy as function of the number of inactive cycles

(30)

Figure 3.5: Convergence of the eigenvalue as function of the number of inactive cycles

As Figure 3.4 shown for the stabilization of the Shannon Entropy, one can observe that obtaining the stable value is delayed by the case with 4 million neutrons per cycle, around 150 cycles, in contrast to the others at 50 cycles. The same trends are seen for the eigenvalue itself (Figure 3.5). After these first 150 cycles, the entropy is still decreasing weakly until 300 cycles and beyond, but this decrease is considered rather insignificant and 300 inactive cycles are considered enough to ensure the convergence of the fission source at the first stage of the simulation. The choice to use 4 million neutrons per cycle is also confirmed in this study because the corresponding relative statistical error of the Shannon entropy (given by Serpent itself) is significantly less important when using 4 million neutrons than in other cases.

Serpent2 has the possibility to reduce the number of inactive cycles used during a burnup calculation.

When the feature is applied, the fission source distribution at the end of one transport calculation is used as the initial source distribution for the next transport calculation. Therefore, the number of inactive cycles used to accomplish convergence of the eigenvalue and the eigenvector can be reduced as it is expected that the source distribution will not significantly change between two burnup steps. Reducing the number of inactive cycles in all but the first burnup step will contribute to the reduction of the total computing time. By looking at both figures (Figure 3.4 and Figure 3.5), 100 inactive cycles are considered enough to ensure the convergence of the fission source distribution during the depletion calculations.

Finally, the choice of using 300 active cycles for this problem was based mainly on earlier experience with an unrodded single-channel system demonstrating the appropriateness of using a few hundreds of active cycles to obtain uncertainties below 10 pcm [19].

(31)

3.2.3 Void Profile

As previously mentioned, a predefined axial coolant density distribution typical for BWR core operation has been imposed over the considered core. In Serpent2, it is possible to change the density of a material using a so-called multi-physics interface file (ifc). In this file, the user may specify the material zones where these changes are applied, and also choose the number of axial bins (i.e., equivalent to axial nodes in POLCA8 model). For each axial bin, values of the coolant density and temperature are provided in this interface file. A detailed description of this multi-physics interface file is found in the Serpent wiki website [17]. To apply the same coolant density and temperature profiles in the POLCA8 simulations, manually edited distribution files with the relevant values for these parameters need to be provided to POLCA8. Also, any thermal-hydraulics feedback needs to be deactivated.

3.3 Quantities of interest and statistics

The most important parameters to be considered in this project are the core eigenvalue and the power distribution at different levels of resolution, i.e., on an assembly node and pin basis. The results of POLCA8 simulations are directly compared, after converting them into the correct unit, to the ones of the Serpent2 reference simulation. All quantities of interest (i.e., eigenvalues, power distributions, etc.) are directly computed by POLCA8 and provided in the corresponding simulator output files. This is in contrast to Serpent2, where only the core eigenvalue is directly obtained in the output file. Consequently, in order to obtain the nodal pin power distribution, the most detailed level of resolution for powers, a mesh detector is used. All other quantities of lower resolution are then derived from these detailed pin-based powers, i.e., the nodal power distribution as well as the assembly radial and axial power distributions.

After obtaining the numerical results for both simulations (Serpent2 and POLCA8), one relevant way to compare their eigenvalue is by using the following definition for the eigenvalue error:

ε_k[pcm] = k − k^{re f} × 10⁵ (3.1)

The corresponding error in powers is given by the following expression:

εpi[%] =

p_i− p^{re f}_i

× 10² (3.2)

It should be recalled that by using an absolute formulation of these power errors in contrast to a relative one (i.e., not dividing the difference by the reference value), an undesirable amplification of the errors in low power regions will be avoided.

3.4 Normalization of the results

Before applying the comparisons between POLCA8 simulations and the 3D Serpent2 reference simulation described in section 3.3, results coming from the different power distributions are normalized according to the described process.

The most detailed level of power distributions, the pin power distributions, are normalized as follow:

∑

i

p^rel_i

N_f = P_n^rel (3.3)

(32)

where p^rel_i is the relative power density of the pin i, P_n^rel is the relative power density at the node n, and N_f is the number of fuel pin at node n. The summation is performed at one precise axial level.

Then, the axial power distribution over one bundle is also normalized as follow:

25 n=1

∑

P_n

P_tot = 1.0 (3.4)

where P_nis the power produced, in Watts, at axial node n and P_tot is the total power produced, in Watts, by the bundle.

Finally, radial power distributions over the mini-core can also be normalized as follow:

4

∑

b=1

P_b^rel

4 = 1.0 (3.5)

where P_b^relis the relative power density at the bundle b.

It should be recall that POLCA8 already returns the results as relative power density.

(33)

Chapter 4

Results

This chapter presents the analysis of obtained numerical results. All results are processed using a Matlab script, designed for this master thesis. A description of this analysis tool is given in Appendix A. The analysis is divided into two different sections, each focusing on a particular aspect of the conducted simulations. The first part consists of a complete analysis at the very beginning of the simulation (i.e., at the zero cycle burnup), to assess and quantify errors induced due to assumptions employed in the two-step core analysis methodology without consideration of any history effects (see the discussion in Section 1.1). The second part considers depleted core conditions to assess the impact of leakage history and CBH on the prediction accuracy of the two-step nodal core analysis methodology. As previously mentioned in Section 2.4 about the cross-section generation, the CBH-tables for Serpent2 are taken from PHOENIX5-calculations due to the lack of such branch capability in the Serpent2/POLCA8 methodology.

Due to the computing power required for other simulations at Westinghouse, it was not possible to run the reference simulation up to the selected burnup point of 10.000 MWd/kg. The simulation was carried out up to a burnup of 9.250 MWd/kg. However, this will not change the outcomes and conclusions of the analysis because the challenging burnup points with regard to CBH at the time point of control rod withdrawal have already been successfully evaluated.

4.1 Analysis at the fresh core conditions

4.1.1 Core multiplication factor predictions

One core condition of interest in this study to be analyzed is the first depletion point, i.e., at fresh core conditions. As the CBH effect occurs during burnup of materials, the results of POLCA8 with the CBH model activated or inactivated should by definition be the same. It should also be recalled that any differences or errors observed in predicted parameters at fresh core conditions will accumulate with time and will be either amplify or diminish depending on the ability of these errors to heal with burnup.

One of the first parameters to look at is the core eigenvalue, presented in Table 4.1 and in Table 4.2.

Table 4.1 presents comparisons of the core eigenvalue between all simulations, as well as the differences, expressed in pcm, between POLCA8 simulations and the reference simulation, using both PHOENIX5 and Serpent2 generated nodal cross section data. In Table 4.2, various reactivity components are shown using either PHOENIX5 or Serpent2 as the lattice physics code in combination with POLCA8.

19