Simplified Simulator for BWR Instabilities

Simplified Simulator for BWR Instabilities

ISTVÁN SOMA KOVÁCS

Master’s Thesis at the Department of Reactor Technology, KTH Supervisor: Jan Dufek

Examiner: Henryk Anglart

TRITA FYS 2017-41

iii

Abstract

The main idea behind this thesis project was to develop a demonstration and visualization tool for the Reactor Dynamics and Stability course at KTH, which can further improve the quality of the course. This means the devel- opment of both a mathematical model and a simulation program, which are able to reproduce certain phenomena related to BWR dynamics and stability, specifically the global and regional oscillations.

Due to their relative simplicity, several simulation programs have been developed to illustrate the dynamics of Pressurized Water Reactors (PWRs).

However, these kind of instability events (mainly due to void feedback) mostly occur in Boiling Water Reactors (BWRs), which are much more difficult to model from a thermal-hydraulic perspective, considering the multi-phase flow present in the core.

Therefore the usual PWR neutron-kinetic model could not be used for our purposes, and a new coupled neutron kinetic and thermal hydraulic mathemat- ical model was developed based on solving the neutron diffusion equation with the finite differences method numerically. This mathematical model was then translated into Python code, using Python 3.5.

The first simulation aimed to reproduce the global and regional instability events and serve as a reference for the future simulations. The objective of the second round was to establish how the reactor power and mass flow affects the longevity and intensity of the oscillations. The third and final round was con- ducted to measure the decay ratio of an induced global oscillation at different operational points and produce a stability map for the model reactor based on the results. Animations of the reactor transients were also recorded both in stable and unstable states, providing excellent demonstration material for the course.

Acknowledgement

I am forever indebted to all the fantastic people who had helped me during my thesis project. First of all, I would like to thank my supervisors Jan Dufek and Henryk Anglart for their guidance and motivation which helped me to finish the project.

Secondly, I would like to thank Vasily Arzhanov, who helped me a great deal during the deveopment of the neutron diffusion solver. Without him, I would be probably still working on it today.

My dear colleagues and flatmates Zsolt, Valerio, Salvatore, Olof and Simon also have my gratitude for being there for me and helping me out whenever I got stuck on something, and for providing me with a productive atmosphere.

Lastly, I want to thank my family for making all of this possible and supporting me throughout these long university years, never doubting me.

Kaposvár, 2016.06.25.

István Soma Kovács

v

Contents vi

1 Motivation and Objectives 1

2 Background Theory 3

2.1 The Neutron Diffusion Equation . . . 3

2.2 Finite Differences Method and Nodalization . . . 5

2.3 Core Thermal Hydraulics . . . 10

2.4 Basics of Reactor Dynamics and Reactivity Feedback . . . 14

3 Model Description 19 3.1 General Structure of the Coupled Kinetic and T-H model . . . 19

3.2 Steady State Iteration Scheme . . . 20

3.3 Transient Calculations . . . 23

3.4 Pressure Drop in Transient Flow . . . 25

4 Results and Discussion 29 4.1 Simulation Objectives and Parameters . . . 29

4.2 The Effect of Increased Reactor Power . . . 31

4.3 The Effect of Reduced Mass Flow Rate . . . 32

4.4 Stability Map . . . 33

5 Conclusions 37

6 Appendix 40

Bibliography 49

vi

Chapter 1

Motivation and Objectives

Good visualization and demonstration are both important elements of a high quality university course. There are many areas in engineering which are difficult to under- stand through purely studying systems of equations and mathematical interactions, and yet they can be shown and demonstrated graphically in a really educative and simple way. Nuclear reactor dynamics and stability and dynamics is definitely one of these fields of study, where sometimes one can feel overwhelmed by the number of connections and interactions between different elements of the system, and therefore any good illustration tool is greatly appreciated.

The main motivation behind this thesis project is exactly such an endeavour:

to develop and produce a software with the purpose of demonstrating and visual- izing certain stability related phenomena for the Reactor Dynamics and Stability course at KTH, greatly improving the course experience for all parties. Naturally, program development has endless possibilities, therefore clear goals had to be set before starting the project:

• the program should be able to simulate spatial effects, and therefore two dimensional model is needed at least,

• both neutron kinetic and thermal-hydraulic phenomena should be modelled with sufficient accuracy,

• the mathematical model and its solution should be simple and efficient to minimize the calculation time,

• not only a steady-state, but a transient model is needed to investigate the model’s response to certain instabilities, especially the so-called global (in- phase) and regional (out of phase) oscillations,

• the results of calculations should be recorded, stored with several visualization options for later use.

1

Chapter 2

Background Theory

2.1 The Neutron Diffusion Equation

This section introduces the one-speed neutron diffusion model and derives the one- speed Neutron Diffusion Equation (NDE). The NDE plays a very fundamental role in the development of reactor analysis tools. It is a great tool for solving neutron physical problems in complex systems where the transport equation would be ex- tremely difficult or costly to solve.

The equation is derived and explained thoroughly in a textbook by Duderstadt.

A shortened version of his description is presented in this section, consequently all statements and derivations in this section are taken from Duderstadt (1976), unless otherwise stated. For the purposes of the thesis project, the one-speed diffusion theory model was used, since according to Duderstadt "it is sufficiently sufficiently simple to allow detailed calculations and also sufficiently realistic to allow us to study many of the more important concepts arising in nuclear reactor analysis".

The neutron diffusion model uses neutron density N(r,t), which gives the num- ber of neutrons per unit volume at position r and time t. Sometimes it is more convenient to use the neutron flux, which is defined as φ = vN(r,t), where v is the average speed of the neutrons.

Neutron flux is very useful when it comes to calculating reaction rates as well, since in a way, it represents the total distance travelled by neutrons in a unit vol- ume per unit of time. Since the unit of reaction cross-sections is 1/cm or barn, multiplying the respective cross-section with the local flux gives the local reaction rate for the specific nuclear reaction.

There are generally three types of reactions we have to consider when it comes to one-group neutron diffusion: diffusion in and out of a volume, absorption inside the volume and neutron production. Therefore the one-speed neutron diffusion ap-

3

proximation takes the following form:

1 v

∂φ

∂t = −∇J − Σ_aφ + S (2.1)

where J is the net neutron current through the boundary surface, Σ_ais the macro- scopic cross-section for neutron absorption, and S is the so-called source term. We have two unknowns in this equations, since we do not know neither φ(r, t) nor J(r, t). There is no exact relationship between the two, but it is well known from other fields that the current density is approximately proportional to the negative gradient of density, or in our case, the flux. This is known as Fick’s law (Hetrick 1971).

J(r, t) ≈ −D(r)∇φ(r, t) (2.2)

where D(r) is the so-called diffusion coefficient. Equation 2.2 is often called Fick’s law. If we substitute Eq. 2.2 into Eq. 2.1, we arrive at the one-speed neutron diffusion equation:

1 v

∂φ

∂t = ∇D(r)∇φ(r, t) − Σ_aφ(r, t) + S(r, t) (2.3)

2.2. FINITE DIFFERENCES METHOD AND NODALIZATION 5

2.2 Finite Differences Method and Nodalization

This section describes the finite differences method for numerical calculations, and in particular its uses in solving the Neutron Diffusion Equation(NDE). The neutronic module of the simulator program is based on the NDE, which can be used to determine the neutron flux distribution in the core. For a realistic neutron diffusion calculation, the heterogeneous nature of the nuclear reactor core has to be considered. Not only there are several different kinds of materials in the core (fuel, cladding, spacers, coolant, control rods), but also there are non-uniformities due to power, temperature and density differences (Duderstadt 1976).

These non-uniformities mean that an analytical solution to the NDE is much less appealing than a direct numerical solution. The general procedure for solving the NDE numerically is to rewrite the equation in finite difference form, and then solve the resulting system of equations. In his textbook, Duderstadt also proposes a way how to efficiently perform these steps, which has been widely used in the past decades, therefore our solution method was also largely based on his, so all derivations are taken to be from his textbook (Duderstadt 1976), unless specified otherwise.

The first step is to properly divide our geometry (be it 1D, 2D or even 3D) into multiple nodes or cells. Inside these cells, all material properties are considered uniform for each point in time, but they can vary from cell to cell. To perform a detailed calculation, a multidimensional geometry is needed (usually 2D or 3D with rectangular cells). For the purposes of the thesis project, a 2D geometry was used, therefore in this section specifically the 2D discretization and solution of the NDE is derived. Initially, we seek the solution for the steady-state NDE, by using Eq.

2.3 as if its left side was zero.

S(r) = −∇D(r)∇φ(r) + Σ_a(r)φ(r) (2.4) This means that neutron production is in balance with the losses, practically in each cell the number of neutrons leaving the cell plus the number of neutrons absorbed inside the cell is equal to the number or neutrons produced. The same equation also has an integral form, where the incoming flux, the absorption and the production is integrated over the specific cell’s area (A_i,j) (Arzhanov 2017):

− D_i,j Z Z

Ai,j

∂²φ

∂x² dA − D_i,j Z Z

Ai,j

∂²φ

∂y² dA + Z Z

Ai,j

Σ_aφ(x, y) dA = Z Z

Ai,j

S(x, y) dA (2.5)

To simplify the equation, uniform rectangular discretization can be assumed, where all cells have the same dimensions (h_x and h_y). For each cell, we get five- point equations: In other words, the flux in each cell is influenced by the flux of the four neighbouring cells and the diffusion coefficients in each direction, as well as the

absorption and production in the specific cell itself. This five-point connection grid is shown on Fig. 2.1.

Figure 2.1. The 5-point stencil

Finite difference approximations can be used to solve the integrals in Eq. 2.5, which yields the following form (Arzhanov 2017):

− 2h_yDi,jDi+1,j

hxDi,j + h_xDi+1,j

(φ_i+1,j − φ_i,j) + 2h_yDi−1,jDi−1,j

hxDi,j+ h_xDi−1,j

(φ_i,j− φ_i−1,j)−

− 2h_xDi,jDi,j+1

h_yD_i,j+ h_yD_i,j+1(φ_i,j+1− φ_i,j) + 2h_xDi,jDi,j−1

h_yD_i,j+ h_yD_i,j−1(φ_i,j− φ_i,j−1)+

+ h_xh_yΣ_a,i,jφ_i,j = h_xh_yS_i,j (2.6)

More on information about the application of the finite differences method can be found in a textbook by Quarteroni (2000). To achieve a more compact form, both sides are divided by h_xh_y and some new notations are introduced:

2.2. FINITE DIFFERENCES METHOD AND NODALIZATION 7

− 2D_i,jD_i+1,j

h²_x(D_i,j+ D_i+1,j) = a_i+1,j

− 2D_i,jDi−1,j

h²_x(D_i,j+ D_i−1,j) = a_i−1,j

− 2D_i,jDi,j+1

h²_y(D_i,j+ D_i,j+1) = a_i,j+1

− 2D_i,jDi,j−1

h²_y(D_i,j+ D_i,j−1) = a_i,j−1

Substituting these new terms back into Eq. 2.6 and dividing both sides with hxhy, we get the following form:

φi,j

Σ_a,i,j− a_i+1,j − a_i−1,j− a_i,j+1− a_i,j−1^+ φ_i+1,jai+1,j+

+ φ_i−1,ja_i−1,j+ φ_i,j+1a_i,j+1+ φ_i,j+1a_i,j+1= S_i,j (2.7) This form is quite expressive and therefore relatively easy to understand. Basi- cally the terms inside the square bracket represent the absorption of neutrons inside the cell and leakage from the cell to the neighbouring ones. The φ_ka_k terms rep- resent the incoming neutron leakage from the respective neighbouring cells, while S(i, j) represents the number of neutrons produced in the cell. To further simplify the equation, we introduce the so-called central term:

ci,j = Σ_a,i,j − a_i+1,j− a_i−1,j− a_i,j+1− a_i,j−1 (2.8) So far our model only considered prompt neutrons, but for a more realistic calculation, delayed neutrons and their precursors also need to be accounted for.

The decay of delayed neutron precursors is simply added to our source term, while we get Eq. 2.10, which describes the decay and production of precursor nuclei (Anglart 2011):

φi,jci,j+ φ_i+1,jai+1,j+ φ_i−1,jai−1,j + φ_i,j+1ai,j+1+ φ_i,j+1ai,j+1=

(1 − β)νΣ_f,i,jφ_i,j+ λC_i,j (2.9)

βνΣi,jφi,j = λC_i,j (2.10)

To find the steady-state solution for this 2D system, these two equations must be solved for all points. To conveniently achieve this, we can organize the different terms into vectors and matrices. We define φ as a column vector, which contains the flux values for every point, i, j double indices are converted into a third index:

k = (N − 1)(j − 1) + i, where N is the number of cells along the x-axis. Σ_a

absorption cross-section, Σ_f fission cross-section, C precursor concentration, D diffusion coefficient vectors are defined in a similar fashion. The leakage terms and the central term can be organized into a matrix, which contains our five-point equations for every point (Duderstadt 1976):

A =



































































c_1,1 a_1,2 0 0 0 · · · a_1,N · · · 0

a_2,1 c_2,2 a_2,3 0 0 . .. . .. a_2,N . .. ... 0 a_3,2 c_3,3 a_3,4 0 . .. . .. . .. a_3,N ... 0 0 a4,3 c4,4 a4,5 . .. a4,N . .. . .. . ..

... . .. ... . .. . .. . .. . .. . .. . .. ... ... . .. ... . .. . .. . .. . .. . .. . .. ... a_k,k−N . .. ... a_k,k−1 c_k,k a_k,k+1 . .. a_k,k+N . .. ... ... . .. ... . .. . .. . .. . .. . .. . .. ... ... . .. ... . .. . .. . .. . .. . .. . .. ... ... . .. ... . .. . .. . .. . .. . .. . .. ... ... . .. ... . .. . .. . .. . .. . .. . .. ... 0 · · · aZ,Z−N . .. . .. aZ,Z−1 cZ,Z



































































Matrix A is always a diagonal, k × k matrix, where k is the number of cells in each direction. The central terms are in the main diagonal for each cell, and there are 4 sub-diagonals which contain the leakage terms for all four directions.

The general structure of the matrix and the mesh is shown on Fig. 2.2. Naturally, if said cell is at one of the boundaries, then in that direction a specific leakage term has to be calculated, based on the local boundary condition. Putting all the different terms into one matrix simplifies the solution of the diffusion problem into a linear matrix equation (Duderstadt 1976):

A × φ = S (2.11)

2.2. FINITE DIFFERENCES METHOD AND NODALIZATION 9

Figure 2.2. The FD mesh and the structure of matrix A

2.3 Core Thermal Hydraulics

In this section the basics of reactor thermal-hydraulics and some related phenom- ena are discussed: the so-called Homogeneous Equilibrium Model, as well as several simplifications and generalizations used in our model. These have been described by many people in the past, a very good source of information on the topic is the Nu- clear Reactor Technology course compendium by Henryk Anglart. Several selected concepts from the book are introduced in this section to provide a solid basis, from which the thermal-hydraulic module of the program was developed. All statements and derivations are taken from (Anglart 2016), unless stated otherwise.

The BWR reactor is the simplest way to harness fission energy and convert it to electricity, but this is not entirely true when one tries to describe the thermal- hydraulics of the BWR. Very simply, the coolant water (which also acts as the moderator) enters the core subcooled, at the inlet at the bottom of the core. Nu- clear fission produces heat in the fuel rods, which is transferred to the coolant. As a result, the water is gradually heated up as it flows through the fuel assemblies.

Assuming there is no reverse flow in the core, the enthalpy of the coolant is in- crementally increased moving from the inlet towards the outlet. Eventually, the water reaches saturation and starts boiling, and by the time it reaches the outlet, a considerable fraction of the initial liquid is converted into steam. Then electricity is produced from the steam using steam turbines (Lewis 2008) .

q"

z

z z=z+dz

W,i_li i_l(z)+di

il(z)

Figure 2.3. A heated channel. Source: (Anglart 2016)

Since power is non-uniformly produced in the reactor, the resulting outlet en- thalpies and properties depend on the radial and azimuthal position of the chan- nel. Additionally, the fuel assemblies usually have a solid wall, which prevents the coolant from mixing between different channels. Each assembly can be considered a separate one-dimensional channel, where there is minimal contact or interaction between the different channels. Such a channel is shown on Fig. 2.3.

2.3. CORE THERMAL HYDRAULICS 11

Let us assume a heated channel with an arbitrary axial distribution of heat flux, q⁰(z), and an arbitrary, axially-dependent geometry, as shown on Fig. 2.3. The coolant flowing in the channel has a constant mass flow rate W . Energy balance for a differential channel length between z and z + dz is as follows:

W il(z) + q”(z)P_H(z)dz = W [(z) + di_l] (2.12)

where q” is the heat flux , P_H is the heated perimeter of the channel and i_l(z) is the specific enthalpy of the coolant. This leads to the following differential equation:

di_l(z)

dz = q”(z)PH(z)

W (2.13)

If we integrate Eq. 2.13 from the inlet of the channel to elevation level z, we get:

i_l(z) = i_l,i+ 1 W

Zz

−H 2

q”(z)P_H(z)dz (2.14)

From Eq. 2.14, we can calculate the enthalpy distribution for each channel re- spectively by either finding the analytical solution or by using numerical integration.

Knowing the enthalpies is only the first step however, other factors of interest are the thermodynamic quality, void fraction and the pressure loss along the different channels. These properties are directly linked to the different feedback mechanisms, which are discussed in the next section.

The Homogeneous Equilibrium model (HEM) is the simplest two-phase flow model, which is built on the assumption that both the liquid and gaseous phases move with the same velocity and they are in a thermodynamic equilibrium with each other. Therefore the local, channel-average void fraction is determined from the corresponding local value of the equilibrium thermodynamic quality. The HEM expression for the void fraction (α) takes the following form:

α =















0 xe≤ 0

1 1 +ρg

ρ_l 1 − x_e

x_e

0 < x_e< 1

1 x_e> 1

(2.15)

Calculating the void fractions for every point makes it possible to get the pres- sure losses in each channel. Pressure loss values are of grave importance from a stability aspect since they influence the mass flow distribution in the channels. In

a BWR channel, the water coolant usually enters subcooled at the inlet. As the water starts heating up, and its enthalpy gradually increases, it reaches a point in the channel, where boiling first appears.

This is called the Onset of Nucleate Boiling (ONB) point. Since the pressure loss is calculated differently for single-phase and two-phase flows, it is essential to know the location of this point in each channel. If subcooled boiling is neglected, one can assume that up to the ONB point, there is a single-phase flow, while after the ONB point a two-phase flow can be observed in the channel.

In general, there are several mechanisms which cause pressure drop along the channel:

• Friction losses from the fuel and control rods (walls)

• Local losses from spacer grids

• Local losses at the core inlet and exit

• Elevation pressure drop

The total pressure drop in a channel with a constant cross-section area can be calculated from the following equation:

− ∆p_tot = −∆p_{f ric}− ∆p_loc− ∆p_elev= 4C_fL D_h +^X

i

ξi

!G²

2ρ + Lρg sin ϕ (2.16) where C_f is the Fanning friction coefficient, L is the length of the channel, G is the mass flux, D_h is the channel hydraulic diameter, ξ_i is the local pressure loss coefficient and ρ is the coolant density. This equation can be used to calculate the pressure drop in the channel up to the ONB point. From then on, the two-phase flow pressure drop is calculated as:

− ∆p = r₂G²

ρ_l + r₃Cf,lo

4H D_h

G²

2ρ_l + r₄Lρlg sin ϕ +

N

X

i=1

Φ²_lo,iξi

!G²

2ρ_l (2.17) Here r₂, r₃ and r₄ are integral two-phase pressure drop multipliers (acceleration, friction and gravitation, respectively). Φ²_lo,i, ξ_i are local the local loss multiplier and the local pressure loss coefficient, respectively, at location i in the channel. The Fanning friction factor C_f,lo is found the same way as for single-phase flow applica- tions.

2.3. CORE THERMAL HYDRAULICS 13

What is important to note here, is the fact that the two-phase pressure drop sources are more or less the same as they are for single phase. However, there is a new effect to consider, i.e. the acceleration loss. r_i drop multipliers basically ad- just the single-phase loss terms to the two-phase values. These multipliers largely depend on the pressure, power distribution and their accurate value generally can only be found with numerical integration along the channels.

From Eq. 2.17 it is also evident, that the pressure drop mainly depends on two factors: the mass flow of the channel and its average thermodynamic quality. The larger the mass flow in the channel, the larger is the pressure drop. Water always flows towards the least resistance, therefore in an equilibrium, every channel has the same pressure drop in the core, although their exit qualities and mass flows may vary. The channels with higher exit quality are going to have a lower flow rate compared to those with lower exit quality.

2.4 Basics of Reactor Dynamics and Reactivity Feedback

This section gives a short introduction to the specific feedback loops present in a Boiling Water Reactor. Fig. 2.4 presents the block scheme of the most notable physical feedbacks. In a well designed reactor, these feedback effects help to control the reactor, but under special circumstances, they can also lead to instabilities.

While in PWR reactors there is virtually no boiling, and therefore the void fractions are negligible (and consequently oscillations seldom occur), in a BWR the coolant is boiling, and at the outlet the steam takes up the majority of the volume in the subchannels. Since void has a very strong effect on neutron moderation, it significantly influences the neutron dynamics of the reactor. There are two types of feedbacks to be mentioned: the first one is the reactivity void feedback which affect the reactivity directly, while the second, thermal hydraulic feedback only im- pacts the void fraction, and therefore affects the reactivity indirectly. An excellent introduction to these feedback mechanisms can be found in a textbook by B-G.

Bergdahl and others. All statements and figures in this section are taken from (B.- G. Bergdahl 2005), unless otherwise stated.

Figure 2.4. The different physical feedback loops in a BWR.

As it can be seen on Fig. 2.4, there are two main closed feedback loops in a BWR, which can lead to unstable oscillations. The red loop (void feedback) marks the void feedback, which influences reactor power directly. If this loop becomes un- stable during operation, it causes a so-called global oscillation in the reactor core, which means neutron flux oscillations in-phase in the complete core. The frequency of this global oscillation is the reactor resonance frequency, which is usually between 0.5-0.6 Hz. This is the most common type of instability in BWRs.

2.4. BASICS OF REACTOR DYNAMICS AND REACTIVITY FEEDBACK 15

Instability in the blue loop can also cause global oscillations. This loop corre- sponds to the so-called hydrodynamic feedback. Instability in this second loop is enough to cause power oscillations due to the void reactivity feedback, even if the red one is stable. This proves that the loops are not independent of each other and both of them need to be stable in order to achieve total stability. In many cases, this fact makes it difficult to decide which of the feedback loops is actually unstable, and to find out what caused the oscillation.

One can get a fairly good understanding on what happens exactly during a global oscillation by following the arrows around the red loop on Fig.2.4. A sudden jump in reactivity causes an increase in power and fuel temperature. This excess heat is conducted to the surface of the fuel rods, and is transferred to the coolant water, resulting in increased enthalpies and void production. This affects the reac- tivity in return, via the negative void feedback effect.

We can also see how the dark blue loop can cause global instability if we trace it around. The coolant is pumped from the lower plenum to the upper plenum through the core, and then it takes a turn into the downcomer. If there is a pump failure in the circuit, the coolant flow is reduced to natural circulation only. The driving force of the flow in this case is the pressure difference between top and bottom of the core. As the coolant passes through the core, it is heated up, and starts boiling, therefore its mean density decreases with the elevation. This builds up a pressure difference, which drives the flow, slowly increasing it. As the flow increases however, the pressure drop also starts to increase in the channels, and therefore the driving forces are eventually equalized. Therefore this is also a negative feedback loop.

It is also important to note here that there is no direct feedback from this loop to the nuclear power, so instabilities in this loop cannot be cancelled out with power regulation only. Oscillations in the hydrodynamic loop are much less likely to ap- pear however, since the feedback is stable as long as the forced cooling of the reactor core is upheld. Therefore if there is a global oscillations occurring, while there is forced cooling, it is always due to instability in the red loop.

Another type of instability worth mentioning is the so-called regional instability.

This type of oscillation can be caused by instability in the dark blue thermal hy- draulic feedback and the red void feedback loops together. It can be distinguished from the global oscillation by the fact that power starts to oscillate out of phase in the different parts of the reactor. The frequency is also around 0.5 Hz, similarly to the previous kind of instability, but here there is a 180⁰ phase difference between the two halves of the reactor. The danger of this instability event is that while power increases in one half of the core, it also decreases in the other half, so the total thermal power output of the reactor is potentially unchanged, and therefore it’s really hard to notice the oscillation from the operator room.

The initiating event of a regional oscillation can be a small fluctuation, when the mass flow slightly increases in one half of the core by a bit, and it simultaneously decreases by the same amount in the other region. This means that both pressure drops over the core and total coolant flow through the core are constant. In the first region, where the mass flow increased, this causes a reduction in void as well.

A reduction in void increases reactivity and power due to the negative feedback, which leads to more boiling and production of void.

Figure 2.5. The main phases of global (above) and regional (below) oscillations.

The void fraction is illustrated by bubbles.

In the second region, the process is very similar, but it starts from the opposite

2.4. BASICS OF REACTOR DYNAMICS AND REACTIVITY FEEDBACK 17

phase. The suddenly decreased flow results in higher void fractions, which lead to a decrease in reactor power. With decreased power, the void starts to also decrease, and then things start from where they did in the first region. It is very important that the total mass flow of coolant is kept constant at all times, but due to the changes in the void fractions, we get differences in flow rates between the two re- gions. Fig. 2.5 shows the main phases of the two kinds of oscillations mentioned before.

The difference between the global and regional oscillations is that a global one has to overcome the friction forces in the loop, while the regional one just needs to redistribute the flow, according to the different pressure losses in each channel.

Therefore the global instability is significantly damped by the thermal hydraulic feedback, but a regional one is strengthened by it. Neutron kinetics and the void feedback work the other way around however, they strongly dampen regional insta- bilities, while supporting global oscillations. In both cases, there are two competing effects, and either can prevail, depending on the design of the reactor and the op- erational point.

Figure 2.6. The decay ratio of an oscillation, and its calculation method.

An important characteristic of any particular oscillation is its decay ratio(DR).

The decay ratio by definition is the ratio between the amplitudes of two subsequent oscillations (Fig. 2.6). If DR is smaller than one, the oscillation is damped, and therefore the system is stable from a dynamics point of view. On the other hand, if the DR reaches one, the oscillation is sustained and it can cause serious damage in the core through material fatigue. If the DR is grater than one, the amplitude increases with every subsequent oscillation and the reactor is unstable. Therefore

the DR is a good indicator of stability and a key factor to look at when testing how different variables influence the behaviour of the reactor. (Anglart 2011)

The operators of a BWR usually have a so-called operational diagram for the reactor, which shows the decay ratio of the oscillations depending on the relative power and mass flow. This helps them to avoid potentially dangerous operational points, and guarantee the safe operation of the reactor (D’Auria 2008). One exam- ple operational diagram is shown on Fig. 2.7.

Figure 2.7. An example power-flow operational diagram for Leibstadt NPP.

Source:(D’Auria (2008))

As it can be seen on the figure 2.7, the most dangerous operational points are those with particularly high power-to-mass flow ratios. Therefore there are strict regulations and protocol how to avoid the unstable and the non-permitted opera- tional regions, especially during reactor startup and shut-down.

Chapter 3

Model Description

3.1 General Structure of the Coupled Kinetic and T-H model

The program created for this project is able to simulate the neutron physical and thermal-hydraulic behaviour of a cubical slab boiling water reactor. The main objective of creating the program was to help visualize the global and regional in- stability events, and make is easier to analyse the effect of the operational point on the evolution of these instabilities. For this purpose however, two dimensions were sufficient, and therefore the cube reactor was converted into a square (a = 1 m), which is reflected from both sides, but not from above and below.

In order to visualize a regional oscillation, at least a 2D solver was needed, since spatial effects are hard to observe in 1D. Some compromises and assumptions were made during development, since the program had to be as realistic and accurate as possible, while also performing the calculations relatively fast. Some of the assumptions and simplifications:

• As mentioned before, the reactor is reflected perfectly from all sides, except from its top and bottom, where we assumed vacuum boundary condition.

• Only a one-group neutron physical model was used, consequently all spectral effects are neglected.

• One group of delayed neutron precursors was used, therefore the different half-lives of delayed neutron precursors are also neglected.

• Mass flows are constant along a single channel, but they can obviously vary from channel-to channel (cells in the y-direction all have the same mass flow, but cells in the x-direction from each other can have different mass flow rates).

Pressure was also assumed constant for the sake of simplified water property calculations.

19

The reactor is divided into 50x50 two-dimensional cells, each with uniform 2x2 cm dimensions. Inside these cells, all properties such as enthalpy, void, neutron flux are considered uniform, but these properties can differ from cell-to-cell.

3.2 Steady State Iteration Scheme

The steady state solver is the tool which allows the program to iterate the initial values of the simulation based on the initial and boundary conditions specified by the user. In a simple neutron kinetic calculation one would only need to solve the steady-state NDE, but since in our model the neutronics are coupled with a thermal- hydraulic model, a complex iteration scheme is needed to ensure convergence to the correct solution. Our scheme is shown on Fig. 3.1, where we can see how the results of the n + 1^th iteration cycle are obtained using the results from the n^th step.

Closed loops can only be solved implicitly, so they all represent another inner iteration cycle. The first such loop is between the mas flow and the pressure drops, the second one is about how enthalpies and mass flow both influence each other.

The third, and most important loops is the one which gives us the new flux values based on the previous void distribution. This loop represents the void feedback described in the previous section.

Figure 3.1. The iteration scheme used to obtain the steady state solution. The different indices mark the results from different iteration cycles.

3.2. STEADY STATE ITERATION SCHEME 21

For the first equation, we must return to Eq. 2.11 in Chapter 2. If the source term is written out properly, we can see that the φ vector appears on both sides.

Consequently, the solution can only be found by iteration. To start the iteration, a starting S, Σ_a, Σ_f are needed, as well an initial flux-distribution, φ⁰. From the results of the n-th iteration, φⁿ⁺¹ is acquired by solving Eq. 3.1, while Cⁿ⁺¹ is calculated by solving Eq. 3.2.

Aⁿ× φⁿ⁺¹= 1

kef fνΣ_f× φⁿ+ λCⁿ (3.1) Cⁿ⁺¹= βνΣfφⁿ⁺¹

λ (3.2)

The next step on Fig. 3.1 is to calculate the resulting enthalpies for each cell from the flux-distribution of the reactor and the mass flow of the respective chan- nels in the previous step. This is done by solving a simple energy balance equation (Eq. 3.3), and then based on the enthalpies, qualities and void fractions can be calculated in accordance to the HEM model(Eq. 3.4).

iⁿ⁺¹_k = i_k−N+dy φⁿ⁺¹_k c_th

W_k (3.3)

α_k =















0 if x_e,k≤ 0

1 1 +ρ_g

ρl

1 − x_e,k xe,k

if 0 < x_e,k< 1

1 if xe,k> 1

(3.4)

where i_k−N is the enthalpy of the cell just below of the current one, c_th is the con- version ratio between neutron flux and thermal power and W_k is the mass flow in the channel.

Based on these results, the new absorption cross sections can be calculated for the whole core. The reactor was assumed to be undermoderated (i.e. the void feedback is negative). Since there is no need for high accuracy, a simplified linear model is used to translate the changes in void fraction (still calculated with the HEM model) to changes in the absorption cross-section. At the beginning of the iteration, each cell has the same base absorption value. The strength of absorption increases with the void, and in the model, the absorption cross section is assumed to be directly proportional to the void fraction:

Σⁿ⁺¹_a,k = K^αⁿ⁺¹_k + 1^Σ^base_a,k (3.5) where Σⁿ⁺¹_a,k is the absorption cross-section in the next step, K is an arbitrary con- stant, αⁿ⁺¹_k is the void fraction calculated from the current enthalpy values and

Σ^base_a,k is the base absorption cross-section.

After the new cross-section are found for each cell, the new matrix (Aⁿ⁺¹) can be constructed. The new pressure losses and mass flow can be iterated using Eq.

2.17 and the fact that in the end, the mass flows are redistributed in such a way, that each channel has the same pressure drop.

3.3. TRANSIENT CALCULATIONS 23

3.3 Transient Calculations

The steady-state iteration is only needed to obtain the necessary initial conditions for the simulation. After the solution has converged, the time-dependent simulation can be started. This is based on solving a number of time-dependent differential equations using the well-known Euler-formula (Eq.3.6).

y(t + ∆t) = y(t) + dy

dt dt (3.6)

The Euler integration is the simplest way to solve these equations, but it is also the fastest for the same reason. It also has a constraint in the fact that there is a maximum for the size of the time step, if we want to preserve mathematical stability (Quarteroni 2000). The limit for this system was found slightly above 10⁻⁴, so the time step was chosen to be 10⁻⁴seconds for the sake of convenienceArzhanov (2017).

Many of the equations used in the steady-state calculations were derived from the actual transient equations, assuming that the time-derivatives were zero. Con- sequently, the calculation steps are more or less similar, as it is shown on Fig. 3.1

The first task is to solve the transient NDE numerically. In order to increase the calculation speed, all equations were organized into a matrix, just like for the steady-state case. The key difference here is that the source term is no longer nec- essarily equal to the number of neutrons leaking out and getting absorbed, and therefore the time-derivative of flux can be non-zero. A new equation matrix (B) had to be constructed, which also contains the source term of each cell in the cell’s central term, while the leakage terms remained unchanged. Substituting matrix B into the the transient NDE yields the following form:

v (B × φ + λC) = dφ

dt (3.7)

This way, the left side of the equation can be calculated in a few easy steps, yielding a derivative-vector ^dφ_dt, which contains the time-derivatives of flux for each individual cell. Naturally, the delayed neutron precursor concentrations also change with time, their time-derivatives^dC_dt ^can be found by solving Eq. 3.8.

dC

dt = βν (φ × Σ_f) − λC (3.8)

Solving these two equations basically concludes the neutron kinetic part of the task, and only the thermal-hydraulic and feedback calculations remain. The simple energy balance equation which was used to obtain the steady-state enthalpies is no longer usable, since incoming and outgoing energy flows are no longer necessarily

equal in a cell. The flow directions also need to be taken into consideration, since the enthalpy of the cell just below our current one affects the enthalpy of the water entering the cell quite significantly. In our model reverse and transverse flow is not enabled, and therefore we can always be sure that the coolant water always enters a cell from below, and not from any other direction. The transient thermal hydraulic model is derived from Eq. 3.9 (M. Ishii 2011):

∂ (ρ_mi_m− p)

∂t + 1

A

∂ (W i_m)

∂y = q

A (3.9)

where the index m means mixture, referring to the mixture of saturated steam and liquid in a cell. A is the flow area of a channel, V is the mass flow, q is the heat flux on the channel wall, i is the enthalpy, ρ is the density and p is the pressure.

Assuming constant pressure in every cell and expressing the two partial derivatives using the finite difference method yields Eq. 3.10.

ρ_m,k di dt+ W

A

i_m,k− i_m,k−N dy = q_k

A (3.10)

From this, the time-derivative ^di_dt can be easily expressed and calculated. To minimize calculation time, ^di_dt vector was also defined, which contains each time- derivative. The equation to find ^di_dt takes the following shape:

di dt = Q

A +W A

i − i_low dy

dt

ρ (3.11)

where Q is the vector holding all current thermal power values, W is the vector containing the mass flow rates for all cells, while i represents the enthalpy values in each cell in the previous timestep. i_low is a special vector, for every k^th cell, it holds the enthalpy of the cell directly below i_k−N. For cells at the bottom of the reactor, this means the specified inlet enthalpy. The changes in void fractions and absorption cross-sections are calculated similarly as in the steady-state calculation (Eq. 3.4, Eq. 3.5).

3.4. PRESSURE DROP IN TRANSIENT FLOW 25

3.4 Pressure Drop in Transient Flow

As we saw in the previous section, some of the equations only needed to be extended by the time-derivative, and this slight modification allowed their use in the transient model. This is not the case however with the pressure loss and mass flow calculation, because it would require simplifications to such a degree, that it would render the whole reactor-dynamic model useless for the observation of dynamic instabilities.

The reason for this is the fact that the flow momentum of the channels cannot be disregarded. This means that if the water in the core moves with certain velocity at a given moment, the velocity in the next moment does not only depend on the momentary pressure loss, but also the previous velocity value as well. Additionally, changes close to the inlet only spread gradually towards the outlet, the speed of this depending on the flow velocity.

To include all these considerations, a completely new model was needed. Gener- ally, there are two kinds of models which can be used: analogue-type representations based on transfer-functions and finite differences numerical solutions of three ba- sic partial differential equations in space and time for continuity, momentum and energy balance(L. S. Tong (1997)). An example for the former is the HYKAMO code developed by Schoneberg in 1968 (Schoeneberg (1968)), an example for the numerical solution of the problem was presented by Turner and Trimble in 1976 (Trimble (1976)). Both these options have their pros and cons:

• The numerical solution is much more demanding from a computational as- pect, but it certainly yields accurate results.

• On the other hand though, an analogue representation of the pressure drop feedback would require orders of magnitude less computational capacity, but naturally it would also be much less accurate.

After much consideration, the choice was made to use an analogue represen- tation for two reasons: Firstly, the simulations were run on ordinary laptops and PCs, so it had to run at a reasonable speed. This meant that the model needed to be as efficient as possible, requiring the minimal amount of calculations. Secondly, since this was to be a demonstration and visualization tool, accuracy of the highest degree was not of paramount importance. This model was not to be representative of any real-world reactor and the calculation results would not be compared to the results of any tests. This meant that the analogue model was a more than accept- able compromise, giving us what we needed, while taking away what we did not really need in the first place.

In order to have a better understanding on what this analogue representation means, in this section, some new concepts are introduced from the field of automatic control theory. In a way, we can think of the pressure loss feedback as some sort of regulator circuit, which controls the flow in each channel, depending on certain flow parameters. A good example of such a control circuit is the so-called tempomat in cars, which controls the amount of fuel the engine receives, depending on how fast our car is moving. If the car is moving faster than the maximum allowed speed, the engine gets less and less fuel and the car starts decelerating until we reach the set speed (B.-G. Bergdahl 2005).

These control circuits can be easily represented by block diagrams, one such diagram on Fig. 3.2 shows an example negative feedback control circuit loop. Each block has a so-called transfer function, which determines how said block modifies the incoming signal. The input signal is always multiplied by the transfer function of all block it goes through. Consequently, the whole system has a transfer function as well, which is a superposition of all the individual block transfer functions. In our case, the output signal is then returned to the input, creating a closed feedback loop (B.-G. Bergdahl 2005).

Figure 3.2. A simple negative feedback control circuit, represented by a block diagram. Source: Anglart (2011)

An analogue representation of the pressure loss feedback is basically a closed feedback circuit which gives some kind of a response depending on the flow proper- ties (i.e. determine the direction and scale of the change in the pressure loss for the next time-step). Naturally, for the model to behave like a real BWR, this feedback loop had to preserve all important characteristics of the real reactors. These main characteristics were:

• Mass flow in each channel depends on the pressure loss of said channel, and flow is always distributed in such a way to make pressure losses equal.

• Any change in the flow is not immediate due to the flow momentum. Therefore the feedback system also needs to be incremental and delayed to represent this

3.4. PRESSURE DROP IN TRANSIENT FLOW 27

feature.

• Pressure loss mainly depends on two factors: the mass flow of the channel and its void content.

• Our artificial feedback loop should be able to maintain a steady state, which means that it should give a zero output if the system reaches the converged steady-state values.

Each one of these requirements need to be translated to the language of math- ematics and control systems, and the transfer function of the feedback loop has to reflect these properties. Fortunately, there are many different base types of con- trollers, which can be combined to do just that (Korondi 2016). Since mass flow distribution is essentially a result of the pressure loss differences between the chan- nels, the feedback circuit was designed to control the pressure loss of a channel depending on the flow properties. The detailed block diagram is shown on Fig. 3.3.

Figure 3.3. Block diagram of the implemented transient pressure loss feedback loop.

Looking at the diagram it becomes apparent, that this analogue representa- tion requires the introduction of several new constants, which influence the relative strength of the different signals or change the delay for example. These have been tuned through vigorous testing, to find a combination of values which ensures that

each component of the systems has a significant contribution to the total effect, and the system is not dominated by only one of the incoming signals. The mass flow and void fraction signals are subtracted from their equilibrium values, which means that a non-zero signal is only generated when the values are different than the equilibrium.

Chapter 4

Results and Discussion

4.1 Simulation Objectives and Parameters

This thesis project stemmed off of an idea to improve the educational materials available for the Reactor Dynamics and Stability course at KTH. To provide a tool for visualization of specific phenomena demonstrate the effect of various parameters on these phenomena.

Multiple series of simulations were carried out in order to visualize and demon- strate the instability events described in Chapter 2.4, and then to identify the factors which influence how these events unfold. Oscillations usually develop due to some kind of disturbance in the system and depending on the DR, the oscillation can be only a temporary short-term phenomenon, or can persist over a longer period. This disturbances are slight fluctuations in mass flow rates: In case of global oscillation the mass flow rate drops slightly all over the core while in case of a regional in- stability, half the core experiences a slight decrease while the mass flow rate in the other half is increased.

To make sure that these disturbances have observable consequences, their mag- nitude was increased, and two kinds of reference disturbances were designed, which could each be introduced to the system by pressing their specific buttons in the simulation program. This also ensured that the disturbances are truly identical in each and every simulation, and the changes in results are only due to the differences in other parameters, not differences in the disturbances themselves. These reference disturbances were:

1) A 30% overall decrease in mass flow in order to induce a global oscillation.

2) A momentary 50% increase in mass flow in the right half of the reactor, and a simultaneous 50% decrease of mass flow in the left half to induce regional instability.

29

First a reference point was established at nominal power and nominal mass flow rate, where the system’s response to the reference disturbances was recorded to pro- vide a basis for future comparisons. The program recorded flux, enthalpy and void fraction values for each time-step, and these were plotted out into sets of images.

Animations were also produced from these image sets to provide a better visual experience and easier interpretation.

The reactor flux changes during the transients are presented on Figures 4.1.

From the sequences of images, it is clear that the reactor is in a stable operational point at nominal power and mass flow, since the induced oscillations are strongly damped, and after several periods, their effects become negligible as the reactor returns to a steady state again.

Figure 4.1. Neutron flux response to the reference disturbances: disturbance 1) induces a global oscillation (top), while disturbance 2) induces a regional oscillation (bottom, the two lines represent the left and right sides of the reactor). An animated version of the oscillations is also available on these links: GlobalRegional

4.2. THE EFFECT OF INCREASED REACTOR POWER 31

4.2 The Effect of Increased Reactor Power

After the reference response was established, the main objective of the simulations was to gather data concerning the stability of the model reactor. The influence of two parameters was tested, first the effect of increased power and the second was the effect of reduced mass flow. Both changes increase the power-to mass flow ratio, and as we saw on Fig. 2.7 in Section 2.4, supposedly the higher this ratio is, the higher the decay ratio becomes, and consequently the reactor drifts toward insta- bility. Fig. 4.2 presents how increased reactor power affects the DR of a regional oscillation:

Figure 4.2. The flux values on both sides of the reactor plotted versus time during regional oscillation. From the top left to bottom right, the different charts show the system response to the reference regional instability initiating event at increased power levels.

Theoretically, increasing the power is equivalent to going upwards gradually on the operational diagram (Fig. 2.7) along a straight vertical line. Therefore our expectation is to see the DR increase, which means that the oscillation takes longer and longer to settle, until we reach a point where there is no visible dampening, and beyond which the oscillation becomes unstable and amplitudes start to increase instead.

This is exactly what we can see on Fig. 4.2, it is already visible at 110% power that the oscillation takes longer to settle, but as we reach 112.5% it is just barely

damped, and at 115% reactor power, the reactor is on the verge of becoming un- stable. As it was mentioned, several datasets were plotted out and animated for even better visualization. Table 4.1 contains link to the animations produced from the results of this first series of simulations. Images of several selected cases (re- gional and global oscillations, stable and unstable states) are also available in the Appendix (Section 6).

Reactor Power Flux Enthalpy Void

100% Link Link Link

110% Link Link Link

112.5% Link Link Link

115% Link Link Link

Table 4.1. Links to the animated .gif versions of the oscillations

4.3 The Effect of Reduced Mass Flow Rate

The second set of simulations was carried out the see how the model behaves with lower initial mass flow, while the power was kept on the nominal level. Fig 4.3 presents the flux plotted versus time with different (lower than nominal) mass flow rates.

As expected, the decay ratio becomes greater and greater with every step.A slight decrease of mass flow does not threaten stability, but gradually a point is reached at 70%, where the oscillation becomes sustained. Going further to 65%

we can observe an unstable state where instead of dampening, the oscillation is strengthened after every period.

Reactor Power Flux Enthalpy Void

80% Link Link Link

75% Link Link Link

70% Link Link Link

65% Link Link Link

Table 4.2. Links to the animations created from the results of the second round of simulations

4.4. STABILITY MAP 33

Figure 4.3. The flux values on both sides of the reactor plotted versus time. The diagrams show the model response to the reference disturbance with 80, 75, 70, and 65% of the nominal mass flow (From the top left to the bottom right) respectively.

These two round of simulations also included recording the system’s response to an induced global oscillation. The results were quite similar: both the trends, and the points where stability is lost matched perfectly. From this fact it can be concluded that the power-to-mass flow ratio affect the stability of both kinds of oscillations in the same way, even though their mechanisms are quite different.

4.4 Stability Map

The next objective was to create a stability map for the reactor, where it would be shown whether an operational point is table or not based on the mass flow and power value. Therefore several additional round of simulations were carried out with various mass flow and power values. To make the decay ratio calculation simpler

and more accurate, a global oscillation was induced in each case.

The decay ratios were calculated for each operational point and a surface was fit to the 16 recorded data points. Study of the data points led to the observations that along the power axis, the best fit is a quadratic polynomial, while along the mass flow axis, an exponential function yielded the best fit. Therefore the following structure was used when S(P, W ) surface function was fitted to the points:

S(P, W ) = a(bP + c)²e^{−f W +d}+ g (4.1) where P is the reactor power, W is the mass flow, while a, b, c, d, f and g are con- stants. Since the decay ratio is calculated by dividing two positive numbers, it always stays positive. This fact helped introducing some constraints for the con- stants, and a good fit was achieved. The resulting stability map for our model reactor is shown on Fig 4.4.

Compared to the operational diagram on Figure 2.7, this map looks very sim- ilar. The most stable operational points are those with the lowest power with the highest mass flow, while the reactor becomes unstable in those points where the power is high, while the mass flow rate is relatively low. Therefore we can con- clude that our model was able to reproduce the characteristics of a BWR, as the changes in mass flow rate and reactor power affected its stability in a similar fashion.

To sum it up, the results are generally in good agreement with our expectation and confirm the theory behind. This implies that they are valid and can be used to demonstrate the simulated stability related phenomena for educational purposes.

Several aspects of reactor stability remain untouched however, so there are several ways the project could be still extended upon.

For example, literature suggests that the stronger the void effect in the reactor, the more prone it becomes for global oscillation, while regional oscillations become heavily damped, and can switch from out-of phase to in-phase oscillations (B.- G. Bergdahl 2005). Furthermore, certain instability types (density wave and dual oscillations) were not explored and could be included in an extended model.